俺たちは何回コインを投げるべきなのか

表と裏が同じ確率で出ることが期待されているコインがあり、コインにイカサマがないことを確認したいとする。なんやかんやあって、統計的手法で、不正の有無(つまり、コインの表が出る確率が50%か否か)に決着をつけることになった場合、コインを何回投げれば十分といえるか


統計的手法にも色々あるかもしれないけど、不正の有無を統計的にチェックしろと言われた場合、オーソドックスに使われるのは仮説検定と思う。表が出る回数は二項分布で表されるが、コインを投げる回数Nが十分大きい時、これは正規分布で近似できるようになって、実際に表が出た回数をx、有意水準は、例えば95%とすると
\frac{|x - Np|}{\sqrt{Np(1-p)}} > 1.96
であれば、コインにイカサマがないという仮説は棄却される(但し、p=1/2で、不正がない時、表が出る確率)。二項分布を正規分布で近似した時の分散σは
\sigma^2 = Np(1-p)
で、おおまかに、イカサマがない時の期待値Npから2σ以上外れていればアウトという判定基準になっている(偏差値でいうと、70以上か30以下)。イカサマを疑われている側からすると、イカサマしてない場合でも、冤罪をかけられる確率が5%あるということなので、もっと厳しくしたいと思うかもしれないが、そこらへんは話し合いで同意が取れたとする。もし、表が出る確率が本当に寸分の狂いもなく1/2だったら、(Nが十分大きい時には)コインを投げる回数が100回だろうが5000兆回だろうが、冤罪発生率が5%であることは変わらない。どんなに公正なコインでも、たまたま偏った出目になる可能性は0でないので、残念ながら、これを0%にすることはできない


逆に、本当にイカサマがあって、表が出る確率が60%くらいあった場合、コインをN=5回しか投げないで、有意差がないという結論になっても、明らかに、Nが小さすぎるので、もっとNを増やす必要がある。こういうのは検出力(statistical power)が足りないと言われる。また、極端な話、表が出る確率が90%であれば、少し試してみれば不正があることを見抜けそうであるけど、表が出る確率が60%などの場合は、それより多くのNを必要とすると予想される。確率的なものだから、表が出る確率が0%でない限り、たまたま全部表が出て、イカサマがあるという結論になってしまう可能性はある。けど、例えば、N=10として、表が出る確率が90%あるコインを投げた時、全部表になる確率は35%くらいあるのに対して、表が出る確率が60%のコインを投げて全部表になる確率は0.6%くらいしかない。


数学ばかり見てないで、現実を見ると、コインは表裏対称でないことが多いし、そもそも、人間が肉眼で区別できないほど表裏対称だったら、コインの表が出たのか裏が出たのか判定できない。表面にマーキングしたり汚れがある場合など、人間にとっては微細に思える差であっても、表が出る確率が、1/2から僅かにずれるには十分と思われる。何にせよ、コインが表が出る確率が完全に1/2であることは多分期待できない。仮説検定は、差があるかないかしか調べないので、表が出る確率が50.0001%とかで、わずかに裏より表が出やすいだけという場合でも、Nが大きければ、有意差ありという結果になる可能性は非常に高くなりうる。実際、最初に書いた仮説検定の式を少し変形すると
 p - 1.96 \frac{\sqrt{p(1-p)}}{\sqrt{N}} < \frac{x}{N} < p + 1.96 \frac{\sqrt{p(1-p)}}{\sqrt{N}}
のようになって、表が出た割合x/Nが許される範囲は、Nが大きくなるに従って狭くなっていく。一方、Nが大きくなれば、x/Nは、大数の法則によって、表が出る真の確率に収束していく。結果として、表が出る確率が寸分の狂いもなく1/2でない限り、有意差ありという結論になる可能性は、Nが大きくなると、どんどん高まっていくことになる。しかし、表が出る確率が、実は50.0001%で、統計的に有意な差があると言われても、常識的な感覚からすると、意味のある差とは思えない


ということを考えると、不正を検出するために必要なNの大きさは、どのくらいの不公平まで許容範囲と考えるかに依存する。つまり、表が出る確率が50±ε(%)の範囲外であればよくないコインであるが、範囲内であれば問題ない水準と考えられるεの大きさを事前に決めておく必要がある。このεの大きさは、有意水準αと同じく客観的な決め方はない。


