scipy.stats: 適合度の検定 chisquare

Help では "Calculates a one-way chi square test" と書いているが，一般的には適合度の検定である。¹

カテゴリーデータの集計の結果，度数分布が特定の分布にしたがっているかどうかの検定である。

特定の分布としては，

一様分布
帰無仮説：各カテゴリーの頻度は同じである $H_0: p_i = 1/k,\ \ i=1,2,\dots,k$
その他の分布
帰無仮説：各カテゴリーの頻度は特定の割合（母比率 $\pi_i$）である $H_0: p_i = \pi_i,\ \ \sum \pi_i = 1,\ \ i=1,2,\dots,k$

from scipy.stats import chisquare

chisquare(f_obs, f_exp=None, ddof=0, axis=0)

f_obs は観察度数
f_exp は期待度数。サンプルサイズ=$n$，カテゴリー数=$k$ のとき，デフォルトで，すべてのカテゴリーの期待度数は $n/k$ に等しい。
ddof 自由度は $k-1$ から更に余分に減少すべき自由度。デフォルトは 0。

1.　一様分布に従うか？

1.1.　カテゴリー数が 3 以上の場合

よくある例で，サイコロを振ったときに 1〜6 の目が同じ頻度で出るかの検定である。

普通のサイコロのデータではありふれているので，ここでは現実世界には存在しない正 10 面体のサイコロを使ってみる。

例題：　正 10 面体のサイコロを 500 回振って，出た目を記録した。1 〜 10 の目は均等に出ているといえるか？

import numpy as np
np.random.seed(123)
x = np.random.randint(1, 11, 500) # 1 ~ 10 の整数乱数を 500 個生成
print(*x)

3 3 7 2 4 10 7 2 1 2 10 1 1 10 4 5 1 1 5 2 8 4 3 5 8 3 5 9 1 8 10 4 5 7 2 6 7 3 2 9 4 6 1 3 7 3 5 5 7 4 1 7 5 8 7 8 2 6 8 10 3 5 9 2 3 2 2 4 6 10 1 9 2 7 4 4 6 10 8 10 3 4 4 4 9 7 10 8 7 4 10 7 7 7 2 4 5 4 2 1 6 9 7 9 10 2 1 4 2 4 5 8 7 2 5 4 4 8 7 9 7 5 5 8 1 1 10 9 9 5 9 7 2 7 9 8 10 2 8 2 8 10 9 8 2 4 2 9 8 6 2 3 6 3 3 10 4 3 7 8 10 2 4 9 4 8 10 10 4 4 6 7 1 9 8 8 5 5 6 1 9 10 3 6 2 6 10 3 5 4 1 4 8 8 3 6 2 8 6 10 2 3 9 6 1 10 4 4 4 2 8 7 4 2 8 3 7 6 10 9 10 5 8 9 6 10 3 5 4 4 2 8 9 6 8 8 8 6 1 1 9 2 1 5 6 5 6 1 8 2 3 3 8 8 4 7 2 9 5 2 7 6 10 1 4 2 4 8 8 8 6 4 7 4 6 2 6 10 2 4 10 2 6 7 2 6 1 8 1 4 1 7 3 8 6 2 1 1 3 7 6 10 2 6 4 7 5 8 2 7 1 4 2 4 1 9 6 2 5 4 3 3 1 3 1 3 8 10 1 4 10 8 1 6 5 7 1 8 9 9 2 1 2 8 6 4 1 5 9 5 10 4 6 2 9 1 10 2 6 9 9 5 7 9 8 10 3 1 1 10 10 6 4 6 6 3 3 8 4 2 4 10 4 7 3 4 2 10 9 1 3 4 8 10 3 8 10 8 2 5 8 5 8 9 1 3 8 5 4 8 4 8 1 8 7 9 7 9 1 6 4 3 5 4 4 10 8 2 9 1 2 5 8 3 5 9 4 5 10 7 3 8 8 8 5 8 9 5 1 8 7 4 5 2 2 8 3 8 9 10 5 5 5 8 8 10 10 5 2 6 6 9 10 6 2 1 7 1 10 5 7 5 6 4 5 7 3 6 3 5 4 3 5 9 1 7 9 9 6 7

集計結果は以下のようになる。

