scipy.stats: 適合度の検定 chisquare

# scipy.stats: 適合度の検定 chisquare

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

[^1]: `scipy.stats.binom_test`, `scipy.stats.power_divergence`, `scipy.stats.chi2_test` との関連についても述べる。

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

特定の分布としては，

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


```python
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 の目は均等に出ているといえるか？


```python
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


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


```python
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]))




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




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




```python
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` が使われている。


```python
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 であるといえるか?


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




    Power_divergenceResult(statistic=4.266666666666667, pvalue=0.03886710381241731)




```python
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)` である。


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

    88.3125 29.4375 29.4375 9.8125



```python
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$ は，データの平均値に一致する。


```python
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!}$
で求める。


```python
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
```

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


```python
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 未満なので合併する。


```python
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]
```

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


```python
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` を指定する。 


```python
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. 無理やり応用して，クロス集計表の独立性の検定

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


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


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


```python
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` で求めて，


```python
from scipy.stats.contingency import expected_freq
```


```python
e = expected_freq(tbl)
e
```




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




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

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


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




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