計算精度

一般的な電卓の計算精度というか扱える数値の最大値は我々が思っているよりは小さい。一般的には 1e100 程度だ。

コンピュータはそれよりは大きいものの限度はある。.Machine によれば，

> .Machine
　　：
$double.xmax
[1] 1.797693e+308

ということなので，例えば，階乗は，170! までしか求められない。

> factorial(170)
[1] 7.257416e+306
> factorial(171)
[1] Inf
 警告メッセージ：
In factorial(171) :  'gammafn' 中の値が範囲を超えています

要するに，必要な計算をする場合に，途中でこの範囲（限界）を超えてはならないということだ。

例えば，2×2分割表の独立性を検定するときに，χ自乗統計量を計算 n*(ad-bc)/(a+b)/(c+d)/(a+c)/(b+d) するわけだが（記号法の説明は省略），その計算順序はどうであってもよいわけではない（普通はどうでもよいのだけど）。

> func1 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n*(a*d-b*c)^2/(ab*cd*ac*bd))
+ }

> func2 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n*(a*d-b*c)^2/ab/cd/ac/bd)
+ }

> func3 <- function(x) {
+     y <- addmargins(x)
+     a  <- y[1,1]; b  <- y[1,2]; ab <- y[1,3]
+     c  <- y[2,1]; d  <- y[2,2]; cd <- y[2,3]
+     ac <- y[3,1]; bd <- y[3,2]; n  <- y[3,3]
+     return(n/ab*(a*d-b*c)/cd/ac*(a*d-b*c)/bd)
+ }

普通の範囲ではどれでもちゃんと計算できる # 1-1, # 1-2, # 1-3

> x <- matrix(1:4, 2)
> func1(x) # 1-1
[1] 0.07936508
> func2(x) # 1-2
[1] 0.07936508
> func3(x) # 1-3
[1] 0.07936508

>しかし，例外的な場合には，計算順序を考えないと計算は失敗する # 2-1, # 2-2 は失敗する。しかし，# 2-3 はちゃんと計算できる。

x <- x*1e62
> func1(x) # 2-1
[1] Inf # 何てっこった
> func2(x) # 2-2
[1] Inf # 何てっこった
> func3(x) # 2-3
[1] 7.936508e+60 # 正しい

途中の計算を log を使って，一番最後に exp を使ってという方策で対処できる場合もあるけど，（アルゴリズムがマズイ場合は問題外だが），それでもダメ（途中の結果が限界を超える）ということもある。

「そのような場合は考慮外」とすましていることもできるけど，それは，やはり恥ずかしいことだよ。

どんな場合にも正しい答えを与えるようにプログラムする事は大変難しい（場合もある。プログラマのスキルにもよる）。

しかし，可能な限り，どんな場合にも正しい答えを出せるようにするのが「職人」だ。悔しければ，やってみな。