カイ二乗分布¶
import matplotlib.pyplot as plt
from scipy.special import comb
カイ二乗検定(Chi-squared test, $\chi^2$ 検定) とは、「帰無仮説が正しければ検定統計量が漸近的にカイ二乗分布に従うような統計的検定法の総称」であり、以下で表される $\chi^2$ 値を用いる。
※ なお、一般に事象の真の期待値はわからないので、「観測値から推定した値」や「理論値」を用いる。
ピンポン玉実験¶
それでは実際に、カイ二乗分布を求めるために、簡単な事象を例に考えてみる。
白とオレンジのピンポン玉をそれぞれ $50$ 個ずつ箱に入れておき、無作為に $10$ 個取り出す。調べたら、また取り出した分を箱の中に戻して、同じように無作為に $10$ 個取り出す。
この時、白とオレンジのピンポン玉数の期待値はそれぞれ $5$ 個であり、それを用いるとカイ二乗値は以下のように計算できる。
def Calc_Chi_squared(n, orange):
"""
関数の概要：カイ2乗値を求める関数
@param n ：取り出すピンポン玉の総数
@param orange：取り出したオレンジ色の球の数
"""
white = n-orange
mean = n/2
return ((mean - orange)**2 + (mean - white)**2) /mean
N = 100 # ピンポン玉の総数
n = 10 # 一度に取り出す玉の総数
Chi_squared = [Calc_Chi_squared(n, num_orange) for num_orange in range(n+1)]
Chi_squared
したがって、以下のような表が書ける。
| 白とオレンジの個数 | カイ二乗値 |
|---|---|
| 5個と5個 | 0 |
| 6個と4個 | 0.4 |
| 7個と3個 | 1.6 |
| 8個と2個 | 3.6 |
| 9個と1個 | 6.4 |
| 10個と0個 | 10.0 |
ここで、左の事象が起こる確率は実際に以下のように計算することができる。
def Calc_Prob(N, n, orange):
"""
関数の概要：事象が起こる確率を求める関数
@param N ：球の総数(オレンジと白は同じ数と仮定)
@param n ：取り出すピンポン球の総数
@param orange：取り出したオレンジ色の球の数
"""
white = n-orange
return comb(N/2, orange, exact=True) * comb(N/2, white, exact=True) / comb(N, n, exact=True)
Prob = [Calc_Prob(N, n, num_orange) for num_orange in range(n+1)]
Prob
よって、先ほどの表が以下のように更新できる。
| 白とオレンジの個数 | カイ二乗値 | 事象が起こる確率 |
|---|---|---|
| 5個と5個 | 0 | 0.259 |
| 6個と4個 | 0.4 | 0.211 |
| 7個と3個 | 1.6 | 0.113 |
| 8個と2個 | 3.6 | 0.0380 |
| 9個と1個 | 6.4 | 0.00724 |
| 10個と0個 | 10.0 | 0.000593 |
以上により、カイ二乗値と確率の関係性が得られた。
※ なお、気をつけるべきことは、「カイ二乗値=0.4となる確率は0.211ではない」ということである。白玉4個とオレンジ玉4個でもカイ二乗値は0.4となるので、2倍する必要がある！！
Chi_squared[int(n/2):]
# 単純に、2倍しているだけ。
Probability = [p*2 if i!=0 else p for i,p in enumerate(Prob[int(n/2):])]
Probability
また、確率密度関数(probability density function, PDF)を求めたいので、確率を足し合わせていく。
Probability_density = [sum(Probability[i:]) for i in range(len(Probability))]
Probability_density
したがって、これらの関係をプロットすると、以下のようになる。
plt.plot(Chi_squared[int(n/2):], Probability_density, label="degree of freedom = $1$")
plt.title("$\chi^2$ distribution")
plt.ylabel("Probability density")
plt.xlabel("$\chi^2$")
plt.legend()
plt.show()
かなりガタガタの分布となっている:(
そこで、同じ操作を $n=50$ で行うと、以下の分布が得られる。
n = 50
Chi_squared = [Calc_Chi_squared(n, num_orange) for num_orange in range(n+1)]
Prob = [Calc_Prob(N, n, num_orange) for num_orange in range(n+1)]
Probability = [p*2 if i!=0 else p for i,p in enumerate(Prob[int(n/2):])]
Probability_density = [sum(Probability[i:]) for i in range(len(Probability))]
plt.plot(Chi_squared[int(n/2):], Probability_density, label="degree of freedom = $1$")
plt.title("$\chi^2$ distribution")
plt.ylabel("Probability density")
plt.xlabel("$\chi^2$")
plt.xlim(0,8)
plt.legend()
plt.show()
そこそこ綺麗に求めることができた。
カイ二乗検定¶
上で求めたカイ二乗分布のカイ二乗値が $0.384$ のところでグラフを区切ると、それよりも左側の面積が $0.95$、右側の面積が $0.05$ となる。
したがって、カイ二乗値が $3.84$ 以下の値は、帰無仮説の下では $0.95$ の確率で起き、カイ二乗値が $0.38$ 以上の値は $0.05$ の確率で起きる、と言える。
一般に $\alpha=0.05$ を有意水準と用いる場合が多い（事前に決めておく）。すると、例えばカイ二乗値 $4.46$ が出てくる確率は帰無仮説が正しいとすると有意水準以下であり、「滅多に起こらない事象」である。
つまり、偶然の誤差ではないと認めることができ、帰無仮説を棄却することができる。
仮説検定（カイ二乗検定）において計算されるP値は、$p=\mathrm{P}(X_{\mathrm{obs}}\leq X)$ のことである。
観察された検定統計量 $X_{\mathrm{obs}}$（カイ二乗値）を用いて確率密度関数を $f_x$ とすると、以下の値がP値となる。
$$\mathrm{P}(X_{\mathrm{obs}}\leq X) = \int_{X_{\mathrm{obs}}}^{\infty}f_x(x)dx$$※ 難しく書いてありますが、確率密度関数を使っているので当たり前です。
自由度¶
ここまで考えたのは自由度1のカイ二乗分布であった。一般に自由度 $k$ のカイ二乗分布に拡張することが可能であり、以下の式で定義される。
確率変数 $Z_1,Z_2,\ldots,Z_n$ が互いに独立であり、それぞれが標準正規分布 $N(0,1)$ に従うとき、 $$\chi^2 = Z_1^2+Z_2^2+\cdots+Z_k^2$$ の $\chi^2$ に従う分布を、自由度 $k$ のカイ二乗分布と言う。
自由度 $\phi$ のカイ二乗分布は、以下で表される。
