統計: ピアソンのカイ二乗検定で標本が理論分布と適合しているか調べる

例えば、ある六面ダイス (サイコロ) に歪みがないことを調べたいとする。もしサイコロに歪みが無いなら、出る目の理論的な度数分布はどれも $\frac{1}{6}$ となるはず。しかし、サイコロの出る目は無限母集団なので、実際にすべてのパターンを試して確認することができない。つまり、全数調査は不可能ということ。そのため、歪みがあるか否かは実際に何度かそのサイコロを振った有限な結果から推測する必要がある。これはつまり、標本から母集団の分布を調べる推測統計になる。

上記のようなシチュエーションでは、今回紹介するピアソンのカイ二乗検定という方法が使える。ピアソンのカイ二乗検定で適合度を調べると、実際に振ってみた結果から母集団が理論分布となっているか否かが判断できる。この検定はノンパラメトリックなので、特定の分布の仕方には依存しないところが便利に使える。

ピアソンのカイ二乗検定の公式 (適合度)

まず、ピアソンのカイ二乗検定では、カイ二乗値という統計量を計算することになる。これは、次のような公式で計算する。

$\chi^{2} = \displaystyle \sum_{i = 1}^{n} \frac{(O_i - E_i)^{2}}{E_i}$

上記で $O$ は実際に観測した値、つまり標本を表している。また、 $E$ は期待値、つまり母集団を表している。変数となっている ${}_i$ は分類、つまりサイコロでいえば出目に相当する。

数式の意味をざっくり説明すると、期待値と観測した値の差を取ってそれを二乗している。ここで、二乗するので得られる値は必ず正の値になることが分かる。そして、その値をまた期待値で割っている。これは、期待値とのズレを二乗した値と、期待値との比率を計算していることになる。そして、その比率をすべての分類ごとに足し合わせたのがカイ二乗値となる。

で、得られた値はどうするかというと、自由度と有意水準にもとづいてカイ二乗分布表と比べる。そして、カイ二乗分布表の値よりも小さければ母集団が理論上の分布に沿っていると判断できる。ただ、ここは今いっぺんに説明するよりも、あとから詳しく書いた方が分かりやすいと思う。

実際に計算してみよう

それでは、実際にサイコロのカイ二乗値を計算してみることにしよう。ただ、本物のサイコロを用意して振るのがめんどくさかったので計算機のシミュレーションで済ませてしまうことにする。これも、モンテカルロ法という名前のついた実験方法のひとつだ。

モンテカルロ法を実施する環境としては Python を使うことにしよう。

$ python --version
Python 3.5.2

まずは Python のインタプリタを起動する。

$ python

REPL が起動したら、サイコロに相当する関数 dice() を定義する。

>>> import random
>>> def dice():
...     return random.randint(1, 6)
...

この関数を実行するとサイコロの出目が 1 から 6 の間でランダムに得られる。

>>> dice()
3
>>> dice()
5
>>> dice()
5

果たして、このサイコロは歪んでいるのか、それとも歪んでいないのか。今回は、それをカイ二乗検定で確かめたい。

ダイスを所定の回数だけ振る

まず、カイ二乗検定を使うときは、次の条件を満たすようにすべきらしい。

$np_i \ge 10$

上記の $n$ は標本数で $p_i$ は分類の理論的な出る確率を表している。

今回、理論的な確率は全ての出目で $\frac{1}{6}$ となることを想定している。

$p_i = \frac{1}{6}$

そのため $n$ について解くと

$n \frac{1}{6} \ge 10$

から

$n \ge 10 \times 6 = 60$

となるため、標本は最低でも 60 以上ほしい。これはつまり、カイ二乗値を計算するのに最低でも 60 回はサイコロを振る必要があるということ。

そこで、ひとまず多めに 120 回振ることにした。 Python には collections.Counter という種別ごとに登場する回数を数えるためのクラスがあるので、それを使う。

>>> from collections import Counter
>>> c = Counter([dice() for _ in range(120)])
>>> c
Counter({1: 22, 3: 22, 2: 21, 6: 21, 5: 20, 4: 14})

上記はサイコロの出目が 1 のときは 22 回あった、というように読める。つまり 1 が出た回数が一番多くて、4 が出た回数が一番少なかったということ。

上記は出た回数が多い順にソートされていて、ちょっと見づらい。そこで、出目ごとにソートすると、こんな感じ。

>>> sorted(c.items(), key=lambda t: t[0])
[(1, 22), (2, 21), (3, 22), (4, 14), (5, 20), (6, 21)]

カイ二乗値を計算する

標本が得られたのでカイ二乗値を計算する準備が整った。先ほどの公式をもう一度見てみよう。

$\chi^{2} = \displaystyle \sum_{i = 1}^{n} \frac{(O_i - E_i)^{2}}{E_i}$

