当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.12の「一般の分布に関する検定法」に関して演習問題を中心に解説を行います。尤度比検定に加え、二項分布やポアソン分布に関する検定や適合度検定などは抑えておくと良いように思われました。
本章のまとめ
演習問題解説
例12.1
・標準正規分布を用いた検定
最尤推定値$\displaystyle \hat{\theta} = \frac{12}{30} = 0.4$を用いて、検定統計量の値は下記のように計算できる。
$$
\large
\begin{align}
\frac{\sqrt{n}(\hat{\theta}-\theta_0)}{\sqrt{\theta_0(1-\theta_0)}} &= \frac{\sqrt{30}(0.4-0.5)}{\sqrt{0.5 \times (1-0.5)}} \\
&= -1.0954…
\end{align}
$$
上記の絶対値が標準正規分布の上側$2.5$％点の$1.96$を下回ることから、帰無仮説は棄却できない。
・尤度比検定
尤度比は下記のように計算できる。
$$
\large
\begin{align}
2n \left( \hat{\theta} \log{\frac{\hat{\theta}}{\theta_0}} + (1-\hat{\theta}) \log{\frac{1-\hat{\theta}}{1-\theta_0}} \right) &= 60 \times \left( 0.4 \log{\frac{0.4}{0.5}} + 0.6 \log{\frac{0.6}{0.5}} \right) \\
&= 1.2081…
\end{align}
$$
ここで自由度$1$の$\chi^2$分布の上側$5$％点がおよそ$3.84$であることから、帰無仮説が棄却できないことがわかる。
例12.2
母比率の検定は、適合度検定においてカテゴリ数を$2$と考える場合に対応する。このとき、帰無仮説の$H_0: p=\theta_0$に基づく期待度数はそれぞれ$n \theta_0$と$n (1-\theta_0)$のように表せる。このとき、$\chi^2$統計量は下記のように表せる。
$$
\large
\begin{align}
\chi^2 &= \frac{(x_1 – n \theta_0)^2}{n \theta_0} + \frac{(n-x_1 – n(1-\theta_0))^2}{n (1-\theta_0)} \\
&= \frac{(1-\theta_0)(x_1 – n \theta_0)^2 + \theta_0(n \theta_0 – x_1)}{n \theta_0 (1-\theta_0)} \\
&= \frac{(x_1 – n \theta_0)^2}{n \theta_0 (1-\theta_0)} \\
&= \frac{n^2(x_1/n – \theta_0)^2}{n \theta_0 (1-\theta_0)} \\
&= \frac{n(\hat{\theta} – \theta_0)^2}{\theta_0 (1-\theta_0)}
\end{align}
$$
上記が$(12.1)$式で用いた統計量の$2$乗に対応するので、母比率に関する両側検定は適合度検定の特殊ケースであることが確認できる。
例12.3
独立性の仮説は$p_{ij} = \alpha_{i}\beta_{j}$で表せる。また、自由度は$d=(I-1)(J-1)$である。
詳細の考え方はワークブックが詳しいのでここでは省略。
問12.1
$[1]$
下記を実行することで$\chi^2$統計量を計算することができる。
import numpy as np
observed = np.array([7., 3., 5., 2., 7., 12., 13.])
estimate = np.repeat(np.mean(observed), observed.shape[0])
np.sum((observed-estimate)**2)/np.mean(observed)・実行結果
> np.sum((observed-estimate)**2)/np.mean(observed)
15.142857142857142$[2]$
$\chi^2_{\alpha=0.05}(7-1)=\chi^2_{\alpha=0.05}(6)=12.59$より、帰無仮説は棄却され。問い合わせの数が曜日によって異なると考えることができる。
問12.2
$[1]$
ポアソン分布は平均と分散が等しい分布であるが、ここで標準偏差の二乗が$1.7^2=2.89$のように平均の$2.99$に近い値を示すことが確認できる。これより、観測値がポアソン分布に従う可能性が示唆されると読み取れる。
$[2]$
ポアソン分布の確率関数を用いることで、下記のような計算を行うことができる。
import math
import numpy as np
def factorial(vec):
res = np.zeros(vec.shape[0])
for i in range(vec.shape[0]):
res[i] = math.factorial(vec[i])
return res
lamb, count = 2.99, 69.
x = np.arange(0,11,1)
p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[10] = 1-np.sum(p_x[0:10])
print("prob: {}".format(p_x))
print("estimate: {}".format(p_x*count))・実行結果
> print("prob: {}".format(p_x))
prob: [ 0.05028744 0.15035944 0.22478736 0.22403807 0.16746845 0.10014614
0.04990616 0.02131706 0.00796725 0.0026469 0.00079142]
> print("estimate: {}".format(p_x*count))
estimate: [ 3.46983313 10.37480107 15.5103276 15.45862651 11.55532331
6.91008334 3.44352487 1.47087705 0.5497403 0.18263594
0.07422687]上記の結果を元に表を作成することができる。
$[3]$
下記の計算に基づいて$\chi^2$統計量が$16.37$であることを示すことができる。
import math
import numpy as np
def factorial(vec):
res = np.zeros(vec.shape[0])
for i in range(vec.shape[0]):
res[i] = math.factorial(vec[i])
return res
lamb, count = 2.99, 69.
x = np.arange(0,11,1)
observed = np.array([4., 7., 17., 18., 12., 7., 3., 0., 0., 0., 1.])
p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[10] = 1-np.sum(p_x[0:10])
estimate = p_x*count
chi2_1 = np.sum((observed-estimate)**2/estimate)
estimate[10] = 0.07
chi2_2 = np.sum((observed-estimate)**2/estimate)
print("chi2-statistic 1: {}".format(chi2_1))
print("chi2-statistic 2: {}".format(chi2_2))・実行結果
> print("chi2-statistic 1: {}".format(chi2_1))
chi2-statistic 1: 15.5647586322
> print("chi2-statistic 2: {}".format(chi2_2))
chi2-statistic 2: 16.3740364649上記のプログラムのメインの処理では$15.56$が得られたが、estimate[10]に$0.07$を代入すると$16.37$が得られることが確認できた。これは誤差の丸め方によって変わる。
また、ここで棄却域を考えるにあたっては自由度が$11-1-1=9$より、$\chi^2_{\alpha=0.05}(9)=16.92$を用いる。統計量を大きくしている理由に関しては、下記のように$10$回以上の期待度数が小さいことに起因する。
print((1.-0.07)**2/0.07)・実行結果
> print((1.-0.07)**2/0.07)
12.3557142857$[4]$
下記を実行することで。$\chi^2$統計量を計算できる。
import math
import numpy as np
def factorial(vec):
res = np.zeros(vec.shape[0])
for i in range(vec.shape[0]):
res[i] = math.factorial(vec[i])
return res
lamb, count = 2.99, 69.
x = np.arange(0,7,1)
observed = np.array([4., 7., 17., 18., 12., 7., 4.])
p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[6] = 1-np.sum(p_x[0:6])
estimate = p_x*count
chi2 = np.sum((observed-estimate)**2/estimate)
print("chi2-statistic: {}".format(chi2))・実行結果
> print("chi2-statistic: {}".format(chi2))
chi2-statistic: 2.27565946483上記は$\chi^2_{\alpha=0.10}(7-1-1)=9.24$よりも小さい。一方で$[3]$の結果の$16.37$は$\chi^2_{\alpha=0.10}(11-1-1)=14.68$よりも大きいので、$P$値で比較すると、上陸数をまとめたときの方があてはまりが良いと考えることができる。
問12.3
$$
\large
\begin{align}
\frac{\hat{\theta}_1-\hat{\theta}_2}{\sqrt{\left(\frac{1}{n_1}+\frac{1}{n_2}\right) \hat{\theta}_{*}(1-\hat{\theta}_{*})}} \geq z_{\alpha} \quad (12.4)
\end{align}
$$
上記で表した$(12.4)$式の左辺に対し、$n_1=114, n_2=107, \hat{\theta}_1=40/114, \hat{\theta}_2=62/107, \hat{\theta}_{*}=102/221$を代入し、計算を行う。
import numpy as np
t_statistic = (40./114.-62./107.)/np.sqrt((1./114.+1./107)*(102./221.)*(1.-102./221.))
print("t_statistic: {:.3f}".format(t_statistic))・実行結果
> print("t_statistic: {:.3f}".format(t_statistic))
t_statistic: -3.406参考
・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1