分散がマイナスになる時

数式上は全く同じ結果になるはずが、計算機誤差のために全く違う結果になる例として、分散の計算をDave Gilesが示している。その内容は概ね以下の通り。


データの平均を
        X* = [Σ(Xi)] / n
とした場合、分散の式として最初に習うのは、以下の(1)式であろう。
        Var(X) = (1 / n) Σ [(Xi - X*)2]             (1)
そして、それと完全に等価な式として、以下の(2)式を習うであろう。
        Var(X) = (1 / n) [Σ (Xi2) - nX*2]            (2)

どちらの式を使うか選択する際に考慮すべき要因としてGilesは、計算速度と正確性を挙げている。そして、正確性が違ってくることを示すために、以下のRのコードを提示している。

> set.seed(100)
> n<- 10^3
> x<-  rnorm(n)
> xbar<- mean(x)
> x2<- sum(x^2)
> var1<- sum((x-xbar)^2)/n
> var2<- (x2 - n*xbar^2)/n
> var1
[1] 1.061049
> var2
[1] 1.061049

このコードの場合は、どちらの式を使っても結果に差が生じない。
だが、データ値に100万を加えると以下のようになる。

> set.seed(100)
> n<- 10^3
> x<-  rnorm(n)+10^6
> xbar<- mean(x)
> x2<- sum(x^2)
> var1<- sum((x-xbar)^2)/n
> var2<- (x2 - n*xbar^2)/n
> var1
[1] 1.061049
> var2
[1] 1.06125

データをシフトしただけなので、本来は分散の値は変わらないはずで、実際、1番目の式は前と同じ結果を返している。しかし2番目の式の結果は、似て非なるものとなっている。


データ値に加える値を1億にするとどうなるか?

> set.seed(100)
> n<- 10^3
> x<-  rnorm(n)+10^8
> xbar<- mean(x)
> x2<- sum(x^2)
> var1<- sum((x-xbar)^2)/n
> var2<- (x2 - n*xbar^2)/n
> var1
[1] 1.061049
> var2
[1] 0
> 

これもとんでもない結果だが、加える値を100億にすると…

> set.seed(100)
> n<- 10^3
> x<-  rnorm(n)+10^10
> xbar<- mean(x)
> x2<- sum(x^2)
> var1<- sum((x-xbar)^2)/n
> var2<- (x2 - n*xbar^2)/n
> var1
[1] 1.061049
> var2
[1] -16777.22
> 

何とマイナスの分散が出てきてしまう。


以上の結果についてGilesは、(1)式の2乗項は10のゼロ乗のオーダーの数値を2乗しているのに対し、(2)式の(Rコード中の変数名で言う)「x2」、および、「xbar」の2乗は、10の20乗のオーダーになっていることを指摘している。それだけ大きな数値の差を取れば、丸め誤差が発生する余地が大いに出てきてしまう、というわけだ。


なお、コメント欄では、Rのvar関数は平均計算にTwo-pass Algorithmsを使っているために上の1番目の式よりさらに正確な結果を返すこと、および、こちらのブログではさらに正確な結果を返す合計のアルゴリズムが提示されていること、が指摘されている。