counts = np.unique(x, return_counts=True)
counts

(array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10]),
 array([48, 56, 42, 62, 48, 45, 44, 65, 43, 47]))

f = counts[1] # 観察度数
f

array([48, 56, 42, 62, 48, 45, 44, 65, 43, 47])

from scipy.stats import chisquare
chisquare(f)

Power_divergenceResult(statistic=11.92, pvalue=0.21785699690374188)

帰無仮説は $H_0: p_i = 1/k,\ \ i=1,2,\dots,k$ である。
すなわち「1〜10 の出方は同じである」ということ。

有意水準 5% で検定するとき，$p$ 値が 0.21786 なので，帰無仮説を棄却できない。すなわち，「出た目の頻度は一様でないとはいえない」ということである。

間違っても「出た目の頻度は等しい」といってはいけない。

なお，分析結果の名前が Power_divergenceResult になっていることからも判るが，実際は　　power_divergence が使われている。

from scipy.stats import power_divergence
power_divergence(f)

Power_divergenceResult(statistic=11.92, pvalue=0.21785699690374188)

1.2.　カテゴリー数が 2 の場合

例えば，コイントスで表と裏が等確率（すなわち 0.5 ずつ）であるかのような場合にも，chisquare を使うことはできるが，二項検定 binom_test によれば，正確な検定（exact test）を行うことができる。

例題：　コインを 60 回投げて，表が 38 回，裏が 22 回出た。このコインの表・裏が出る確率はそれぞれ 0.5 であるといえるか?

from scipy.stats import chisquare
chisquare([38, 22], [30, 30])

Power_divergenceResult(statistic=4.266666666666667, pvalue=0.03886710381241731)

from scipy.stats import binom_test
binom_test(38, 60)

0.051893895928921435

漸近検定である chisquare では，有意水準 5% で帰無仮説は棄却された（$p = 0.0389$）。しかし，正確な検定である binom_test では帰無仮説は棄却されなかった（$p = 0.519$）。

なお，カテゴリー数が 2 の場合には二項分布に基づく正確な二項検定があるが，カテゴリー数が 3 以上の場合にも多項分布に基づく正確な検定がある(scipy.stats には用意されていない）。

2.　特定の理論分布に従うか。

これもよくある例で，交配実験結果がメンデルの法則にしたがうかというのがある。メンデルの法則を知っている必要はない。

例題：　4 つのカテゴリーの頻度が，9:3:3:1 になるはずのところ，観察度数が 91, 25, 34, 7 であった。これは理論分布にしたがっているといえるだろうか？

サンプルサイズは　sum([91, 25, 34, 7])=157，期待度数は 157 * np.array([9, 3, 3, 1])/16) である。

f = [91, 25, 34, 7]
e = 157 * np.array([9, 3, 3, 1])/16
print(*e)

88.3125 29.4375 29.4375 9.8125

chisquare(f, e)

Power_divergenceResult(statistic=2.2639773531493277, pvalue=0.5194569931782831)

帰無仮説は確率でいえば $H_0: p_1 = 9/16,p_2=3/16, p_3=3/16, p_4=1/16$ である。
あるいは頻度でいえば $H_0: f_1 = 157\cdot9/16,f_2=157\cdot3/16, f_3=157\cdot3/16, f_4=157\cdot1/16$ である。
有意水準 5% で検定するとき，$p$ 値が 0.51946 なので，帰無仮説を棄却できない。すなわち，「頻度は理論分布にしたがっていないとはいえない」ということである。

間違っても「理論分布にしたがっている」といってはいけない。

3.　データから理論分布のパラメータを推定して適合度の検定を行う

たとえばポアソン分布にしたがうと推定されるが，ポアソン定数はデータから推定する必要がある場合である。

例題：　観察度数は，階級 0 〜 10 に対して， [27, 61, 77,71,54, 35, 20, 11, 6, 2, 1] であった。ポアソン分布にしたがっているといえるだろうか?

ポアソン定数 $\lambda$ は，データの平均値に一致する。

import numpy as np
i = np.arange(0, 11)
frequency = np.array([27, 61, 77,71,54, 35, 20, 11, 6, 2, 1])
λ = sum(i * frequency) / sum(frequency)
print(λ)

2.9917808219178084

