円周率1,000,000桁で各数字が等確率で出現するか問題

math
  • 0から9の各数字が、円周率1000桁に等確率で出現するかを仮説検定してみる

一般に、

ある属性Aによって、n個の個体がk種のカテゴリー A_1, A_2, ..., A_k へ分類されるとき
各カテゴリーへ属する観測度数が f_1, f_2, ..., f_k であるとする
これが、各カテゴリーの理論確率 p_1, p_2, ..., p_k に適合するかを見るには、これが正しければ生じるであろ
う理論度数 np_1, np_2, ..., np_k を、観測度数と比べ、K.ピアソンの適合度基準
\chi^2 = \Sigma^k_{i=1} \frac{(f_i-np_i)^2}{np_i}
で判断すればよい。
この適合度のχ^2統計量χ^2はnが大きいとき、自由度 k-1 の χ^2 分布 \chi^2(k-1) に従う。

帰無仮説を
H_0: P(A_1)=p_1, ..., P(A_k)=p_k
とするとき、
\chi^2 > \chi^2_\alpha(k-1)
ならば、"観測度数は理論確率分布 p_1, ..., p_k に適合している" という仮説H_0は、有意水準\alphaで棄却される。ここで、自由度=カテゴリー数(k)-1である。

円周率1,000桁でやってみる

計算

# observed frequency of 0, 1, ..., 9 (これは適当にスクリプトでカウントした)
freq <- c(93,116,103,102,93,97,94,95,101,106)

# expected probability
p <- 0.1

# number of observations
n <- 1000 

np <- n * p

chi <- sum((freq - np)^2 / np)
# => 4.74 

qchisq(0.95, 9)
# => 16.91898

chiの値(4.74) < \chi^2(9)(16.91898) なので、有意水準5%で棄却されない(0, ..., 9 は等確率で出現している)。

1,000,000桁でもやってみる

# これは適当にスクリプトでカウントした
freq <- c(99959,99758,100026,100229,100230,100359,99548,99800,99985,100106)
p <- 0.1
n <- 10^6
np <- n * p

chi <- sum((freq - np)^2 / np)
# => 5.50908

qchisq(0.95, 9)
# => 16.91898

\chi^2(5.50908) < \chi^2_{0.05}(9)(16.91898) なので、棄却されない。

同じ事だけどp値の計算でやると

pchisq(5.50908, 9, lower.tail=F)
# => 0.7878669

p値(0.7878669) > 0.05 (有意水準) となるので帰無仮説は棄却されない。

参考文献