上記で $O_i$ が観測値なので、これはさっきサイコロを振って得られた。

そして $O_i$ は期待値となっている。これは $np_i$ で計算できるので、こうなる。 $n$ が標本数で $p_i$ が確率だった。確率は、すべての出目で $\frac{1}{6}$ となるのは前述した通り。

$E_i = np_i = 120 \times \frac{1}{6} = 20$

つまり期待値として使う値はすべての出目で 20 にすれば良い。

ということで、実際に数式に値を代入してみると、こうなる。

$\chi^{2} = \displaystyle \sum_{i = 1}^{n} \frac{(O_i - E_i)^{2}}{E_i} = \frac{(22 - 20)^{2}}{20} + \frac{(21 - 20)^{2}}{20} + \frac{(22 - 20)^{2}}{20} + \frac{(14 - 20)^{2}}{20} + \frac{(20 - 20)^{2}}{20} + \frac{(21 - 20)^{2}}{20}$

上記で使った値として、それぞれの出目の回数を見に戻るのが大変だろうから再掲しておく。

>>> o = [i[1] for i in sorted(c.items(), key=lambda t: t[0])]
>>> o
[22, 21, 22, 14, 20, 21]

これは、計算してみると次のようになる。

$\chi^{2} = \frac{4 + 1 + 4 + 36 + 0 + 1}{20} = \frac{46}{20} = \frac{23}{10} = 2.3$

Python で計算するまでもないけど、一応確かめておこう。すでに標本の入った変数 o はあるので、あとは期待値の入った変数 e を用意しよう。これは、すべてに 20 が入った要素数 6 のリストになる。

>>> e = [20 for _ in range(6)]
>>> e
[20, 20, 20, 20, 20, 20]

順を追って計算していこう。まずは標本から期待値を引いて二乗しよう。これで公式の分母部分が計算できる。

>>> [(oi - ei) ** 2 for oi, ei in zip(o, e)]
[4, 1, 4, 36, 0, 1]

それらを足すと、こうなる。

>>> sum([(oi - ei) ** 2 for oi, ei in zip(o, e)])
46

そして、それを期待値で割る。

>>> sum([(oi - ei) ** 2 for oi, ei in zip(o, e)]) / 20
2.3

ちゃんと 2.3 になった。

標本の自由度を計算する

次に、カイ二乗分布表を参照するために標本の自由度を計算しよう。自由度というのは、標本でいくつまで種別の値が埋まると、種別全体の値が確定するかという値になっている。

今回の例なら標本数が 120 というのは、あらかじめ分かっている。そして種別は 6 種類あって、各種別が何回登場するかは 5 種類の値が分かれば確定する。これはつまり 1 から 5 までの出目の回数が分かってしまえば、6 が何回出たかは自然と分かるということ。ということで、今回使う自由度 $df$ はサイコロの面の数である 6 から 1 を引いた値で 5 になる。

$df = k = 1 = 6 - 1 = 5$

ただし、これには注意点があって、使う自由度の数は母分布によって異なる。サイコロであれば全ての種別が同じ確率で出るので一様分布になる。一様分布であれば自由度は種別の数から 1 を引いたもので良い。これがピアソン分布や二項分布であれば 2 を、正規分布であれば 3 を引いて使う。

カイ二乗分布表を参照する

自由度が求まったので、それにもとづいてカイ二乗分布表を参照しよう。通常であれば、あらかじめ求めてあるものを用いる。例えば、次のサイトなどにあるようなもの。

χ2分布表

なぜ、ここでカイ二乗分布表というものが出てくるのか疑問に思うかもしれない。実は、先ほど計算したカイ二乗値というものは、自由度によってどのような値が出るのかが既に分かっている。これがカイ二乗分布表というもので、ようするにどんな値がどれくらいの確率で出るか書いてある。

表の見方としては、まず先ほど求めた自由度 (df) にもとづいて 5 の行を参照する。ひとつの行には色々な数値が並んでいるが、これは先ほど計算したのと同じカイ二乗値になっている。それぞれのカイ二乗値は何が異なるかというと、その値になる確率が異なっている。

上にある列を見ると有意確率が書いてある。例えば有意確率 0.99 の列で自由度が 5 の項目を見ると 0.55 となっている。これは、自由度 5 のカイ二乗値は 99% の確率で 0.55 以上となることを表している。

有意水準を決める

カイ二乗検定というのは、仮説検定の一種となっている。仮説検定の考え方では、ある統計量を計算したとき、その値がどれくらいの確率で出るのかを考える。

