適合度の検定

# 適合度の検定

キーワード：名義尺度変数，$\chi^2$ 分布，漸近検定，一様性の検定

## 例題

Julia の乱数発生関数 `rand()` で，1 〜 $m(=10)$ の整数乱数を $n(=200)$ 個発生させた。1 から $m$ の整数が均等に発生されたといえるか。

整数 | 個数
:----:|:----:
1 | 24
2 | 25
3 | 17
4 | 14
5 | 14
6 | 22
7 | 26
8 | 21
9 | 18
10 | 19
合計 |200

常識的に考えてそれぞれの数値が 20 回程度は出てくると思うであろう。

所見では，4 とか 5 が 14 回しか出ていない。7 は 26 回も出るのか，1 とか 2 も多いしと思うであろ。

これを，統計学的に見るとどうだろうか...という話。

このデータは以下のように生成したものである。$n$，$m$ を変えれば，いろいろなデータを生成できる。


```julia
using Random; Random.seed!(123)
n = 200
m = 10
x = rand(1:m, n)
println(x)
```

    [6, 6, 9, 2, 6, 4, 1, 10, 6, 4, 6, 9, 1, 5, 4, 2, 8, 6, 7, 3, 3, 10, 9, 3, 4, 2, 1, 4, 5, 8, 3, 1, 6, 4, 2, 3, 4, 1, 1, 10, 4, 8, 4, 9, 3, 10, 9, 8, 2, 6, 9, 3, 9, 4, 7, 10, 8, 7, 9, 5, 5, 2, 6, 6, 5, 9, 6, 9, 6, 8, 1, 10, 9, 8, 8, 7, 5, 1, 7, 8, 10, 2, 2, 1, 1, 7, 7, 3, 7, 4, 8, 9, 7, 1, 2, 2, 8, 9, 2, 8, 6, 1, 6, 2, 7, 10, 3, 10, 7, 6, 3, 3, 8, 7, 8, 1, 5, 1, 1, 1, 10, 7, 7, 8, 8, 1, 7, 10, 1, 1, 6, 6, 7, 9, 2, 3, 10, 4, 5, 2, 3, 2, 7, 5, 1, 6, 2, 10, 9, 10, 5, 7, 7, 7, 5, 6, 8, 8, 10, 8, 2, 7, 2, 10, 5, 2, 1, 4, 7, 3, 3, 10, 2, 3, 2, 7, 6, 4, 1, 9, 6, 9, 2, 10, 8, 9, 10, 2, 7, 7, 2, 1, 1, 8, 7, 5, 6, 5, 2, 3]


集計すると以下のようになる。

度数分布は `FreqTables.freqtable()` で行う。


```julia
using FreqTables
f = freqtable(x)
```




    10-element Named Vector{Int64}
    Dim1  │ 
    ──────┼───
    1     │ 24
    2     │ 25
    3     │ 17
    4     │ 14
    5     │ 14
    6     │ 22
    7     │ 26
    8     │ 21
    9     │ 18
    10    │ 19



度数の合計は $m$ になる。


```julia
sum(f)
```




    200



度数分布図を描いてみると，かなり凸凹しているのが実感される。

これでも，1 〜 10 が均等に出現しているといえるのか?


```julia
using Plots
pyplot(grid=false, tick_direction=:out, label="");
bar(f, xticks=1:m, ylabel="個数")
```




    
![output_7_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/e141d26a-1300-8e90-b58f-b7e687ed2862.png)

    



## 検定手順

帰無仮説と対立仮説をたてる。

帰無仮説 $H_0$: 1 〜 $m$ の整数の出る確率はすべて等しい。 すべての $i$ において，$p_i = 1/m$

対立仮説 $H_1$: 1 〜 $m$ の整数の出る確率はすべて等しいのではない。どれかの $i$ で $p_i \neq 1/m$

対立仮説は，ある意味自分が興味を持っていて（極端にいえば）証明したいこと。今の場合でいえば，「ちょっと凸凹していて一様ではないんじゃないか?」ということ（作業仮説に似ている）。

帰無仮説はそれに対して，「マアマア，穏便に。そう騒ぎ立てなさるな。凸凹のように見えるが大した差ではないでしょう。」というようなこと。