更に冤罪発生率を0%にできないのと同様、コインに不正がある場合でも、たまたま、表と裏が同じ回数ずつ出てしまう可能性は0にはならない。つまり、不正がある時に不正を100%検出することはできない。とはいえ、Nを大きく取れば、表が出る確率が真に1/2でない限り、このようなことが起きる確率は減っていくと考えられる。というわけで、不正がある時に、正しく不正が検出できる確率を、有意水準と同様に設定する必要がある。これが検出力の定義で、1-βで表されることが多い。というわけで、α、β、εの3つのパラメータを設定すると、不正検出のためにコインを投げる回数Nを決められそうな気がする。


もう少し定量的に評価することを考える。さっきの式に戻ると、有意水準α=0.05と設定して、コインをN回投げた時、表が出た回数xが
 p - 1.96 \frac{\sqrt{p(1-p)}}{\sqrt{N}} < \frac{x}{N} < p + 1.96 \frac{\sqrt{p(1-p)}}{\sqrt{N}}
を満たせば、コインは公平であると帰結される(前回同様、p=1/2)。

コインに不正があって、表が出る確率がp+εであったとする。このとき、コインが表が出る回数は、やはり二項分布で記述され、上の区間に入る確率を計算することができる。この確率は、表が出る確率がp+εであるにも関わらず、不正はないと結論づけてしまう確率とみなせる。この二項分布は、Nが十分大きい時、平均N(p+ε)、分散N(p+ε)(1-p-ε)の正規分布で記述されることに注意すれば、不正検出ミスの確率は簡単に計算できる。表が出る確率が、例えばp+2εだとか、p+2.4εだというケースもありえるけど、検出力が最も低くなるのは、本当の表が出る確率がp+εか、p-εの場合なので、この2点のみチェックすればいい。更に、p+εの場合と、p-εの場合で、結果は同じことになるので、片方だけ見ればいい。この理由から、以下、εは正としておく


