当記事は「基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)」の読解サポートにあたってChapter.$5$の「適合度検定」の章末問題の解説について行います。
基本的には書籍の購入者向けの解答例・解説なので、まだ入手されていない方は下記より入手をご検討ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。
Contents
章末の演習問題について
問題5.1の解答例
それぞれの分類の型に対応する三項分布の係数$\displaystyle \frac{9!}{a!b!c!}$と、型を与える観測値の個数$n_{(a,b,c)}$を求め、それぞれの積の和が$3261$に一致することを示せば良い。また$\chi^2$値は$\chi_{(a,b,c)}^2$のように表記することとする。以下それぞれの分類の型に対し、計算を行う。
・$(8,1,0)$
$$
\begin{align}
\frac{9!}{8!1!0!} &= 9 \\
n_{8,1,0} &= 3 \\
\chi_{(8,1,0)} &= \frac{(8-3)^2}{3} + \frac{(1-3)^2}{3} + \frac{(0-3)^2}{3} \\
&= \frac{38}{3}
\end{align}
$$
・$(7,2,0)$
$$
\begin{align}
\frac{9!}{7!2!0!} &= \frac{9 \cdot 8}{2} \\
&= 36 \\
n_{8,1,0} &= 3! \\
&= 6 \\
\chi_{(7,2,0)} &= \frac{(7-3)^2}{3} + \frac{(2-3)^2}{3} + \frac{(0-3)^2}{3} \\
&= \frac{26}{3}
\end{align}
$$
・$(7,1,1)$
$$
\begin{align}
\frac{9!}{7!1!1!} &= \frac{9 \cdot 8}{1} \\
&= 72 \\
n_{7,1,1} &= 3 \\
\chi_{(7,1,1)} &= \frac{(7-3)^2}{3} + \frac{(1-3)^2}{3} + \frac{(1-3)^2}{3} \\
&= \frac{24}{3} \\
&= 8
\end{align}
$$
・$(6,3,0)$
$$
\begin{align}
\frac{9!}{6!3!0!} &= \frac{9 \cdot 8 \cdot 7}{3 \cdot 2} \\
&= 3 \cdot 4 \cdot 7 \\
&= 84 \\
n_{6,3,0} &= 6 \\
\chi_{(6,3,0)} &= \frac{(6-3)^2}{3} + \frac{(3-3)^2}{3} + \frac{(0-3)^2}{3} \\
&= \frac{18}{3} \\
&= 6
\end{align}
$$
・$(6,2,1)$
$$
\begin{align}
\frac{9!}{6!2!1!} &= \frac{9 \cdot 8 \cdot 7}{2} \\
&= 9 \cdot 4 \cdot 7 \\
&= 252 \\
n_{6,3,0} &= 6 \\
\chi_{(6,2,1)} &= \frac{(6-3)^2}{3} + \frac{(2-3)^2}{3} + \frac{(1-3)^2}{3} \\
&= \frac{14}{3}
\end{align}
$$
・$(5,4,0)$
$$
\begin{align}
\frac{9!}{5!4!0!} &= \frac{9 \cdot 8 \cdot 7 \cdot 6}{4 \cdot 3 \cdot 2} \\
&= 9 \cdot 2 \cdot 7 \\
&= 126 \\
n_{5,4,0} &= 6 \\
\chi_{(6,2,1)} &= \frac{(5-3)^2}{3} + \frac{(4-3)^2}{3} + \frac{(0-3)^2}{3} \\
&= \frac{14}{3}
\end{align}
$$
上記より、$\displaystyle \chi^2 \geq \frac{14}{3}$となる確率の分子は下記のように計算できる。
$$
\begin{align}
1 \cdot 3 &+ 9 \cdot 6 + 36 \cdot 6 + 72 \cdot 3 + 84 \cdot 6 \\
&+ 252 \cdot 6 + 126 \cdot 6 = 3261
\end{align}
$$
上記より、$\displaystyle \chi^2 \geq \frac{14}{3}$となる確率は、$\displaystyle \chi^2 \geq \frac{3261}{3^9}=0.166$となる。
問題5.2の解答例
$\chi^2$値は下記の計算を実行することで得ることができる。
import numpy as np
observed = np.array([60., 208., 655., 1826., 3650., 5909., 7512., 7737., 6063., 3641., 1820., 683., 188., 48.])
prob = np.array([0.00135, 0.00486, 0.01654, 0.04406, 0.09185, 0.14988, 0.19146, 0.19146, 0.14988, 0.09185, 0.04406, 0.01654, 0.00486, 0.00135])
estimate = prob*40000
print(estimate)
print((observed-estimate)**2/estimate)
print(np.sum((observed-estimate)**2/estimate))・実行結果
> print(estimate)
[ 54. 194.4 661.6 1762.4 3674. 5995.2 7658.4 7658.4 5995.2
3674. 1762.4 661.6 194.4 54. ]
> print((observed-estimate)**2/estimate)
[ 0.66666667 0.95144033 0.06584039 2.29514299 0.15677735 1.23939819
2.79862112 0.80669069 0.7667534 0.29640719 1.88252383 0.69220073
0.21069959 0.66666667]
> print(np.sum((observed-estimate)**2/estimate))
13.4958291239上記より$\chi^2 \simeq 13.5$であり、$\chi^2 < 22.36… = \chi^2_{\alpha=0.05}(13)$より帰無仮説は棄却できない。
問題5.3の解答例
$\chi^2$分布において上側確率が$100\alpha$％となるパーセント点に対応する$\chi^2$の値を$\chi^2_{\alpha}$とする。(5.13)式を元に$\chi^2$の値の計算を行う。
$$
\begin{align}
\chi^2 = \sum_{i}^{a} \sum_{j}^{b} \frac{(f_{ij}-f_{i \cdot}f_{\cdot j}/n)^2}{f_{i \cdot}f_{\cdot j}/n}
\end{align}
$$
以下、表の数字を元に上記を用いて計算する。
$$
\begin{align}
\chi^2 &= \sum_{i}^{a} \sum_{j}^{b} \frac{(f_{ij}-f_{i \cdot}f_{\cdot j}/n)^2}{f_{i \cdot}f_{\cdot j}/n} \\
&= \frac{(7142-17884 \cdot 7232/18101)^2}{17884 \cdot 7232/18101} + \frac{(3021-17884 \cdot 3081/18101)^2}{17884 \cdot 3081/18101} + \frac{(1841-17884 \cdot 1879/18101)^2}{17884 \cdot 1879/18101} \\
&+ \frac{(5880-17884 \cdot 5909/18101)^2}{17884 \cdot 5909/18101} + \frac{(90-217 \cdot 7232/18101)^2}{217 \cdot 7232/18101} + \frac{(60-217 \cdot 3081/18101)^2}{217 \cdot 3081/18101} \\
&+ \frac{(38-217 \cdot 1879/18101)^2}{217 \cdot 1879/18101} + \frac{(29-217 \cdot 5909/18101)^2}{217 \cdot 5909/18101} \\
&= 50.4733…
\end{align}
$$
上記が自由度$(2-1)(4-1)=3$の$\chi^2$分布に従うため、有意水準5％で片側検定するにあたっては$\chi^2_{\alpha=0.05}(3)=7.815$と比較すればよい。
このとき、$\chi^2=50.4733…>7.815=\chi^2_{\alpha=0.05}(3)$のため、独立を仮定した帰無仮説は棄却される。（独立性が成立しないと考える方が妥当である）
問題5.4の解答例
i)
試合数$z$の確率分布$p(z)$は$z-1$回目までに勝利チームが$3$回勝利する確率と考えることができる。ここで「勝利チームが$z-1$回目までに$3$回勝利する確率」としているので、$z$回目の分岐は生じない。（両チームの勝率が同じため分岐が生じると考えて、チームの重複を考えるという方法でも良い）
$p(z)$は下記のように計算できる。
$$
\begin{align}
p(z) &= {}_{z-1} C_{3} \left( \frac{1}{2} \right)^3 \left( 1 – \frac{1}{2} \right)^{z-1-3} \\
&= {}_{z-1} C_{3} \left( \frac{1}{2} \right)^3 \left( \frac{1}{2} \right)^{z-4} \\
&= {}_{z-1} C_{3} \left( \frac{1}{2} \right)^{z-1} \\
&= {}_{z-1} C_{3} 2^{-z+1}
\end{align}
$$
ⅱ)
$\chi^2$分布において上側確率が$100\alpha$％となるパーセント点に対応する$\chi^2$の値を$\chi^2_{\alpha}$とする。
$$
\begin{align}
\chi^2 = \sum \frac{(O-E)^2}{E}
\end{align}
$$
観測度数を$O$、理論度数を$E$とした際に、上記で得られる$\chi^2$の値を元に$\chi^2$検定を行えばよい。理論度数は「試行回数×理論確率」で計算できるので、この問題において$\chi^2$は下記のように計算できる。
$$
\begin{align}
\chi^2 &= \sum \frac{(O-E)^2}{E} \\
&= \frac{(5-42p(z=4))^2}{42p(z=4)} + \frac{(8-42p(z=5))^2}{42p(z=5)} + \frac{(15-42p(z=6))^2}{42p(z=6)} + \frac{(14-42p(z=7))^2}{42p(z=7)} \\
&= \frac{(5-42{}_{4-1} C_{3} 2^{-4+1})^2}{42{}_{4-1} C_{3} 2^{-4+1}} + \frac{(8-42{}_{5-1} C_{3} 2^{-5+1})^2}{42{}_{5-1} C_{3} 2^{-5+1}} + \frac{(15-42{}_{6-1} C_{3} 2^{-6+1})^2}{42{}_{6-1} C_{3} 2^{-6+1}} + \frac{(14-42{}_{7-1} C_{3} 2^{-7+1})^2}{42{}_{7-1} C_{3} 2^{-7+1}} \\
&= \frac{(5-42{}_{3} C_{3} 2^{-3})^2}{42{}_{3} C_{3} 2^{-3}} + \frac{(8-42{}_{4} C_{3} 2^{-4})^2}{42{}_{4} C_{3} 2^{-4}} + \frac{(15-42{}_{5} C_{3} 2^{-5})^2}{42{}_{5} C_{3} 2^{-5}} + \frac{(14-42{}_{6} C_{3} 2^{-6})^2}{42{}_{6} C_{3} 2^{-6}} \\
&= \frac{(5-42 \cdot 2^{-3})^2}{42 \cdot 2^{-3}} + \frac{(8-42 \cdot 4 \cdot 2^{-4})^2}{42 \cdot 4 \cdot 2^{-4}} + \frac{(15-42 \cdot 10 \cdot 2^{-5})^2}{42 \cdot 10 \cdot 2^{-5}} + \frac{(14-42 \cdot 20 \cdot 2^{-6})^2}{42 \cdot 20 \cdot 2^{-6}} \\
&= 0.93333…
\end{align}
$$
上記が自由度$4-1=3$の$\chi^2$分布に従うため、有意水準5％で片側検定するにあたっては$\chi^2_{\alpha=0.05}(3)=7.815$と比較すればよい。
このとき、$\chi^2=0.93333…<7.815=\chi^2_{\alpha=0.05}(3)$のため、帰無仮説は棄却されない。（得られた結果は妥当と考える方が良い）
問題5.5の解答例
i)
平均$\bar{x}$、不偏分散$s^2$はそれぞれ下記のように求めることができる。
$$
\begin{align}
\bar{x} &= \frac{1}{228}(1 \cdot 24 + 2 \cdot 16 + 3 \cdot 8 + 4 \cdot 3 + 5 \cdot 2) \\
&= 0.44736… \\
s^2 &= \frac{1}{227} \left( 175(0-\bar{x})^2 + 24(1-\bar{x})^2 + 16(2-\bar{x})^2 + 8(3-\bar{x})^2 + 3(4-\bar{x})^2 + 2(5-\bar{x})^2 \right) \\
&= 0.93554…
\end{align}
$$
問題5.6の解答例
下記を実行することで$\theta$の推定を行うことができる。
import numpy as np
K = 100.
y = np.array([[0., 15., 15., 15., 13., 16.], [11., 0., 14., 15., 14., 17.], [11., 12., 0., 14., 17., 13.], [11., 11., 12., 0., 13., 19], [13., 12., 9., 13., 0., 17.], [10., 9., 13., 7., 9., 0.]])
theta = np.repeat(100., 6)/6.
for epoch in range(10):
for i in range(6):
r_theta = 0
for j in range(6):
if i != j:
r_theta += (y[i,j]+y[j,i])/(theta[i]+theta[j])
theta[i] = np.sum(y[i,:])/r_theta
theta = K*theta/np.sum(theta)
print(theta)・実行結果
[ 20.63301132 19.08005757 17.19619821 16.75537645 15.90772652
10.42762994]また、$\chi^2$検定は上記で計算を行なったthetaを用いて、下記を実行することで行うことができる。
from scipy import stats
chi2 = 0
expected_y = np.zeros([6,6])
for i in range(6):
for j in range(6):
if i != j:
expected_y[i,j] = theta[i]*(y[i,j]+y[j,i])/(theta[i]+theta[j])
chi2 += (y[i,j]-expected_y[i,j])**2 / expected_y[i,j]
print("・expected_y")
print(expected_y)
print("・chi^2 test of goodness of fit")
if chi2 > stats.chi2.ppf(1.-0.05,10):
print("chi^2: {:.2f}, P_value: {:.2f}, reject H_0 and expected_y is not good.".format(chi2,stats.chi2.cdf(chi2,10)))
else:
print("chi^2: {:.2f}, P_value: {:.2f}, accept H_0 and expected_y seems to be good.".format(chi2,stats.chi2.cdf(chi2,10)))・実行結果
・expected_y
[[ 0. 13.50835655 14.18106011 14.34825962 14.68110186
17.27132064]
[ 12.49164345 0. 13.67510197 13.84332325 14.17870579
16.81194085]
[ 11.81893989 12.32489803 0. 13.1687899 13.50598629
16.18534372]
[ 11.65174038 12.15667675 12.8312101 0. 13.33736688
16.02618127]
[ 11.31889814 11.82129421 12.49401371 12.66263312 0. 15.70515631]
[ 8.72867936 9.18805915 9.81465628 9.97381873 10.29484369 0. ]]
・chi^2 test of goodness of fit
chi^2: 6.84, P_value: 0.26, accept H_0 and expected_y seems to be good.問題5.7の解答例
$\chi^2$分布において上側確率が$100\alpha$％となるパーセント点に対応する$\chi^2$の値を$\chi^2_{\alpha}$とする。
一様分布を想定した際の$10^i$桁までに対応する$\chi^2$適合度統計量を$\chi_i^2$とおき、$\chi_4^2$から$\chi_6^2$までを計算する。（$\chi_7^2$〜$\chi_9^2$の値は書籍に記載があり、途中式がわかれば十分と思われるので、$\chi_7^2$〜$\chi_9^2$は取り扱わないものとする。）
$$
\begin{align}
\chi_4^2 &= \frac{(968-1000)^2}{1000} + \frac{(1026-1000)^2}{1000} + \frac{(1021-1000)^2}{1000} + \frac{(974-1000)^2}{1000} + \frac{(1012-1000)^2}{1000} \\
&+ \frac{(1046-1000)^2}{1000} + \frac{(1021-1000)^2}{1000} + \frac{(970-1000)^2}{1000} + \frac{(948-1000)^2}{1000} + \frac{(1014-1000)^2}{1000} \\
&= 9.318 \\
\chi_5^2 &= \frac{(9999-10000)^2}{10000} + \frac{(10137-10000)^2}{10000} + \frac{(9908-10000)^2}{10000} + \frac{(10025-10000)^2}{10000} + \frac{(9971-10000)^2}{10000} \\
&+ \frac{(10026-10000)^2}{10000} + \frac{(10029-10000)^2}{10000} + \frac{(10025-10000)^2}{10000} + \frac{(9978-10000)^2}{10000} + \frac{(9902-10000)^2}{10000} \\
&= 4.093 \\
\chi_6^2 &= \frac{(99959-100000)^2}{10000} + \frac{(99758-100000)^2}{10000} + \frac{(100026-100000)^2}{10000} + \frac{(100229-100000)^2}{10000} + \frac{(100230-100000)^2}{10000} \\
&+ \frac{(99548-100000)^2}{10000} + \frac{(100359-100000)^2}{10000} + \frac{(99800-100000)^2}{10000} + \frac{(99985-100000)^2}{10000} + \frac{(100106-100000)^2}{10000} \\
&= 5.50908
\end{align}
$$
上記はそれぞれ自由度$10-1=9$の$\chi^2$分布に従うため、有意水準5％で片側検定するにあたっては$\chi^2_{\alpha=0.05}(9)=16.919$と比較すればよい。
$\chi_4^2<\chi^2_{\alpha=0.05}(9)$、$\chi_5^2<\chi^2_{\alpha=0.05}(9)$、$\chi_6^2<\chi^2_{\alpha=0.05}(9)$より、等確率(一様分布)で分布すると下帰無仮説は棄却されない。（一様分布にしたがっていると考える方が妥当である）
また、書籍より$\chi_7^2$〜$\chi_9^2$についても同様であることが確認できる。
まとめ
Chapter.$5$の「適合度検定」の演習問題について確認を行いました。様々な問題のパターンはある一方で、基本的には$\displaystyle \chi^2 = \sum \frac{(O-E)^2}{E}$を用いて$\chi^2$検定を行うだけではあるので、解法の整理はしやすいと思います。
[…] ここで得たXに対して問題5.2と同じ要領で適合度検定を行うことを考える。 […]