対立仮説は複雑で扱いづらいが，帰無仮説は単純である。そこで，帰無仮説が正しいと仮定して，凸凹の程度を評価するために計算式により検定統計量を計算する。

観察された場合の検定統計量より極端な値を取る確率を求める。これを $p$ 値と呼ぶ。

観察された事象が起こりにくいのか，そうでもないのかの判断基準を有意水準 $\alpha$ という。普通は（慣例的に） $\alpha = 0.05$ とかが使われる。つまり，あることが起こる確率が 0.05 以下なら，そのようなことは起こりにくいと判断するということである。有意水準は検定を行う人（筆者）が決めるが，読者（第三者）が独自に設定して運用することもできる。つまり，筆者は有意水準 5% で有意といっているが，自分は有意水準 1% で判断する（その結果有意とは認めないよ）ということもありうる。このようなことができるためにも，検定結果は $p$ 値とともに提示されなければならないのである。

帰無仮説のもとで観察された事象が起こりうる確率（$p$ 値）が設定された有意水準より小さければ，その前提「帰無仮説が正しい」ということが疑わしくなる。つまり，「帰無仮説が正しいなら，このような観察事例は起こりにくい。しかし，実際に起きている。そもそも $p$ 値を計算したのは帰無仮説が正しいと仮定してのことである。つまり帰無仮説が間違っている。」と結論づける。この過程を「背理法」として説明する人がいるが，それは正しくはない。統計学検定は「帰無仮説が間違っている」と言い切ることができない。帰無仮説が正しいこともありうる。帰無仮説が正しいのに帰無仮説が間違っていると判断する確率は有意水準 $\alpha$ に等しい。この誤りのことを第一種の過誤 type I error または $\alpha$ errorと呼ぶ。

検定の結論は以下のようになる。

1. $p \lt \alpha$ なら，帰無仮説を棄却する。すなわち，対立仮説を採択する。
2. $p \gt \alpha$ なら，帰無仮説は棄却できない（対立仮説を採択するのではない）。

「等号はないのか?」あってもよいが，$p = \alpha$ の確率は 0 なので，実質的に考慮する必要はない。$p \leqq \alpha$ とすることが多いようではある。

さて，実際に検定を進めよう。

サンプルサイズ $n$，度数分布表の階級数 $m$，観察度数 $o_i$ とする。


```julia
n = 200;
m = 10;
o_i = [24, 25, 17, 14, 14, 22, 26, 21, 18, 19];
o_i |> println
```

    [24, 25, 17, 14, 14, 22, 26, 21, 18, 19]



```julia
p_i = repeat([1/m], inner=m);
p_i |> println
```

    [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]


期待度数 $e_i$ を求める。

$\hspace{2cm}e_i = n \ p_i,\ \ i = 1, 2, \dots, m$


```julia
e_i = n * p_i;
e_i |> println
```

    [20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0]


検定統計量 $\chi^2_o$ を求める。

$\hspace{2cm}\displaystyle \chi^2_o = \sum_{i = 1}^{m} \frac{(o_i - e_i)^2}{e_i}$


```julia
chisq_o = sum((o_i - e_i).^2 ./ e_i)
```




    8.4



$\chi^2_o$ は 自由度 $df = m - 1$ の $\chi^2$ 分布にしたがう。


```julia
df = m - 1
```




    9



$p$ 値は以下のようにして求める。

`Distributions.ccdf(Chisq(df), chisq_o)` は，自由度が `Chisq(df)` の $\chi^2$ 分布の上側確率 $P(x \geqq \chi^2_o)$ を求める関数である。 


```julia
using Distributions
p_value = ccdf(Chisq(df), chisq_o)
```




    0.4943916864821227



別のデータも評価できるように，図示する関数を定義する。

右側の水色で塗りつぶした部分の面積（確率）が $p$ 値である。




```julia
using Plots
pyplot()
x = range(0, 30, step=0.01);
y = pdf.(Chisq(df), x);
plot(x, y, xlabel=raw"$\chi^2$", ylabel=raw"$f(\chi^2, m-1)$")
x2 = range(chisq_o, 30, step=0.01);
y2 = pdf.(Chisq(df), x2);
plot!(vcat(x2, 30, chisq_o, chisq_o...), vcat(y2, 0, 0, y2[1]...),
    seriestype=:shape, linecolor=1, fillcolor=:tomato)
hline!([0], linecolor=1)
plot!([20, 12], [0.06, 0.025], arrow=1)
annotate!(20, 0.06, raw" $p$ value", :left)
```




