歪んだ硬貨があって,投げると確率 p で表(おもて)が出て,確率 1 - p で裏が出るとします。この硬貨を n 回投げて,表が r 回出る確率を求めてください。
これが2項分布(binomial distribution)の問題です。答えは
です。nCr は n 個から r 個を選ぶ組合せ(combination)の数で,英語ではよく n choose r と読みます。階乗(!)を使えば
と表せます。
Rでは,例えば 10C3 は choose(10,3)
と打てば求められます。
> choose(10,3) [1] 120
投げると確率 0.4 で表が出る硬貨を10回投げて表が3回出る確率
10C3 0.43 0.67
は dbinom(3,10,0.4) で求められます。
> dbinom(3,10,0.4) [1] 0.2149908 > choose(10,3)*0.4^3*0.6^7 [1] 0.2149908
投げると確率 0.5 で表が出る硬貨を10回投げて表が0〜10枚出る確率を全部出力するには次のようにします。
> dbinom(0:10,10,0.5) [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 [5] 0.2050781250 0.2460937500 0.2050781250 0.1171875000 [9] 0.0439453125 0.0097656250 0.0009765625
これをグラフにするには次のように打ち込みます。
barplot(dbinom(0:10,10,0.5),names.arg=0:10)

ただし,デフォルトではy軸の文字が横向きになります。上の図のようにするには,あらかじめ
par(las=1)
と打ち込んでおきます。par() は図を描く際のオプション(parameter)を指定するコマンドです。las は label style の意味です。詳しい説明はRに
?par
と打ち込んで出るヘルプを見てください。
これを順に加えていった累積確率は次のようにして求められます。
> pbinom(0:10,10,0.5) [1] 0.0009765625 0.0107421875 0.0546875000 0.1718750000 [5] 0.3769531250 0.6230468750 0.8281250000 0.9453125000 [9] 0.9892578125 0.9990234375 1.0000000000
以下では「p 値」という言葉が別の意味で登場しますので,2項分布のパラメータとしての確率は以下では θ という文字で表すことにします。
硬貨を10枚投げて表が2枚しか出ませんでした。この硬貨は歪んでいるでしょうか。
ある母集団から10人をランダムに選んで聞いたところ,賛成2人,反対8人でした。母集団全体でも反対のほうが多いと言えるでしょうか。
この問いについて考えるために,仮に硬貨は歪んでいない(あるいは母集団全体では賛否が等しい)というモデル(帰無仮説,null hypothesis)を立てます。そして,この帰無仮説が正しかった場合に,実際に観測された以上の外れ方(2:8,1:9,0:10,そして通常はさらにそれをひっくり返した8:2,9:1,10:0)が生じる確率の合計を求めます。この確率の合計を p 値(ピーち,p-value)といいます。p 値が非常に小さければ,実際に起きた事象はこのモデルでは説明しにくいので,たぶん硬貨は歪んでいる(あるいは賛否は等しくない)と推測します。p 値が大きければ,これだけのデータでは何も言えないということがわかるだけです。
p 値が大きいか小さいかの境界(有意水準)を仮に 0.05 として,p ≦ 0.05 であれば帰無仮説からの外れが「統計的に有意」(statistically significant)である,あるいは「帰無仮説は棄却(reject)された」ということがあります。0.05 という値に特に意味はありませんが,伝統的によく使われています(物理学ではもっともっと厳しい条件を課します)。
さて,硬貨を投げて表の出る枚数の分布は2項分布と考えられますので,表も裏も 1/2 の確率で出るとすれば,表が r 枚出る確率は 10Cr (1/2)10 です。表裏が 0:10,1:9,2:8,8:2,9:1,10:0 である確率はそれぞれ
> dbinom(c(0,1,2,8,9,10),10,0.5) [1] 0.0009765625 0.0097656250 0.0439453125 0.0439453125 [5] 0.0097656250 0.0009765625
で,この合計,すなわち p 値は
> sum(dbinom(c(0,1,2,8,9,10),10,0.5)) [1] 0.109375
になります。同じことが
> pbinom(2,10,0.5)*2 [1] 0.109375
でも求められます。また,後で詳しく述べますが,binom.test()
という関数でも2項検定ができます。
> binom.test(2,10,0.5)
Exact binomial test
data: 2 and 10
number of successes = 2, number of trials = 10, p-value = 0.1094
...
したがって,有意水準を 0.05 とすれば,表裏の差は統計的に有意ではありませんし,アンケートであればこんなに少人数の結果から「賛成が少ない」という結論を導いてはいけないということになります。
表が1枚(賛成が1人)なら,p 値は 0.02 ほどになり,水準 0.05 で有意になります。
これが,フィッシャー(R. A. Fisher,1890〜1962年)が「有意性の検定」(tests of significance,significance tests)と呼んだ方法の考え方です。
帰無仮説は現象の理解を助けるために設定した一つのモデルです。そのモデルと現実のデータとの整合性の度合を確率のことばで表したのが p 値です。p 値が非常に小さければ,モデルとデータは両立しにくいので,より良いモデルを考えなければならないことになります。
p 値は「帰無仮説が正しい確率」ではありません。どんな硬貨でも少しは歪んでいます。まったく賛否の割合が等しいことなどありえません。これらの場合,帰無仮説が正しい確率(というものが定義されるなら)は 0 であるというのが適当です。
ただし,Bayes(ベイズ)流(Bayesian,ベイジアン)の立場では,確率の概念を拡張して,仮説が正しい確率を求めようとします。
p 値は,標本の大きさ(調べた個数)によって大きく変化します。例えば歪んだ硬貨を調べる場合,投げた回数が非常に多ければ,検定の「感度」がうんと高くなるので,実用上無視できるほどのちょっとした歪みでも p 値が非常に小さくなることがあります。逆に,かなり違いがありそうでも,投げた回数が少なければ p 値はなかなか小さくなりません。いずれにしても,p 値が十分小さくなければ(帰無仮説が棄却されなければ)表が出やすいとか賛成が多いとか言うことはできません。
p 値が(有意水準より)小さくならなかったとき「帰無仮説を採択(accept)する」ということがありますが,これは誤解しやすい言葉です。「帰無仮説を棄却できない」あるいは「わからない」というほうが正しいでしょう。差がなかったのではなく,測定の感度が悪くて差が見つけられなかっただけなのです。
フィッシャーの「有意性の検定」は,後にネイマン(Jerzy Neyman, 1894〜1981)とピアソン(Egon Pearson, 1895〜1980)の「統計的仮説検定」(statistical hypothesis testing)の数学理論へと発展します。ただ,ネイマンたちの問題の設定のしかたを科学の研究にそのまま適用しようとすると,実験前に有意水準を例えば 0.05 と決めたなら,必ず p = 0.049 で帰無仮説を棄却して p = 0.051 で棄却しないという杓子定規な判断をすることになり,科学の研究者には違和感があります。
科学の歴史の中には,p 値が非常に小さかったことが後でひっくり返ることがいくらでもあります。p 値を求めるのは必要ですが,p ≦ 0.05 かどうかで機械的に判断せず,小さい p 値が一貫して再現できることをもって判断するのが科学の態度です。一方で,工場の品質管理の現場では,p 値によって自動判断することが必要とされます。
ニュートリノという素粒子の基礎研究で小柴昌俊氏が2002年にノーベル物理学賞を受賞しました。
そのニュートリノは今では質量を持つということで決着しています。過去の新聞報道を調べると,ニュートリノが質量を持つ確率は99%,99.99%と,次第に増えています。正しくは,ニュートリノが質量を持たないという帰無仮説から計算すると,実際の実験結果以上にニュートリノ質量の存在を支持するデータが偶然に生じるのは0.01%以下であるということになります(ここで「何々以上に」というときは「何々」を含みます)。

右の図は,分散がいずれも 1 で平均値が 0.2 だけ離れている二つの正規分布 N(-0.1, 1),N(0.1, 1) です。これらの群の平均値の違いを検出したければ,各群200個ずつ調べる必要があります。それだけ集めて調べれば統計的に有意な差が得られる公算が高いのですが,逆に,それだけ集めないと違いがわからないほどであれば実質的な差はほとんどないとも言えます。実際,各群からランダムに1個ずつ取り出したとき,左の群のものが右の群のものより大きくなる確率は44%以上もあります。
「実質的な差」がなくても,科学的に意味がないわけではありません。ノイズに埋もれて,非常に大きなサンプルを調べなければ統計的に有意な結果が得られないような現象が,科学的には非常に重要であることは,よくあることです。「どんなに p 値が小さくても,例えば相関係数が 0.1 といった値であれば,実質的な意味はない」はあくまで一般論で,ケースバイケースで考えるべきことです。
例えば,ニュートリノの質量が 0 でないという発見は,その質量がどんなに小さくても,ノーベル賞級の意味があります。また,血液型と性格に関連があることが示されれば,その関連がどんなに小さくても,性格の生理学的起源に関する重要な発見であると思われます(未だ一貫した結果は得られていません)。
なお,平均値の小さな違いでも,極端な値の割合に大きな違いが出ることがあります。さきほどの二つの正規分布で x ≧ 3 の確率を求めると,
> pnorm(-3.1) [1] 0.0009676032 > pnorm(-2.9) [1] 0.001865813
のように,千個に1個の割合が,千個に2個の割合に増えます。この違いが実質的な意味を持つことは十分考えられます(これらの数値は分布が正確な正規分布であることに強く依存するので,単なる例としての意味しかありません)。
右上の図は次のようにして描きました:
curve(dnorm(x-0.1),xlim=c(-4,4),frame.plot=FALSE,xlab="",ylab="") curve(dnorm(x+0.1),add=T)
各群からランダムに1個ずつ取り出したとき,左の群のものが右の群のものより大きくなる確率が44%以上もあることは,次のような数値積分で調べられます:
integrate(function(x){dnorm(x+0.1)*pnorm(x-0.1)},-Inf,+Inf)
あるいはシミュレーションで次のようにして調べられます:
mean(rnorm(1000000,mean=0.1) < rnorm(1000000,mean=-0.1))
7個の変数についてデータを収集しました。それらの二つずつについて相関係数を調べ, 7C2 = 21 個の相関係数を得ました。うち一つが統計的に有意でした。これで論文を書いていいでしょうか? ―― 駄目です。21個の中で偶然に5%水準で有意になる相関係数は,期待値として1個以上あります。こういうのを多重比較(multiple comparisons)といって,注意が必要です。
ポーカーをやったらツーペアが出ました。ツーペアが出る確率は約 1/21 なので,5%水準で有意です(何が???)。
血液型と性格の関係を調べるために,A/B/AB/O,男/女,いろいろな性格検査の項目をごちゃごちゃやっているうちに,いくつか統計的に有意な結果が見つかりました。 ―― こういうのを data dredging(データの浚渫(しゅんせつ))といいます。いろいろやっていれば,20回に1回は有意な結果が出るものです。
5%水準で有意な結果を載せる論文誌の内容は,20個に1個はウソです。本当? ―― ウソです。みんなが data dredging をしているなら,大部分がウソかもしれません。
じゃあ,統計はウソをつくための道具? ―― そうではありません。有意性検定の関門を設けることにより,ウソを減らそうとしているのです。でもウソを生産しようとする力も強いのが問題です。Data dredging の危険性を十分認識し,有意な結果でも鵜呑みにせず,追試により一貫して同じ結果が出ることを確かめることが大切です。
検定の概念が理解できれば,次に,モデルのパラメータ θ をいろいろ変えてみて,データがどの範囲のモデルと両立し得るかを調べてみましょう。

2項分布
で,以下では n = 10 の場合を考えます。このとき,例えば θ = 0.5 という帰無仮説とデータが水準 0.05 で両立し得る(帰無仮説が棄却されない)のは,r が 2 以上 8 以下のとき,つまり,割合 r/n が 0.2 以上 0.8 以下のときです。θ をいろいろ変えると,この範囲は右図のように階段状に変わります。色を付けた部分が,横軸 θ に対応する帰無仮説が水準 0.05 で棄却されない縦軸 r/n の範囲です。
右図は次のようにして描きました(色はPhotoshopで付けました):
par(las=1) curve(qbinom(0.025,10,x)/10,xlim=c(0,1),n=1001,asp=1) curve(qbinom(0.975,10,x)/10,n=1001,add=T)
ここで,逆にデータ r/n(縦軸)が与えられたときに図を横に見ていけば,そのデータと両立する θ の範囲がわかります。この範囲のことを,θ の「95%信頼区間」(confidence interval,CI)と呼びます。つまり,
と定めます。
「つまり θ の真の値は95%の確率でこの区間に入る」―― と言いたいところですが,そもそもモデル自体が正しいかどうかわかりませんし(というかモデルは常に近似に過ぎませんし),θ の「真の値」が存在したとしても,それは定数であって,確率変数ではありませんので,「θ が何々の確率でこの区間に入る」という意味がはっきりしないと思います。次のように考えるとわかりやすいでしょう。
まず,何でもいいから θ の値を一つ固定します。例えば θ = 0.5 としましょう。すると,2項分布モデルから計算して,95%以上(正確には後で述べるように約97.85%)の確率で r/n は 0.2 以上 0.8 以下になりますが,上の図をよく見て考えれば,r/n = 0.2 から 0.8 までに対応する信頼区間のどれにも θ = 0.5 が含まれます(当然!)ので,確かに95%以上の確率で,r/n(これが確率変数です)に対応する区間(これも確率変数になります)に,θ = 0.5(これは定数です)が含まれることになります。このことは θ をどんな値に固定しても言えることです。
さらに言い換えれば,モデル(ここでは2項分布)さえ正しければ,何度も何度も実験を繰り返して,毎回得られた r/n に対応する95%信頼区間を列挙し続ければ,そのうち95%(以上)が正しいということができます。こちらのほうを信頼区間の定義とすることもできます(その場合には「以上」を除くこともできますが,離散分布の場合に「以上」を除くのは簡単ではありません)。
上に示した97.85%という値は次のようなさまざまな方法で求められます。
> pbinom(8,10,0.5) - pbinom(1,10,0.5) [1] 0.9785156 > 1 - 2 * pbinom(1,10,0.5) [1] 0.9785156 > sum(dbinom(2:8,10,0.5)) [1] 0.9785156
以上が(2項分布のような離散的な場合も含めた)古典的な信頼区間の考え方です。より詳しくは次の文献をご覧下さい。
有意性の検定に基づく信頼区間の定義はシンプルですが,2項分布のような離散分布の場合,同じことをずっと続けた場合に θ が区間内に入る確率は95%以上になり,ぴったり95%になりません。原理的にどんな場合にもぴったり95%にすることは離散分布である限り無理なのですが,できるだけ95%に近づけようという研究がいろいろあります。詳しくは,例えば次の論文とそれに続くコメント集,特に Casella のコメント,そこで言及されている論文をご覧ください。
データから信頼区間を決めるやりかたそのものにランダムさを導入すれば,信頼区間に入る確率をぴったり95%にすることが可能ですが,データから一意的に定まらない信頼区間は奇異ですし,そこまでする必要はないでしょう。
Rには上の考え方に従って2項分布の検定と信頼区間を求める関数 binom.test()
があります。例えば実際に得られた割合が10個のうち2個だったとすると,
> binom.test(2,10,0.5)
Exact binomial test
data: 2 and 10
number of successes = 2, number of trials = 10, p-value = 0.1094
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.02521073 0.55609546
sample estimates:
probability of success
0.2
ということで,θ = 0.5 という帰無仮説についての p 値は 0.1094 で,θ の95%信頼区間は 0.025 から 0.556 であることがわかります(この範囲に θ = 0.5 が含まれることからも θ = 0.5 という帰無仮説が棄却されないことがわかります)。
もっと大きな場合にやってみましょう。
新聞社などの行う世論調査は,今では機械でランダムに電話をかけるRDDという方式がよく使われます。回答数は1000〜2000人程度です。
例えば1000人のうち200人が現内閣を支持すると答えたとすると,
> binom.test(200,1000) # 第3引数を省略すると0.5とされてしまう
Exact binomial test
data: 200 and 1000
number of successes = 200, number of trials = 1000, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1756206 0.2261594
sample estimates:
probability of success
0.2
というわけで,θ = 0.5 という帰無仮説についての p 値が 2.2×10-16 より小さいということはあまり役に立ちませんが,95%信頼区間が 0.176 から 0.226 であることがわかります。つまり,「20%」と言っても,実際は18%かも22%かもしれないので,小数点以下まで表示する意味はありません。
標本サイズ n が大きくなると,ギザギザな図は滑らかになり,2項分布は正規分布で近似できるようになるので,2項分布の分散が nθ(1 - θ) であることを使えば,観測された割合を r/n とすると,θ の95%信頼区間はほぼ r/n ± 1.96×sqrt((r/n)(1 - r/n) / n) となります。1000人のうち200人がYesと答えた場合には,この近似では θ = 0.20 ± 0.025 となり,さきほどの厳密な結果とほぼ一致します。
めでたしめでたし,という話にしたかったのですが,重箱の隅に,まだまだいろいろな魑魅魍魎がいます。
Rの binom.test() で p 値が 0.05
より小さくなっても,真の確率が95%信頼区間に含まれるということがありえます。どちらが正しいのでしょうか。
Rの binom.test() は信頼区間を上の Clopper and Pearson
の方法で求めています。この方法は,95%信頼区間なら,区間の下限は上側2.5%点が観測値になるような母数,区間の上限は下側2.5%点が観測値になるような母数をとるものです。これは両側 p
値を片側の2倍としていることに相当します。一方で,binom.test()
の出力する p 値は,下の Sterne の方法で求めています。これは,Fisher
の正確確率検定と同様に,与えられたモデルによる確率分布の,観測値での値を超えないすべての確率を合計する方式です。
さらに別の方式として,下の Blaker の方法は,もう一方の側の確率として,片側 p 値を超えない範囲で実際に得られる最大の値を使います。
一般に Sterne や Blaker の方法のほうが p 値が小さくなり,したがって信頼区間が狭くなるのですが,たとえば標本の大きさを増やして,p 値が減るはずの結果が出ても,実際には p 値が増えてしまうといった好ましくない現象も生じることがあるようです(下の Vos and Hudson 参照)。
Last modified: 2009-05-04 16:23:54