適合度検定

適合度検定と一言に言っても多様な状況と多様な検定方法がある。
今回はその方法と実例と証明を順におっていこう。
複数のパターンについて、別の記事で紹介していく。

多項分布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$なら、正規分布と仮定して問題ない。

最尤法による検定

今後の課題とする