期待値 expected はポアソン分布の確率関数
$\displaystyle p(x) = \frac{\lambda^x\ e^{-\lambda}}{x!}$
で求める。

from math import exp, factorial
k = len(frequency)
probability = np.array([λ**z * exp(-λ) / factorial(z) for z in range(k)])
probability[k-1] = 1 - sum(probability[0:k-1])
n = sum(frequency)
expected = n * probability

　ここまでを表にしてみよう。

print("number of class =", k, ",  total sample size =", n)
print('\n i  frequency  probability    expected')
for i in range(k):
    print(f"{i:2d}  {frequency[i]:9d}  {probability[i]:11.7f}  {expected[i]:10.7f}")

number of class = 11 ,  total sample size = 365

 i  frequency  probability    expected
 0         27    0.0501980  18.3222567
 1         61    0.1501813  54.8161761
 2         77    0.2246548  81.9989922
 3         71    0.2240393  81.7743374
 4         54    0.1675691  61.1627236
 5         35    0.1002660  36.5970927
 6         20    0.0499957  18.2484133
 7         11    0.0213680   7.7993219
 8          6    0.0079910   2.9167327
 9          2    0.0026564   0.9695805
10          1    0.0010805   0.3943730

適合度検定の際に，下側，上側のいくつかの階級の期待値が 1 未満の場合にはそれぞれ下側，上側から期待値が1以上になるように階級を合併する。

i = 9, 10 の階級の期待度数が 1 未満なので合併する。

frequency[k-2] += frequency[k-1]
probability[k-2] += probability[k-1]
expected[k-2] += expected[k-1]
k = k-1
frequency = frequency[:k]
probability = probability[:k]
expected = expected[:k]

ここまでを表にしてみよう。

print(' i  frequency  probability    expected')
for i in range(k):
    print(f"{i:2d}  {frequency[i]:9d}  {probability[i]:11.7f}  {expected[i]:10.7f}")

 i  frequency  probability    expected
 0         27    0.0501980  18.3222567
 1         61    0.1501813  54.8161761
 2         77    0.2246548  81.9989922
 3         71    0.2240393  81.7743374
 4         54    0.1675691  61.1627236
 5         35    0.1002660  36.5970927
 6         20    0.0499957  18.2484133
 7         11    0.0213680   7.7993219
 8          6    0.0079910   2.9167327
 9          3    0.0037369   1.3639535

どの階級も，期待値が 1 以上になったので，検定を実行できるようになった。

自由度の設定において，データから $\lambda$ を推定したので，自由度は余分に 1 少なくなるので ddof=1 を指定する。

from scipy.stats import chisquare
chisquare(frequency, expected, ddof=1)

Power_divergenceResult(statistic=14.143749368122329, pvalue=0.0780947329266041)

$p > 0.05$ なので，帰無仮説は棄却できない（$p = 0.078$）

4. 無理やり応用して，クロス集計表の独立性の検定

例題：　以下のようなクロス集計表で，独立性の検定を行う。

tbl = np.array([[24, 13, 7], [5, 30, 19]])

普通は， `scipy.stats.chi2_contingency` を使う。

from scipy.stats import chi2_contingency
chi2_contingency(tbl)

(23.93649410229202,
 6.342439710555313e-06,
 2,
 array([[13.02040816, 19.30612245, 11.67346939],
        [15.97959184, 23.69387755, 14.32653061]]))

chisquare で行う場合は，tbl の期待値を scipy.stats.contingency.expected_freq で求めて，

from scipy.stats.contingency import expected_freq

e = expected_freq(tbl)
e

array([13.02040816, 19.30612245, 11.67346939, 15.97959184, 23.69387755,
       14.32653061])

`ddof` を適切に指定すれば，`chi2_contingency` と同じ結果を得ることができる。

`tbl` と `e` を `np.ravel()` で一次元配列（ベクトル）にして使っているところに注目。

k, l = tbl.shape
chisquare(np.ravel(tbl), np.ravel(e), ddof=k+l-2)

Power_divergenceResult(statistic=23.93649410229202, pvalue=6.342439710555313e-06)

scipy.stats.binom_test, scipy.stats.power_divergence, scipy.stats.chi2_test との関連についても述べる。 ↩