カイ二乗値を計算する

　カイ二乗値を計算して、カイ二乗検定とはそういうことだったのかと今さら納得した。

まず、偏りのないサイコロを振ってカイ二乗値を計算してみる。

import bisect

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

np.random.seed(1226)

# create a dice with given frequency
def create_dice(freq):
    cdf = np.cumsum(freq)
    def dice():
        p = np.random.uniform()
        return bisect.bisect_left(cdf, p)
    return dice

# calculate chi-squared value
def calc_chi_sq(freq, sample_size):
    dice = create_dice(freq)
    sample = np.bincount([dice() for _ in xrange(sample_size)])
    populatoin = [sample_size * f for f in freq]
    return np.sum([(x-y)**2/y for (x, y) in zip(sample, populatoin)])

freq = [1./6, 1./6, 1./6, 1./6, 1./6, 1./6]

# plot dice histogram
dice = create_dice(freq)
plt.hist([dice() for _ in range(3000)], bins=range(7), normed=True)
plt.show()

# plot chi-squared
chi_sqs = [calc_chi_sq(freq, 3000) for _ in xrange(3000)]
plt.hist(chi_sqs, bins=50, normed=True)
x = np.linspace(0,35,50)
y = chi2.pdf(x, df=5)
plt.plot(x, y, color='r')
plt.show()

実行したら以下の２つのグラフが出力される。
上がサイコロの目の分布。下がカイ二乗値の分布。

カイ二乗値はカイ二乗分布に従っていることが分かる。

続いて、激しく偏った分布を持つサイコロを振って、同じことをしてみる。
上のコードの

freq = [1./6, 1./6, 1./6, 1./6, 1./6, 1./6]

の部分を

freq = [3./12, 1./12, 4./12, 2./12, 1./12, 1./12]

に変えて実行してみると、

なるほど、母集団の分布に関係なくカイ二乗値はカイ二乗分布に従うと。

ということで、カイ二乗検定が何をやってるのか今さらながら納得できた。。。