Rで cut()による数値のカテゴリー化の時、 seq()を使うと計算誤差が問題になる

データ解析をするとき、連続変数を一定の値範囲ごとにカテゴリー化することがたまにありますが、その時には cut() 関数を使うことが多いと思います。 ただ、cutのbreaks引数に、区切りが一定だからといってseq()を使う、かつその時にseqの増分(by引数) に小数を使うと、思わぬトラブルに遭うことがある、という小ネタです。

たとえばこうすると、irisのSepal.Length は, 2.5, 2.6, 2.7,…を区切りとして、right=FALSEですから左側が閉じた区間となるようにカテゴリー化されるはずです。

cut(iris$Sepal.Length, breaks=seq(2.5, 6.5, by=0.1), right=FALSE)

ところで、iris$Sepal.Lengthには、値が4.8である個体が5例含まれています。

> sum(iris$Sepal.Length == 4.8)
[1] 5

iris$Sepal.Length には小数点以下第2位以下まであるデータはないように見えますから、[4.8,4.9) も5例あるはずですね。

しかし、実際にカウントすると…

> sum(cut(iris$Sepal.Length, breaks=seq(2.5, 6.5, by=0.1), right=FALSE)=="[4.8,4.9)",na.rm=TRUE)
[1] 0

0例となってしまいました。

これは、cut関数やseq関数のバグ、というわけではなく、小数のコンピューター内での表現の精度の問題です。 こうすると、可視化できます。

> print(seq(2.5, 6.5, by=0.1))
 [1] 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9
[36] 6.0 6.1 6.2 6.3 6.4 6.5
> print(seq(2.5, 6.5, by=0.1), digits=10)
 [1] 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9
[36] 6.0 6.1 6.2 6.3 6.4 6.5
> print(seq(2.5, 6.5, by=0.1), digits=20)
 [1] 2.5000000000000000000 2.6000000000000000888 2.7000000000000001776 2.7999999999999998224 2.8999999999999999112 3.0000000000000000000
 [7] 3.1000000000000000888 3.2000000000000001776 3.2999999999999998224 3.3999999999999999112 3.5000000000000000000 3.6000000000000000888
[13] 3.7000000000000001776 3.7999999999999998224 3.9000000000000003553 4.0000000000000000000 4.0999999999999996447 4.2000000000000001776
[19] 4.2999999999999998224 4.4000000000000003553 4.5000000000000000000 4.5999999999999996447 4.7000000000000001776 4.8000000000000007105
[25] 4.9000000000000003553 5.0000000000000000000 5.0999999999999996447 5.2000000000000001776 5.3000000000000007105 5.4000000000000003553
[31] 5.5000000000000000000 5.5999999999999996447 5.7000000000000001776 5.8000000000000007105 5.9000000000000003553 6.0000000000000000000
[37] 6.0999999999999996447 6.2000000000000001776 6.3000000000000007105 6.4000000000000003553 6.5000000000000000000

digitsはminimum number of significant digits to be printed in values を示す引数とされているので、10桁までの表示であればseqが出力する4.8 は本当に4.8なのですが、20桁まで表示すると実は 4.8000000000000007105 となっています。

一方、

> print(4.8, digits=20)
[1] 4.7999999999999998

ですから、4.8とただ書いただけでも厳密には4.8としては格納されていないのですが、比較自体も完全に厳密に実施されるわけではないので、通常問題にはなりません。しかし…

> 4.7 == seq(2.5, 6.5, 0.1)[23]
[1] TRUE
> 4.8 == seq(2.5, 6.5, 0.1)[24]
[1] FALSE
> 4.9 == seq(2.5, 6.5, 0.1)[25]
[1] TRUE

というように、seq(2.5,6.5,0.1) で数列を作ったときには、小数の加算が繰り返されて誤差が蓄積するためか、4.7と4.9は問題ないのですが、4.8のところでこの誤差が比較で同一とみなされる範囲を超えてしまうようで、この数列の24番目の値と4.8が同一とみなされなくなってしまいます。

有効数字が20桁ぐらいまでいかないと顕在化しない問題なので、計算結果自体がこれでおかしくなることはあまりないでしょうが、プログラミングとして考えた時には予期しない結果を招くと思われますので、注意しましょう。対策として、 ***cut関数に対するbreaks引数としては、seq関数で生成した小数は使わない *** ことが必要かと思います。