Pythonプログラミング（ステップ７・G検定・適合度検定）

このページでは、ピアソンのカイ二乗検定に代わって使われるようになっているG検定について考える。

分布の違いを測る指標

$n$種類の事象がある場合に、モデルや観測から得られる確率分布が、一方は $$ P = \{p_1, p_2, \cdots, p_n\} $$ 他方が $$ Q = \{q_1, q_2, \cdots, q_n\} $$ であった（添字は事象を表すとする）とすると、$P$を基準として$Q$との分布の「差異」を表す指標としてKullback–Leibler divergenceが用いられる： $$ D_{\textrm{KL}}(P \Vert Q) = \sum_{i=1}^n p_i \log\left( \frac{p_i}{q_i} \right) $$

ふたつの分布が等しい場合は、対数の中が1になることから、$D_{\textrm{KL}}=0$となることは明らかである。 分布に違いがあると、$D_{\textrm{KL}}$は正値を取る。

適合度検定（G検定）

実験等によって得られる総数$N$サンプルが$n$個のクラスに分類できるような場合を考えよう。 例えば、「赤」か「白」の玉の入った箱の中から無作為に$N$個取り出す際は、クラスは「赤」と「白」の２種類ということになる。

実際には、「赤」が6割、「白」が4割だったとし、これを $$ \{ p_1 = 0.6, p_2 = 0.4 \} $$ と表すことにする。 一方で、赤と白に偏りが無いという帰無仮説を採用するとし、これを $$ \{ q_1 = 0.5, q_2 = 0.5 \} $$ と表現すれば、KL divergenceは $$ 0.6 \log_2\left( \frac{0.6}{0.5} \right) + 0.4 \log_2\left( \frac{0.4}{0.5} \right) \approx 0.029 \, \textrm{[bit]} $$ となる。

ここで、実際の計測を想定して、$i$番目のクラスの計測回数を$o_i$、 モデル（帰無仮説）から期待される回数を$e_i$とする。 $p_i = o_i/N, \, q_i = e_i / N$である。 そして、期待される回数と実測とのずれを $\delta_i = o_i - e_i$とし、相対的なずれ($\delta_i/e_i$)が十分に小さいと仮定して対数関数を展開すれば、 $$ \begin{eqnarray} D & = & \sum_i \left( \frac{o_i}{N} \right) \log\left( \frac{o_i/N}{e_i/N} \right) \\ & = & \sum_i \left( \frac{e_i + \delta_i}{N} \right) \log\left( \frac{e_i + \delta_i}{e_i} \right) \\ & = & \frac{1}{N} \sum_i (e_i + \delta_i) \log\left( 1 + \frac{\delta_i}{e_i} \right) \\ & \approx & \frac{1}{N} \sum_i (e_i + \delta_i) \left( \frac{\delta_i}{e_i} - \frac{1}{2} \left(\frac{\delta_i}{e_i}\right)^2 + \cdots \right) \\ & = & \frac{1}{N} \sum_i \left( \delta_i + \frac{1}{2} \frac{{\delta_i}^2}{e_i} + \cdots \right) \\ & = & \frac{1}{2N} \sum_i \frac{(o_i - e_i)^2}{e_i} + \cdots \end{eqnarray} $$ となる。 ここで、$\sum_i \delta_i = \sum_i ( o_i -e_i ) = N-N = 0$を用いた。

以上から、KL divergenceは、ピアソンのカイ二乗統計量と $$ D_{\textrm{KL}}( O \Vert E) \approx \frac{1}{2N} \sum_i \frac{(o_i - e_i)^2}{e_i} $$ のように関係づけられる。

そこで、改めて、 $$ G = 2 N D_{\textrm{KL}}( O \Vert E) = 2 \sum_i o_i \log\left( \frac{o_i}{e_i} \right) $$ を定義すると、この$G$は、ピアソンのカイ二乗統計量と同様な統計（分布と自由度）に従うと期待できるため、 この$G$を使って適合度検定を行うことができる。 （サンプル数が少ない場合でも、むしろ、カイ二乗統計量よりもカイ二乗分布に近い振る舞いをすると言われている）。

比較実験

最近の文献では、カイ二乗検定ではなくてG検定を適合度検定に用いることが推奨されているが、実際のところ、両者はどの程度異なるのか、 シミュレーションで様子を比較してみよう。

赤と白の玉がそれぞれ500個ずつ(合計1000個）入った袋から30個の玉を取り出す試行を20000回繰り返し、 カイ二乗値とG値の分布をプロットするコードの例を以下に示す （玉の取り出しになぞらえて、冗長ではあるが、bagというリストを敢えて用意している）。 このコードでは、自由度１のカイ二乗分布も併せてプロットしている（赤色の線）。 グラフを見る限り、このような状況では、両者に大きな差異はなさそうに思える。 

# coding: utf-8

import math
import random
from scipy.stats import chi2
import numpy as np
import matplotlib.pyplot as plt

reds = ['r']*500
whites = ['w']*500
bag = reds + whites

x2list=[]
glist=[]

N=30 # サンプル数

cnt=0 
while cnt<20000:
    # sample = random.choices(bag,k=N) # 重複を許す場合
    sample = random.sample(bag,k=N) # 重複を許さない場合
    n_white = sum(x=='w' for x in sample)
    n_red = sum(x=='r' for x in sample)
    x2 = (n_red-N*0.5)**2/(N*0.5) + (n_white-N*0.5)**2/(N*0.5)
    x2list.append(x2)
    g = 2 * (n_red * math.log(n_red/(N*0.5)) + n_white * math.log(n_white/(N*0.5)))
    glist.append(g)
    cnt += 1

df=1
rv = chi2(df)
x = np.linspace(0, 15, 100)

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.hist(glist, bins=30, range=(0,15), density=True, color='green')
ax1.plot(x, rv.pdf(x), 'k-', lw=2, color='red')
ax1.set_xlabel('G')
ax1.set_ylabel('Freq')
ax1.grid(True)

ax2.hist(x2list, bins=30, range=(0,15), density=True, color='blue')
ax2.plot(x, rv.pdf(x), 'k-', lw=2, color='red')
ax2.set_xlabel('X^2')
ax2.grid(True)

plt.show()

icon-pc　練習： 自由度による違いの確認

上のコードを参考に、玉の種類がより多い場合について、比較実験を行ってみなさい。

また、サンプリングの際の重複の効果についても確認しなさい。