![output_20_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/c4f8ef3b-a01c-9604-05d9-2d437c92b603.png)




帰無仮説が正しい場合には，定義上からいって検定統計量 $\chi^2_o$ はあまり大きい値にはならない。したがって，$p$ 値はあまり小さな値にはならない。つまり，そんな極端なことは起きにくいということである。

今の場合は，$p \gt 0.494 \gt \alpha = 0.05$ つまり $p \gt \alpha$ なので，帰無仮説が正しいとした場合には，このような事が起きるのはそう珍しいことではない。ということで，有意水準 5% で帰無仮説は棄却できないと判断する。

すなわち，「1 〜 10 の数値の出現比率に差があるとはいえない」という結論になる。

### 汎用プログラム

同じような問題をたくさん解決するための汎用関数を定義しよう。

引数に度数分布ベクトルをあたえると，$\chi^2_o$ と自由度，$p$ 値を返す。


```julia
using Distributions
function fitnessofgoodnesstest(o)
    m = length(o)
    p_i = repeat([1/m], inner=m)
    e_i = n * p_i
    chisq_o = sum((o_i - e_i).^2 ./ e_i)
    df = m - 1
    p_value = ccdf(Chisq(df), chisq_o)
    chisq_o, df, p_value
end
```




    fitnessofgoodnesstest (generic function with 1 method)



例題を検定してみよう。正しい結果が出ているようだ。


```julia
fitnessofgoodnesstest([24, 25, 17, 14, 14, 22, 26, 21, 18, 19])
```




    (8.4, 9, 0.4943916864821227)



### 応用例

20000 個の整数乱数を発生させて，同じように検定を行うとどうなるか？


```julia
using Random; Random.seed!(123)
n = 20000
m = 10
x = rand(1:10, n)
o_i = freqtable(x)
```




    10-element Named Vector{Int64}
    Dim1  │ 
    ──────┼─────
    1     │ 1954
    2     │ 1990
    3     │ 1901
    4     │ 1950
    5     │ 1954
    6     │ 2011
    7     │ 2106
    8     │ 2061
    9     │ 2051
    10    │ 2022




```julia
pyplot(grid=false, tick_direction=:out, label="");
bar(o_i, xticks=1:10, ylabel="個数")
```




    
![output_28_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/dcb7410a-b9ce-94c4-88f5-7a9b466122a0.png)

    



分布状況を見ると， $m=200$ だったときよりもむしろ均一に見える。

しかし，検定をしてみると。


```julia
chisq_o, df, p_value = fitnessofgoodnesstest(o_i)
```




    (17.398, 9, 0.04283610712658917)



のようになり，有意水準 5％ で帰無仮説は棄却される（つまり，一様分布ではない）ということになる。


```julia
x = range(0, 30, step=0.01);
y = pdf.(Chisq(df), x);
plot(x, y, xlabel=raw"$\chi^2$", ylabel=raw"$f(\chi^2, m-1)$")
x2 = range(chisq_o, 30, step=0.01);
y2 = pdf.(Chisq(df), x2);
plot!(vcat(x2, 30, chisq_o, chisq_o...), vcat(y2, 0, 0, y2[1]...),
    seriestype=:shape, linecolor=1, fillcolor=:tomato)
hline!([0], linecolor=1)
plot!([22, 19], [0.03, 0.005], arrow=1)
annotate!(22, 0.03, raw" $p$ value", :left)
```




    
![output_32_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/e10b9958-a265-b7ba-c272-4c2f5b566428.png)

    



なぜか？

すべての検定は，

1. 帰無仮説との乖離の程度が大きい場合
2. サンプルサイズが大きい場合

に有意な結果になりがちという性質を持つ。

 1. の場合は望むところであるが， 2. の場合は熟慮が必要だ（本当に乖離が大きいのか，それとも単にサンプルサイズが大きいせいなのか）。
 
サンプルサイズが大きいときには，検定結果が有意になっても，嬉しがらないように注意。