観測データがポアソン分布に従うかを検定（Pythonによるポアソン分布の適合度の検定）


# 問題設定

あるスーパーマーケットにおいて、**（1日あたりの）顧客からのクレーム数**を、$365$ 日分集計したところ、次のような結果となった。

|クレーム数|0件|1件|2件|3件|4件|5件|6件|7件|
|:-------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:----:|
|観測度数<br>（日数）|22|49|82|86|63|39|16|8|

このデータは、**ポアソン分布に従っている**と考えてよいか？有意水準 $5\%$ で検定せよ。

# 可視化

ポアソン分布に従うと仮説を立てて計算した「**期待度数**」と「**観測度数**」をグラフにしてみます。

「**期待度数**」は次のように求められます。まず「クレーム数」と「**観測度数**」を掛けて、観測日数 $365$ で割ると、**1日あたりの平均クレーム数** $\mu=2.931$ が得られます。そして、平均 $2.931$ の**ポアソン分布**についての確率質量（$k=0..7$）を求め、観測日数 $365$ を掛けると「**期待度数**」が得られます。


```py:観測度数と期待度数の可視化
import numpy as np
import scipy.stats as st
import japanize_matplotlib
import matplotlib.pyplot as plt

n = np.arange(0,7+1)
f = np.array((22,49,82,86,63,39,16,8)) # 度数

mu = (n*f).sum()/f.sum() # 平均クレーム数を計算 2.931

Po = st.poisson(mu=mu) # μ=2.931 のポアソン分布
xp = Po.pmf(k=n)       # k=1...7 の確率質量  
xp[-1] = 1 - st.poisson.cdf(n[-2],mu=mu) # k>=7の確率質量をk=7に統合
fp = xp * f.sum()      # 期待度数

plt.figure(figsize=(6,3),dpi=120,facecolor='white')
plt.bar(n,f,width=0.9, label='観測度数')
plt.plot(n,fp,'o--', c='tab:orange', label=f'期待度数 $Po\,(\\lambda={mu:.2f})$')
plt.tick_params(axis='x',length=0)
plt.legend(framealpha=1, fancybox=False)
plt.gca().set_axisbelow(True)
plt.grid(axis='y')
plt.show()
```

実行結果は次のようになります。<br>![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/562e96a7-464e-3264-20a7-19f6107247f2.png)


|クレーム数|0件|1件|2件|3件|4件|5件|6件|7件|
|:-------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:----:|
|**観測度数**<br>$f$|22|49|82|86|63|39|16|8|
|**期待度数**<br>$f_{p}$|19.5|57.0|83.6|81.7|59.9|35.1|17.2|11.0|

# 適合度の検定

「対象がポアソン分布に従う」という仮説（帰無仮説）は正しいといえるのか？を**統計的検定**により評価します。

この仮説が正しい場合、**観測度数** $\boldsymbol{x}$ と **仮説に基づく期待度数** $\boldsymbol{x}_p$ より計算する次の **統計量** $t$ は（近似的に）**カイ2乗分布**（自由度 $n-2$）に従う、という性質を利用して検定を行ないます。

$$ t = \sum \frac{(x-x_p)^2}{x_p} $$

現在の問題では、クレーム数が $0$ 件から $7$ 件までの $8$ 区分あるので、**自由度**は $n-2=8-2=6$ になります。


自由度 $6$ のカイ2乗分布の確率密度関数を描いてみます。

```py:自由度6のカイ2乗分布
import numpy as np
import scipy.stats as st
import japanize_matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0,25,200)
p = st.chi2.pdf(x, df=6)

x_p95 = st.chi2.ppf(0.95, df=6)

y_range = (0.00,0.16)
plt.figure(figsize=(6,3),dpi=120,facecolor='white')
plt.plot(x,p)
plt.vlines(x_p95,*y_range,lw=1,ls=':',color='tab:red')

x2 = np.linspace(0,x_p95,200)
plt.fill_between(x2,st.chi2.pdf(x2,df=6),np.zeros(len(x2)),facecolor='tab:blue',alpha=0.3)

plt.text(5,0.05,'0.95',ha='center',va='center')
plt.text(x_p95,-0.004,f'{x_p95:.2f}',c='tab:red', ha='center',va='top')

plt.text(x_p95,0.08,f' → 棄却域',c='tab:red')

arrowprops=dict(arrowstyle='->',connectionstyle='angle,angleA=0,angleB=60, rad=1')
kw = dict(xycoords='data',textcoords='data',ha='left', arrowprops=arrowprops)
plt.gca().annotate('0.05', xy=(13.3, 0.003), xytext=(15.5,0.016), **kw)

plt.xlim(0,25)
plt.ylim(*y_range)

plt.xlabel('統計量 $t$')
plt.show()
```

実行結果は次のようになります。<br>![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/ff150697-a5cc-6ba0-d5bb-c976114629e7.png)


実際に、統計量 $t$ と、それに対応する $p$値 を求めます。

```py:検定統計量とp値の計算
t = ( ((f-fp)**2)/fp ).sum()        # 統計量 t
p = 1 - st.chi2.cdf(t,df=len(n)-2)  # 自由度6(8-2)のカイ2乗分布におけるtに対応するp値

print(f'統計量 t = {t:.2f}、p値 = {p:.2f} ',end='')
if p >= 0.05 :
  print('（ >= 0.05 ) \n帰無仮説（ポアソン分布に従う）は棄却されない')
else :
  print('（ < 0.05 ) \n帰無仮説（ポアソン分布に従う）は棄却される')
```

実行結果は次のようになります。

```
統計量 t = 3.22、p値 = 0.78 （ >= 0.05 ) 
帰無仮説（ポアソン分布に従う）は棄却されない
```