これは ZOZO Advent Calendar 2022 カレンダー Vol.4の15日目の記事です。昨日の投稿は @saitoshi さんの 「openapi-generatorをバージョンアップしたらimportMappingsが動かなくなった」でした。
概要
数値計算の誤差によって負の分散が発生することが知られています[1]。そこで、本記事では、分散の説明をした後に Python の numpy を使って負の分散が発生する具体例を紹介します。最後に安易に独自実装をしない方が良いことについても触れます。
分散の定義と式変形
データサイエンスでよく利用される分散は数値のばらつきを表す指標で、以下の数式で定義されます。
v = {\frac{1}{n}} {\sum_{i=1}^n}(x_i - {\bar{x}})^2 \tag{1}
$x_i$ は i 番目のデータの数値、 ${\bar{x}}$ は$x_i$ の平均値 ${\bar{x}} = {\frac{1}{n}} {\sum_{i=1}^n}x_i$ です。一方で、この式(1)を式変形して以下の式(2)を得ることができます。
v = {\frac{1}{n}} {\sum_{i=1}^n}x_i^2 - {\bar{x}}^2 \tag{2}
数式上では式(1), (2)は完全に等価ですが、数値計算では異なった値を示すことがあります。本記事では、その極端な例として式(2)の数値計算において負の分散が現れることを示します。何故、負の分散が極端な例かというと、式(1)の定義から分散は自乗の和なので常に正であることがわかり、数式上で等価な式(1)、(2)はともに正の値を示すはずだからです。
負の分散が現れる numpy を使った具体例
上記の分散の式(1)、(2)を numpy を使って実装し、結果がどうなるか実験します。
データ生成
以下のコードで分散を計算するための数値データを正規分布から乱数サンプリングします。
mu_true = 10**9 # 真の平均値
v_true = 3 # 真の分散
sample_size = 10**5 # 生成するサンプル数
# 正規分布から一様乱数でサンプルリング
x = numpy.random.normal(
loc=mu_true,
scale=v_true**0.5, # 真の標準偏差
size=sample_size,
)
ちなみに、正規分布の真の平均 $\mu_{\rm{true}}$ を$10^9$、真の分散 $v_{\rm{true}}$ を3としたのはわざと数値誤差を発生しやすくするためで後で説明します。
式(1)での数値計算
式(1)を numpy で実装すると以下になります。
v_cal1 = numpy.sum((x - numpy.mean(x)) ** 2) / sample_size
この計算結果の分散 $v_{\rm{cal1}}$ は 2.9992106767486497 となり、真の分散 $v_{\rm{true}}=3$ に近い値を得ました。ただし、今回数値データを乱数で生成しているので計算毎に数値が僅かに変わることがあります。
式(2)での数値計算
式(2)を numpy で実装すると以下になります。
v_cal2 = numpy.sum(x**2) / sample_size - numpy.mean(x) ** 2
この計算結果の分散 $v_{\rm{cal2}}$ は -384.0 となり負の分散が現れました。当然この結果は真の分散$v_{\rm{true}}$、計算結果の分散$v_{\rm{cal1}}$ともに大きく異なる値です。ただし、今回数値データを乱数で生成しており、この値は数値誤差の影響が大きいために実行毎に値は大きく変わります。
負の分散が現れる理由
今回、正規分布からサンプリングしたデータは大きな値$10^9$オーダーの数値に小さな値$10^1$以下のオーダーの数値が足されたものになっています。例えば、以下のような数値です。
x_0=1000000001.2040771=10^9+1.2040771
一方で、今回の数値計算では numpy のデータ型 float64（倍精度浮動小数点数型）が使われていて、上記のような数値を10進数の有効桁数15〜17桁[2]まで表現できます。
ここで話しを変えて、式(1)、(2)の自乗を行うタイミングに注目します。
式(1)では、 $x_i - {\bar{x}}$ の差分を取ってから自乗を行なっています。つまり、 $x_i - {\bar{x}}$ によって大きな値$10^9$オーダーの桁は打ち消されて小さな値$10^1$以下のオーダーの数値（例えば、式(3)）に対して自乗を行なうことになります。そして、この値の和を取ることで分散を算出します。
x_0 - {\bar{x}} = 10^9+1.2040771 - 999999999.9920104 = 1.2120667695999146 \tag{3}
式(2)では、 $x_i, {\bar{x}}$ に対してそのまま自乗を行なっています。つまり、大きな値$10^9$オーダーの自乗である$10^{18}$オーダーの数値に小さな値$10^1$以下のオーダーの数値を足した数値（例えば、式(4)）を扱うことになります。そして、この値の和を取ることで分散を算出します。
x_0^2= (10^9+1.2040771)^2 = 1.0000000024081542 {\times} 10^{18} + 1.4498016627444097 \tag{4}
上記のことから、今回の数値計算ではfloat64 の10進数の有効桁数15〜17桁を超えた桁数の数値は扱えないために、式(2)で扱う式(4)のような数値の第2項の小さな値は桁落ちでゼロと見做されて消えます。この桁落ちの数値誤差に伴って、データによっては式(2)の ${\frac{1}{n}} {\sum_{i=1}^n}x_i^2 > {\bar{x}}^2$ の大小関係が破れ、負の分散が発生することになります。
安易に独自実装しない方が良い
実は numpy には分散を計算してくれる var というメソッドが存在し、以下のように簡単に実行することができます。
v_cal3 = numpy.var(x)
この結果の分散$v_{cal3}$は2.9992106767486497となり、真の分散$v_{\rm{true}}=3$と近い値を示していることがわかります。このことから分散を求めるという簡単な実装であっても安易に独自実装せずに既存ライブラリを利用した方が良いと言えます。
参考文献
[1] 伊理 正夫(著), 藤野 和建(著) 『数値計算の常識』 共立出版(1985)
[2] Oracle® Solaris Studio 12.4: 数値計算ガイド