例えば、その値が出る確率が 5% 未満や 1% 未満であれば、とても起こりづらいことは明らかだ。そんなとき、仮説検定では「珍しいことが起こった」のではなく「本来起こるべきでないことが起きた」と考える。例えば、サイコロでいえば仮に歪みがなかったとしても偶然偏った値になるのはありえることだろう。しかし、そんなとき仮説検定の考え方では「偶然偏った値になった」のではなく「サイコロ自体に歪みがあった」という結論にする。

もちろん、この考え方は 5% 未満や 1% 未満など、ある水準で間違った結論を導いてしまう恐れはある。これを統計では過誤という用語で呼ぶ。しかし、その過誤が特定の確率で起こることは折り込んだ上で、そのように考えるらしい。この 5% や 1% といった基準のことを有意水準と呼んで、記号では $\alpha$ が使われることが多い。仮説検定では、この有意水準にもとづいて統計量を考える。

有意水準は一般的に 5% (0.05) もしくは 1% (0.01) を採用することが多い。今回は試しに 5% (0.05) を使って考えてみることにしよう。これは 5% の確率で過誤が起こりうる、という意味でもある。

カイ二乗値を比べる

それでは、カイ二乗分布表をもう一度見てみよう。有意確率が 5% (0.05) で自由度が 5 の項目の値は 11.07 となっている。

カイ二乗分布表の値は、それ以上の値が出る確率がどれだけあるか、というものだった。つまり、自由度が 5 のカイ二乗値が 11.07 以上になる確率は 5% ということだ。ようするに計算したカイ二乗値が 11.07 以上になると、サイコロが歪んでいると判断できるわけ。

では、先ほど計算したカイ二乗値はいくつだっただろうか。その値は 2.3 で、あきらかに 11.07 と比べると小さいことが分かる。

$\chi^{2} = 2.3 \lt 11.07$

つまり、このカイ二乗値は全然ありうる値で、サイコロは歪んでいないという結論に至った。

SciPy で計算してみる

ちなみに、上記の計算は手作業が多くてなかなか面倒くさい感じだと思う。特にカイ二乗分布表を参照することなんかは間違えやすそうだ。そこで Python の SciPy というパッケージを使った方法も紹介しておくことにする。

SciPy はサードパーティ製のパッケージなので、まずは pip を使ってインストールしておこう。

$ pip install scipy

そして、改めて Python の REPL を起動する。

$ python

先ほどの出目と同じ状況を用意する。変数 o が標本で変数 e が期待値となる。

>>> o = [22, 21, 22, 14, 20, 21]
>>> e = [20 for _ in range(6)]

続いてカイ二乗検定をする関数を SciPy からインポートしよう。これには scipy.stats.chisquare() を使う。

>>> from scipy.stats import chisquare

あとは、この関数に先ほどの変数 o と e を渡すだけだ。

>>> chisquare(o, e)
Power_divergenceResult(statistic=2.2999999999999998, pvalue=0.80626686988512852)

するとカイ二乗値が 2.23 で、p-値というのが 0.80 と計算されているのが分かる。

p-値というのは、それがどれくらい出やすいのかを表した値になっている。上記であれば 0.80 なので、この値が出る確率は 80% となって全然珍しくないことが分かる。

サイコロが歪んでいるかは、有意水準よりも小さいかで判断する。例えば有意水準が 5% であれば次のようにすれば良いというわけ。

>>> chisquare(o, e).pvalue < 0.05
False

5% の確率で間違っている恐れはあるものの、このサイコロは歪んでいないことが分かった。

コクがありそうなサイコロを使ってみる

先ほど試したのは一様乱数にもとづいたサイコロだった。次は、あえて歪んだサイコロを振って試してみよう。

今度は、次のようにして正規乱数にもとづいたサイコロを使ってみる。このサイコロは、小さな値ほど出やすいことになる。

>>> import random
>>> def dice():
...     return int(abs(random.normalvariate(0, 1) * 6)) % 6 + 1
...

今度は SciPy を使って一気に計算させるのでサイコロを振る数も 1200 回に増やしてみる。

>>> n = 1200
>>> c = Counter([dice() for _ in range(n)])
>>> o = [i[1] for i in sorted(c.items(), key=lambda t: t[0])]
>>> e = [n / 6 for _ in range(6)]

この結果、カイ二乗値は 69.33 が得られた。 p-値は $1.4e^{-13}$ というとんでもなく小さな値になっている。

>>> from scipy.stats import chisquare
>>> chisquare(o, e)
Power_divergenceResult(statistic=69.330000000000013, pvalue=1.4126539924787112e-13)

これを有意水準 5% で検定すると、このサイコロは歪んでいることが分かった。

>>> chisquare(o, e).pvalue < 0.05
True

まとめ

カイ二乗検定を使うと理論分布と適合しているかを調べることができる

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る

CUBE SUGAR CONTAINER

技術系のこと書きます。