自然言語で書くと状況は複雑だけど、数式で書くと、表が出る真の確率がp+εの時に
\int_{p-u}^{p+u} \frac{1}{\sqrt{2\pi\sigma^2}} e^{-(t-\mu)^2/2\sigma^2} dt < \beta
という条件を満たすということになる。但し、p=1/2で
u = 1.96 \frac{\sqrt{p(1-p)}}{\sqrt{N}}
\mu = p+\epsilon
\sigma^2 = (p+\epsilon)(1-p-\epsilon)/N
とする(積分変数のtはコインを投げて表が出た"割合"を表すので、平均と分散を、それぞれN、N^2で割っていることに注意)。で、積分の変数変換をすると
\frac{1}{\sqrt{\pi}} \int_{(-u-\epsilon)/\sqrt{2\sigma^2}}^{(u-\epsilon)/\sqrt{2\sigma^2}}e^{-x^2} dx
のようになる。erf関数
erf(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt
を使うと、最終的に
\frac{1}{2} \left( erf((u-\epsilon)/\sqrt{2\sigma^2}) - erf((-u-\epsilon)/\sqrt{2\sigma^2}) \right) < \beta
となる。とりあえず左辺を評価するためにεを決める必要がある。まぁ、ε=0.01としておく。


数値的に、いくつか計算してみる。pythonで以下のようなコードを実行すると

import math

for N in [10000,15000,20000,25000,30000,35000,40000,50000]:
   p,eps = 0.5,0.01
   u = 1.96/math.sqrt(N/p/(1-p))
   sigma2 = (p+eps)*(1-p-eps)/N
   print(N , 1.0-(math.erf((u-eps)/(math.sqrt(2*sigma2))) - math.erf((-u-eps)/(math.sqrt(2*sigma2))))*0.5)

以下のような結果を得た。

10000 0.5159939775494847
15000 0.6877923080943419
20000 0.8074680951250419
25000 0.8854187382561318
30000 0.9337611517760266
35000 0.9626265172738802
40000 0.9793451544602965
50000 0.9940083977523118

大体
N=20000回で、検出力0.8ちょっと
N=25000回で、検出力88.5くらい
N=50000回で、検出力0.994...
などとなる。つまり、20000回くらいやれば、コインが表が出る確率が51%の時、80%以上の確率で有意差がある(有意水準は5%)という検定結果を得られる計算。まぁ、数万回くらいはコインを投げることになりそう。一日が86400秒なことを考えると、実行するには、くそだるい数字である。


検出力は、0.8とか0.9にすることが多いらしいが、αやε同様、客観的に決める基準はないので慣習的なものに過ぎない。また、有意水準を厳しくすると、同じ検出力を得るには、より多くの試行回数を必要とする。例えば、3σを要求すると

10000 0.15860713527390746
15000 0.29094698846888134
20000 0.43187317605341713
25000 0.5644691796250857
30000 0.6787457862331097
35000 0.7708974857449786
40000 0.841393149895479
50000 0.9295476650393496

のようになった。εが厳しすぎる気もするので、5%くらいまで緩めると、数千回程度でも十分な検出力を期待できるようになる。大雑把に言って、有意区間の幅は、Nの平方根に反比例するので、(αやβを変えずに)εを1/10にしたければ、Nを100倍にする必要があり、逆にεを10倍にしてよければ、Nは1/100程度まで小さくできる。このおかげで、今の場合、(例えば、濡れ衣を着せてやろうと思って)Nを大きくして小さな差でも有意になるようにしようとしても、実際に実行するのはかなり難しくなる。例えば、εを0.001%にして(上で計算したεが1%の時と)同じ有意水準、検出力で検定しようと思えば、100万倍ものデータ量が必要となる。



元ネタ(検定力と検出力は同じもの):
【翻訳】ダメな統計学 (3) 検定力と検定力の足りない統計
http://id.fnshr.info/2014/12/17/stats-done-wrong-03/

これは本にもなってる。Webページには、コインの例が載っているのだけど、上に書いたような計算は載っていなかった(数式を載せるような本ではないせいだろうと思う)ので、計算してみた。適切な統計学の教科書には例題として載っていそうな気もするけど、私が大学で受けた授業では、検出力/検定力とか聞かなかった気がする。



演習問題:各面が同じ確率で出ることが期待されている(6面)サイコロに不正がないかを統計的に調べたいとすると、俺たちは一体、何回サイコロを投げるべきか?

方針:教科書的には、カイ二乗(適合度)検定を行うことが第一歩となるだろうと思う(他の検定手法もあるかもしれないが)。

次に、コインの場合のεに相当する量を何らかの形で定義する必要がある。これに正しい方法というのがあるわけではないけど、サイコロの各面の期待される出現確率をp_i(i=1,...6)として、実際の出現確率をp_i'とする時、次のような量を定義すると便利
w = \sqrt{ \sum_{i=1}^6 \frac{(p_i - p_i')^2}{p_i} }
各面が同じ確率出ることを期待されているという条件は、勿論
p_i = \frac{1}{6}
を意味する。wは効果量(effect size)と呼ばれる。例えば、1〜3が出る確率が1/6+1/100で、4〜6が出る確率は1/6-1/100である(各面1%くらいの誤差)とすると、wは0.06となる。実際の研究では、wは、0.1,0.3,0.5などの値が選ばれることが多いらしい


効果量の定義の仕方は、検定の手法ごとに異なりうるし、検定の手法を決めても、必ずしも標準的なものが一つあるというわけでもないらしいので、Cohen's wという名前で呼ばれることもあるらしい
Cohen's w
https://en.wikipedia.org/wiki/Effect_size#Cohen.27s_w


それで、サイコロを投げる回数をNとして、wを使うと、カイ二乗値は、"平均的には"
\chi^2 = \sum_{i=1}^6 \frac{(O_i - E_i)^2}{E_i} = N w^2
となり、Nに比例して大きくなることが分かる(w=0でも、期待値は自由度kに等しくなるので、本当の意味での平均ではない)。カイ二乗値は、帰無仮説が正しい時(つまり、p_i=p_i'が成立する時)には、近似的にカイ二乗分布に従う(今の場合、自由度は5)というのが、カイ二乗検定の要旨だったけど、p_iとp_i'がずれている時には、非心カイ二乗分布で近似されることは直感的に分かる。検出力は、この非心カイ二乗分布積分して計算できるけど、非心度λを知る必要がある。ところで、非心カイ二乗分布の期待値が自由度kとλの和で書けることに注意すれば、λと、上で計算したカイ二乗値の値Nw^2とが関係することは意外ではない。etc. 具体的な計算は面倒なので省略