適合度検定
適合度検定と一言に言っても多様な状況と多様な検定方法がある。
今回はその方法と実例と証明を順におっていこう。
複数のパターンについて、別の記事で紹介していく。
- 多項分布1元データ
- 多項分布2元データ(確率検定、独立性検定)
- 群データ
- 質的データ
まずは、多項分布1元データパターンについて追っていく。
多項分布1元データパターン
各種類の出現確率がわかっている状況、ガチャと同じ。回した回数と確率を掛けることで理論度数がでる。
検定方法
\sum_{i=1}^{5} \frac{(x_i -n_i p_i)^2}{n_i p_i} \sim \chi(4)
ということがわかっていることから検定が可能。（中心極限定理で正規分布に近似しているので、目安として$n_i p_i > 5$は必要）
検定に必要な条件から、出現確率$1$%なら$500$回は回す必要がある。
証明
に書いてある。要点をまとめる。
ある種類aの度数$x_a$は、出現確率$p$と非出現確率確率の2項分布であることから
x_a \sim B(n,p)
と言える。すべての種類に同じことを言うと、どの種類でも出現しない、同時に$a,b$どちらにも該当する可能性がある。
つまり、すべてをまとめて考える場合、二項分布であると考えるのは間違いであるということである。
今回のモデルは多項分布と呼ばれる。
今回のデータ
についていれば
{}_{25}C_{3}\;
{}_{22}C_{4} \;
{}_{18}C_{5} \;
{}_{13}C_{6} \;
{}_{7}C_{7}\;
p_a^3\;
p_b^4\;
p_c^5\;
p_d^6\;
p_e^7\;
が今回の出現確率である。(25回中3つ出現する場所を選ぶ)(aが3回出現する確率)（残りの22回から4回出現する場所を選ぶ）・・・・
という考え方です。
出現回数を$q_i$として、一般的に書く場合、$q_{1}\cdots q_n$が出現する確率は
f(q_{1}\cdots q_n) = \frac{n!}{q_1! \cdots q_n!} p_1^{q_1} \cdots p_n^{q_n}
になります。
同時に$a,b$どちらにも該当する可能性はないこれは、$q_i,q_j$２つには相関があることを意味する。
証明
x_i = \frac{(q_i-np_i)}{\sqrt{np_i}}
を新たに定義する。期待値は0になる。なぜならば、各種類一つのみに注目する場合は2項分布であるからである。分散は$(1-p_i)$となる。
束縛条件がある
\sum_{i=1}^{n} \sqrt{p_i} x_i = \sum \frac{(q_i-np_i)}{\sqrt{n}} = \frac{(\sum q_i-n \sum p_i)}{\sqrt{n}} =0
$x_i$を独立な変数$z_i$に変換する必要がある。
正規直交行列$A= (a_{ij})$について
a_{nj} = \sqrt{p_j}\\
とする。
つまり、
こんな感じで想定する。
\begin{eqnarray}
z_i &=& \sum a_{ij} x_j\\
\therefore
\begin{pmatrix}
z_1\\
\vdots\\
z_n
\end{pmatrix}
&=&A
\begin{pmatrix}
x_1\\
\vdots\\
x_n
\end{pmatrix}\\
\therefore
z_n &=& \sqrt{p_1} x_1 + \cdots +\sqrt{p_n} x_n = 0 \\
z_i(i=1\cdots n-1 )&=& \sum _{j = 1} ^{n-1} a_{ij} x_j + a_{in} x_n\\
&=& \sum _{j = 1} ^{n-1} a_{ij} x_j + a_{in} (0 -\sqrt{p_1} x_1 - \cdots -\sqrt{p_{n-1}} x_{n-1} ) / \sqrt{p_{n}} \\
&=& \sum _{j = 1} ^{n-1} (a_{ij} - a_{in} \frac{\sqrt{p_j}}{\sqrt{q_n}}) x_j
\end{eqnarray}
前述の通り$x_i$の期待値は0,分散は$(1-p_i)$となる。よって、$z_i$の期待値も0になる。
$z_i$がそれぞれ、独立で分散が1であることを証明できれば、中心極限定理で$z_i$は標準正規分布であると言える。
最初の式
\begin{eqnarray}
\sum_{i=1}^{5} \frac{(q_i -n_i p_i)^2}{n_i p_i}
&=& \sum_{i=1}^{5} x_i^2 =
\begin{pmatrix}
x_1\\
\vdots\\
x_n
\end{pmatrix}^T
\begin{pmatrix}
x_1\\
\vdots\\
x_n
\end{pmatrix}\\
&=&
\begin{pmatrix}
z_1\\
\vdots\\
z_n
\end{pmatrix}^T
(A^T)^{-1}
A^T
\begin{pmatrix}
z_1\\
\vdots\\
z_n
\end{pmatrix}\\
&=&
\begin{pmatrix}
z_1\\
\vdots\\
z_n
\end{pmatrix}^T
\begin{pmatrix}
z_1\\
\vdots\\
z_n
\end{pmatrix}\\
&=& \sum_{i=1}^{5} z_i^2(z_5 = 0)
\sim \chi^2(4)
\end{eqnarray}
となるため、カイ二乗分布に従うことが示せる。
Aの消去はAの正規直交性を使った。
ziの独立性
\begin{eqnarray}
E[x_i x_j]&=& E[\frac{(q_i-np_i)}{\sqrt{np_i}} \frac{(q_j-np_j)}{\sqrt{np_j}} ]\\
&=& \frac{1}{n^2 p_i p_j } (E[q_i q_j]) + n^2 p_i p_j -2n^2p_i p_j)
\end{eqnarray}
である。
\begin{eqnarray}
E[q_i q_j] &=& n(n-1)p_i p_j
\end{eqnarray}
らしい。後で証明する。
\begin{eqnarray}
E[x_i x_j]&=& \frac{1}{n \sqrt{p_i p_j} } (n(n-1)p_i p_j + n^2 p_i p_j - 2n^2p_i p_j) \\
&=&\frac{1}{n \sqrt{p_i p_j} } (-np_i p_j)\\
&=& -\sqrt{p_ip_j}
\end{eqnarray}
となる。
目標である。$z_i$については
\begin{eqnarray}
E[z_i z_k] &=& E[\sum_{i=1}^n a_{ij} x_j \sum_{m=1}^n a_{km} x_m]\\
&=& E[\sum_{j=1}^n \sum_{m=1}^n a_{ij} a_{km} x_j x_m]\\
&=& -\sum_{j=1}^n \sum_{m=1}^n a_{ij} a_{km} \sqrt{p_j p_m}\\
&=& -\sum_{j=1}^n \sum_{m=1}^n a_{ij} a_{km} a_{nj} a_{nm}(\because Aの定義)\\
&=& -\sum_{j=1}^n a_{ij} a_{nj} \sum_{m=1}^n a_{km} a_{nm}\\
&=& 0(i \neq n , k \neq n)
\end{eqnarray}
最後の式展開はAの直交性を使った。細かく書くと
(AA^*)_{ij} = \sum_{k = 1}^n a_{ik} a_{jk} =
\left\{
\begin{array}{ll}
0&(i \neq j)\\
1&(i = j)\\
\end{array}
\right.
である。
$i=n,k=n$の場合は$z_n =0 $であることから$E[z_n z_k] = 0$になる。
よって、$z_i$は独立である。
ziの分散
分散も求めてみよう。
\begin{eqnarray}
V[z_i] &=& V[\sum_{j=1}^n a_{ij} x_j]\\
&=&\sum_{j=1}^n a_{ij}^2 V[x_j] + \sum_{k \neq m} a_{ik} a_{im} E[x_k x_m]\\
&=&\sum_{j=1}^n a_{ij}^2 (1- p_j) - \sum_{k \neq m} a_{ik} a_{im}\sqrt{p_k p_m}\\
&=&\sum_{j=1}^n a_{ij}^2 - \sum a_{ik} a_{im}\sqrt{p_k p_m}(\because シグマ内の除外部分がちょうど移動できた)\\
&=&\sum_{j=1}^n a_{ij}^2 - (\sum_k a_{ik} \sqrt{p_k})^2\\
&=&\sum_{j=1}^n a_{ij}^2 - (\sum_k a_{ik} a_{nk})^2\\
&=&
\left\{
\begin{array}{ll}
0&(i \neq j)\\
1&(i = j)\\
\end{array}
\right.
\end{eqnarray}
となる。
度数の共分散
\begin{eqnarray}
E[q_i q_j] &=& n(n-1)p_i p_j
\end{eqnarray}
を計算の途中で利用した。この証明をすれば証明完了である。
多項分布を考えるときの分散・共分散の計算は二項分布の場合と同じように$(x+y+z)^n$のような式の展開を考えると良い
(x+y+z)^n = \sum_{a,b} {}_nC_a \; {}_{n-a}C_b x^a y^b z^{n-a-b}
となる。
$q_i,q_j$をとる確率は
f(q_{i},q_{j}) = {}_n C_{q_{i}} \; {}_{n-q_i} C_{q_j} p_i^{q_i} p_j^{q_j} {(1-p_i-p_j)}^{n-q_{i}-q_{j}}
である。
ｘ，ｙ，ｚの式に$x = p_i,y=p_j,z=1-p_i-p_j,a=q_i,b=q_j$を代入するとシグマの中身が$f(q_{i},q_{j})$になる。よって
\sum_{i,j} f(q_{i},q_{j}) = (x+y+z)^n = 1
が証明できる。
xy\frac{\partial^2}{\partial x \partial y}(x+y+z)^n = xy n(n-1)(x+y+z)^{n-2}\\
となることを利用して、
\begin{eqnarray}
E[q_i q_j] &=& \sum_{q_i,q_j} q_i q_j f(q_{i},q_{j}) \\
&=&\sum_{q_i,q_j} q_i q_j \; {}_n C_{q_{i}} \; {}_{n-q_i} C_{q_j} p_i^{q_i} p_j^{q_j} {(1-p_i-p_j)}^{n-q_{i}-q_{j}} \\
&=&\sum_{a,b} ab \; {}_n C_{a} \; {}_{n-a} C_{b} x^{a} y^{b} z^{n-a-b} \\
&=&\sum_{a,b} \; {}_n C_{a} \; {}_{n-a} C_{b} \frac{\partial^2}{\partial x \partial y} \{x^{a} y^{b} z^{n-a-b} \}xy (\because 微分で減った分のx,yの次数を別に掛けることで補填する)\\
&=&xy \frac{\partial^2}{\partial x \partial y} \{ \sum_{a,b} \; {}_n C_{a} \; {}_{n-a} C_{b} x^{a} y^{b} z^{n-a-b} \}\\
&=&xy \frac{\partial^2}{\partial x \partial y} \{ (x+y+z)^n \}\\
&=&xy n(n-1) (x+y+z)^{n-2}\\
&=&p_i p_j n(n-1)
\end{eqnarray}
証明終了。
数値実験
次に、Pythonで多項分布を模擬して、カイ二乗検定量を計算して、分布を求めてみよう。
総度数$100$
試行回数$1000$回
実験結果
わかること
種類aの度数分布の平均はちゃんと確率×総度数=50に来ている。
各種類の度数分布の分散は異なる値になっている。
数式的にいれば、二項分布の分散$np(1-p)$に一致するので
標準偏差の値とグラフの６８％点は近い気がする
試行回数のみを1000回から10000回に十倍すると
ほとんど変わっていないことがわかる。試行回数を増やしても、出現度数$q_i$の分布は変わらないので当たり前ではある。
作成したプログラム
import random
import numpy as np
import matplotlib.pyplot as plt
p = np.array([1/2,1/4,1/5])
sum = 0
for i in p:
sum = sum + i
if sum > 1:
print("error")
else:
p = np.append(p,1-sum)
n = 100
sample = 10000
result = np.zeros((p.size,sample))
for sample_num in range(sample):
for j in range(n):
d = random.uniform(0, 100)
sum = 0
rank = 0
for i in p:
sum = sum + i
if d/100< sum:
break
rank = rank + 1
result[rank,sample_num] = result[rank,sample_num] + 1
#各列で平均と分散を計算
# 箱ひげ図
fig, ax = plt.subplots()
points=[]
for i in range(result.shape[0]):
points.append(result[i])
bp = ax.violinplot(points)
plt.grid() # 横線ラインを入れることができます。
# 描画
plt.show()
print("p----------")
print(p*n)
print("result-----")
print(result)
カイ二乗検定量を計算
カイ二乗検定量を計算する。
ちゃんとカイ二乗分布してるように見える。
確率を他のパターンにして見た結果
カイ二乗分布してる。
作成したプログラム
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
p = np.array([1/2,1/4,1/5])
sum = 0
for i in p:
sum = sum + i
if sum > 1:
print("error")
else:
p = np.append(p,1-sum)
n = 100
sample = 100000
result = np.zeros((p.size,sample))
kai = np.zeros(sample)
for sample_num in range(sample):
for j in range(n):
d = random.uniform(0, 100)
sum = 0
rank = 0
for i in p:
sum = sum + i
if d/100< sum:
break
rank = rank + 1
result[rank,sample_num] = result[rank,sample_num] + 1
for sample_num in range(sample):
kai_num = 0
for i in range(p.size):
x = result[i,sample_num]
riron = n * p[i]
kai_num = kai_num + (x - riron)**2 / riron
kai[sample_num] = kai_num
#各列で平均と分散を計算
# 箱ひげ図
fig, ax = plt.subplots()
points=[]
for i in range(result.shape[0]):
points.append(result[i])
#ax.violinplot(points)
#ax.boxplot(points,labels = ['a','b','c','d'])
# 描画
bin_num = 100
ax.hist(kai,range = (0,30),bins = 50,color = "skyblue",density = True)
#カイ二乗グラフ
x = np.linspace(0,30,300)
y = chi2.pdf(x,p.size-1)
ax.plot(x,y,color = "pink")
#ax.grid() # 横線ラインを入れることができます。
plt.show()
print("p----------")
print(p*n)
print("result-----")
print(result)
中心極限定理による近似の限界について
x_i = \frac{(q_i-np_i)}{\sqrt{np_i}} \\
z_i = \sum a_{ij} x_j
おいて、$q_i$は二項分布であったが、$z_i$を正規分布として近似して以降を議論した。
中心極限定理は母集団がどんな分布でも（平均$\mu$、分散$\sigma^2$）、$n$回の平均値の分布は正規分布$N(\mu,\sigma^2 /n)$に収束する定理である。
二項分布はベルヌーイ分布＊試行回数の分布である。つまり、平均値であるので回数が増えるほど正規分布に収束していく。
の証明でも問題ない。
nがどれくらいまでなら近いと言えるか
ピンクの線が近似した正規分布$N(n p_i ,np_i (1-p_i))$、青線が二項分布（$p = 2/100,n=100$）に従う確率変数$x$を試行回数$10000$回で計測したときの密度関数
二項分布（$p = 10/100,n=100$）に変更して実験した結果
二項分布（$p = 50/100,n=100$）に変更して実験した結果
近いかどうかというより、正規分布は負の値も取るが二項分布は正の値しか取らない。
結果、0以下の部分が見切れることになる。
見きれないようにするには、$n p_i$の値がある程度以上大きくなる必要がある。$3\sigma$の値は
n p_i -3 \sqrt{n p_i (1-p_i)}
である。
によれば、３シグマで99.7％はいる、2シグマで95.4％入る、1シグマで68.3%はいる。
$p_i$が小さければ$1-p_i \approx 1$といえることを利用して、各シグマ範囲を計算すると
つまり、$np_i \geq 4$であれば、$95.4$%のデータは収まっていることになる。
若干、ヒストグラムが右側に振れているように見えるが、これは横軸$5$に対する縦軸の値は$5 \leq x < 6$に該当するデータの数を表してるためである。
作成したプログラム
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
p = np.array([50/100])
sum = 0
for i in p:
sum = sum + i
if sum > 1:
print("error")
else:
p = np.append(p,1-sum)
n = 100
sample = 10000
result = np.zeros((p.size,sample))
kai = np.zeros(sample)
for sample_num in range(sample):
for j in range(n):
d = random.uniform(0, 100)
sum = 0
rank = 0
for i in p:
sum = sum + i
if d/100< sum:
break
rank = rank + 1
result[rank,sample_num] = result[rank,sample_num] + 1
#各列で平均と分散を計算
# 箱ひげ図
fig, ax = plt.subplots()
points=[]
range1 = (0,100)
ax.hist(result[0],range = range1,bins = range1[1],color = "skyblue",density = True)
#正規分布グラフ
mu = n * p[0]
sigma = n*p[0]*(1-p[0])
x = np.linspace(0,100,300)
y = np.exp(-(x-mu)**2 / 2 / sigma) / np.sqrt(2 * np.pi*sigma)
ax.plot(x,y,color = "pink")
ax.grid() # 横線ラインを入れることができます。
plt.show()
print("p----------")
print(p*n)
print("result-----")
print(result)
まとめ
二項分布を正規分布として仮定する場合、データが見切れる場合がある。
$np_i \geq 4$なら、正規分布と仮定して問題ない。
最尤法による検定
今後の課題とする