Hatena::ブログ(Diary)

Haskell雑記帳 RSSフィード Twitter

2012-03-04

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

  • 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桁でやってみる

  • 観測数 n=1000
  • 帰無仮説 H_0: P(1)=P(2)=...=P(9)=0.1
  • 円周率ここにあるものを参照する

計算

# 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 (有意水準) となるので帰無仮説は棄却されない。

参考文献

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/epsilon_delta/20120304/1330832580
リンク元
Connection: close