最終更新日: 　counter: (2000/5/24からの累積)

母比率の信頼区間の求め方　(２項分布)

１　力による計算
２　Ｆ分布、ベータ分布関数使用による計算
　(1)Ｆ分布を使用する方法
　(2)ベータ分布を使用する方法
《引用文献》

母比率の信頼区間を求める方法としてＦ分布を使用するものがある。ところがその根拠をきちんと書いてあるものがなかなか見つからない。

問題をきちんと述べると、「Ｘ＝x が観測されたときの、２項分布の母数p に対する信頼区間１－αの上側信頼限界と下側信頼限界、または両側信頼限界を求める」方法である。竹内・藤野(1981:p158)ではこれらを「正確な信頼限界」または「Clopeer & Pearson(1984)の信頼限界」と呼ぶといっている。

１　力による計算

母比率の信頼区間については、Clopper & Peason(1934) にグラフが示されているので、ある意味では簡単な問題なのだろう。

いろいろチェックした結果、Blyth(1986:p843)には素朴に求める方法を記載していた。

n=5, α=.05 の上側信頼限界を求める
Ｘ＝０　p^(.95)(0) =.45072=1-(.05)^1/5
　　　　　　これは(1-p)⁵=.05 を解いたもの

Ｘ＝１　p^(.95)(1) =.65741
　　　　　　これは(1-p)⁵+5p(1-p)⁴=.05 を解いたもの

Ｘ＝２　p^(.95)(2) =.81075
　　　　　　これは10p³(1-p)²+5p⁴(1-p)+p⁵=.95 を解いたもの

Ｘ＝３　p^(.95)(3) =.92356
　　　　　　5p⁴(1-p)+p⁵=.95　を解いたもの

Ｘ＝４　p^(.95)(4) =.98979=(.95)^1/5
　　　　　　p⁵=.95 を解いたもの

Ｘ＝５　p^(.95)(5) =1

下側信頼限界は次のようにして上側信頼限界から求めることができる。
p_(1-α)(x)=1-p^(1-α)(n-x)　 for x= 0, 1,...., n.

２　Ｆ分布、ベータ分布関数使用による計算

もちろん、実際の計算ではこの式を正直に解くのではなく、２項分布とＦ分布の関係を使うことができる。Blyth(1986)ではＦ分布の表を使うのが大変であるとしているが、今はexcel のような表計算ソフトにもＦ分布関数をもっているので簡単である。実際は２項分布との関係はＦ分布よりもベータ分布のほうが単純に表すことができる。同じく excel ではベータ分布関数があり、簡単に適用できる。SASのような高価なパッケージを使う必要もないし、大きなプログラムをつくる必要もない。

母比率の信頼区間の求め方にはいくつか近似法がある(Blyth,1986)が、簡単にパソコンを使える環境においてわざわざ近似法を用いる必要もないだろう。

(1)Ｆ分布を使用する方法

Ｆ分布を使用する実際の計算はZar(1996:p136), Blyth(1986:p845-846)に丁寧な説明がある。内田(1996:p93)にはexcel の計算法を示している。

Ｆ分布との関連づけの導出は廣津(1982:p4-)にあるが、蓑谷（1998)で少し補充する必要があるだろう。

excel を使って信頼区間を求める
試行数,n, 100
生起数,x, 7
α,alpha,=0.05
下側信頼限界, ,
自由度１,dfl_1,= 2*(n-x+1)
自由度２,dfl_2,= 2*x
PL,PL,"=dfl_2/(dfl_1*finv(alpha/2,dfl_1,dfl_2)+dfl_2)"
上側信頼限界,,
自由度１,dfu_1,= 2*(x+1)
自由度２,dfu_2,= 2*(n-x)
PU, PU,"=dfu_1*finv(alpha/2,dfu_1,dfu_2)/(dfu_1*finv(alpha/2,dfu_1,dfu_2)+dfu_2)"
これをベータ分布で行うと、,,
下側信頼限界,PBL,"=betainv(alpha/2,x,n-x+1)"
上側信頼限界,PBU,"=betainv(1-alpha/2,x+1,n-x)"

excel へ移植法: (1)範囲指定して、copy(ctrl+c); (2)excel に貼付(ctrl+v); (3)データ→区切り位置（区切り位置ウィザード起動）→次へ→区切り文字を「カンマ」に→次へ→次へ
　　（これでデータがセルにおさまり計算されている）; (4)２列目が名前になるように、範囲指定→挿入→名前→作成→（左端列にチェックマーク）→ok; 　範囲指定は実験できるように、右および下を広くとっておく。; (4)の操作なしでも計算結果はしているけれどその計算は違っている場合があるので、かならず(4)を実行すること！！

excel 計算結果

試行数	n	100
生起数	x	7
α	alpha	0.05
信頼下限
自由度１	dfl_1	188
自由度２	dfl_2	14
PL	PL	0.028605222
信頼上限
自由度１	dfu_1	16
自由度２	dfu_2	186
PU	PU	0.138919802
これをベータ分布で行うと、
下限	PBL	0.028605223
上限	PBU	0.13891983

上の表のようにF分布を使った計算(PL,PU)でもベータ分布を使った計算(PBL,PBU)でも答えは(小数点以下7桁の範囲内=計算誤差の範囲内で)一致している。

ベータ分布のルーティンのほうが簡略なので誤差の累積は少ないはずだし、ベータ分布を使用した方がはるかに計算の指定が簡単であることからも、ベータ分布を使用するほうがいい。

実際の適用では注意が必要である。なお、次の式になることの証明はBlyth(1986:p845)にある。

上側信頼限界は excel の関数で示せば、=invbeta(1-alpha/2,k+1,n-k) となる。

上側信頼限界PBUは、計算上は下側確率alpha/2 を求めるので、こうなる。この文学的表現でいいのかな？

下側信頼限界PBL は、計算上は上側確率を求める。しかも、これが下側信頼限界であるから確率は(1-alpha/2)ではなく、alpha/2 になる。パラメータk,n-k+1のベータ分布から求められたp が下側信頼限界となる。エクセルの関数では上にあるように、=betainv(alpha/2,k,n-k+1)で p を求めることができる。文学的表現は危ないから無視して上のBlyth の引用している証明を読んでください。

下側信頼限界PBL は、X=kのときにこの点kが上側確率Ｐ（例えば0.025)となるように、母比率pをベータ分布から求める。

(2)ベータ分布を使用する方法

２項分布とベータ分布との関係　（蓑谷(1998:p172)を参考にした）

nを試行回数、pを母比率（成功確率）, xを 0 以上 n 以下の任意の整数(成功回数）とする。

Ｘ～b(n,p)　（Ｘはパラメータ n, p の２項分布に従う確率変数）,
Ｙ～BETA(x,n-x+1) （Ｙはパラメータx,n-x+1のベータ分布に従う確率変数）とする。

このときＸの上側確率Ｐ(Ｘ＞＝x)を、次のようにしてベータ分布から求めることができる。

　　　Ｐ（Ｘ>=x)＝Ｐ(Ｙ＜＝p)

xのときの上側確率αの母比率をp^(α)(x)

xのときの下側確率αの母比率をp_(α)(x)
と表す。

p_(1-α)(x)=p^(α)(x-1)
　　　　　for x=1,2,...,n.

p_(1-α)(x)=1-p^(1-α)(n-x)
　　　　　for x=0,1,...,n.

この２つの式から次の式が導きだされる。

p^(1-α)(x)=1-p^(α)(n-x-1)
　　　　　for x=0,1,...,n-1.

２項分布の累積確率は母比率の分かっている場合なので、上の関係を素直に使うことができます。

ただし、計算はちょっと面倒になります。
0.025の下側信頼限界を求める場合、母比率 p が固定で、k を 0から１ずつ増加させて累積が0.025 になるところを求めます。

２項分布の信頼限界の値を出力する関数がexcel にはあるのでこれを使うと簡単に求めることができます。

=CRITBINOM(n,母比率,alpha/2) で下側信頼限界の個数k
=CRITBINOM(n,母比率,1-alpha/2) で上側信頼限界の個数k

を求めることができます。ただし、この求め方では母比率によって異なりますが、
n=1200 程度で限界になって、オーバーフローします。
青木繁伸氏の「二項分布の確率」を求めるプログラムもだいたい同じです。（まったく同じかをチェックしていません）
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/binomial.html

組み合わせの計算を工夫している岡本安晴さんの『delphi で学ぶデータ分析法』ＣＱ出版社　の２項分布のプログラムでは大幅に限界が緩和されています。 n=16,000程度あたりまで求めることができます。普通はこれで十分だと思います。

東京大学教養部統計学教室編『基礎統計学１　統計学入門』東京大学出版会 1991　p170-171

中心極限定理の応用の二項分布の正規分布の近似の節において

40000回コインを投げて,20,400回以上、あるいは19,600以下表がでることは、どの程度の確率であろうか。

という例を、コンピュータで計算することが不可能な問題としてあげています。40,000回コインを投げて見ろと言いたくなりますが、残念なことにコンピュータ実験ではこれくらいの回数は簡単です。では、コンピュータで確率を求めるのが困難かというとそうでもありません。ベータ分布を使って２項分布の確率を求めれば、それよりも１０倍多い、n=400,000 でもへっちゃらです。

excel で関数を作ってみました。ダウンロード

使い方はexcel の ex のついているcritbinomex,binomdistexはそれぞれ，critbinom,binomdistと同じです。
binominv は次のようになってます。
=binominv(試行回数, 成功率_母比率, 確率)
'返値：成功回数の比率（母比率より上と下で違ってくる。'
' 母比率より下では(1-確率)の下側確率を確保
' 母比率より上では(確率)の上側確率を確保

binominvを直接使う必要性はあまりないでしょう。

Function critbinomex(試行回数, 成功率_母比率, 確率)
'返値：限界値
'=CRITBINOM(100,0.5,0.05/2)
n = 試行回数
p0 = 成功率_母比率
p = 確率
'前処理
If (p0 <= 0) Or (p0 >= 1) Then
critbinomex = "#NUM!:母比率は０と１の間です"
Exit Function
End If
If (p < 0) Or (p > 1) Then
critbinomex = "#NUM!:上側累積確率は０以上１以下の値です"
Exit Function
End If
n = Int(n) 'nは整数。小数点以下は切り捨て
If n <= 0 Then
critbinomex = "#NUM!試行回数は１以上の数です"
Exit Function
End If
'本処理
tmp = binominv(n, p0, p)
critbinomex = Int(tmp * n) '母比率より高いときは小数点以下切り捨て
If (tmp < p0) And (tmp <> critbinomex) Then '母比率より低いときは小数点以下切り上げ
critbinomex = critbinomex + 1
End If
End Function

Function binominv(試行回数, 成功率_母比率, 確率)
'返値：成功回数の比率（母比率より上と下で違ってくる。'
' 母比率より下では(1-確率)の下側確率を確保
' 母比率より上では(確率)の上側確率を確保
'
'=betainv(alpha/2,k,n-k+1)
n = 試行回数
p0 = 成功率_母比率
p = 確率

esp = 0.000000000001: '1e-12
If (p0 >= 1#) Or (p0 <= 0#) Or (p >= 1#) Or (p <= 0#) Then
binominv = "#NUM!"
Exit Function
End If
xl = 0#
xr = n
k = 0

Do
k = k + 1
xm = (xl + xr) * 0.5
If p0 * n >= xm Then
temp = Application.BetaInv(1 - p, xm + 1, n - xm):
Else
temp = Application.BetaInv(1 - p, xm, n - xm + 1)
End If
If (temp < p0) Then
xl = xm
Else
xr = xm
End If
If xm = 0 Then Exit Do
If k >= 100 Then Exit Do
Loop Until ((Abs(xr - xl) / n) < esp)
binominv = xm / n
If binominv < esp Then binominv = 0
If (1 - binominv) < esp And xr = n Then binominv = 1
If k >= 100 Then
If xl = 0 Then
binominv = 0
Else: If xr = n Then binominv = 100 Else binominv = "収束しませんでした:k>=100:xl" + Str(xl) + ":x=" + Str(xr)
End If
End If
End Function

Function binomdistex(成功数, 試行回数, 成功率, 関数形式)
'BINOMDIST(成功数, 試行回数, 成功率, 関数形式)
'BETADIST(x, α, β, A, B)
x = 成功数
n = 試行回数
p0 = 成功率
ftype = 関数形式
If p0 > 1 Or p0 < 0 Or n < x Or n < 1 Then
binomdistex = "#NUM"
Exit Function
End If
If n < x Or n < 1 Then
binomdistex = "#NUM:試行回数は1以上かつ成功数以上でなけばなりません"
Exit Function
End If

binomdistex = 0
If x = n Then
binomdistex = 1
Else
binomdistex = 1 - Application.BetaDist(p0, x + 1, n - x)
End If

If Not (ftype) Then
If x = 0 Then Exit Function
y = 1 - Application.BetaDist(p0, x, n - x + 1)
binomdistex = binomdistex - y
End If
End Function

コラム
excel BINOMDIST(成功数, 試行回数, 成功率, 関数形式)関数の計算のおかしなところ

=binomdist(7,10,0,true)
=binomdist(7,10,1,true)

どちらも 0 を返してきます(trueは累積を求める）。母比率0.0 のときは７のときの累積は1.00にならないといけない。

母比率の信頼区間を求める時、そのkのときに両側でα=0.05 になるようなつまり片側2.5％の危険率になるような母比率を求めています。

《引用文献・参考文献》

Blyth,C.R.(1986) Approximate binomial confidence limits. Journal of the American Statistical Asscoiation, 81, 843-855.
Clopper,C.J., and Pearson,E.S.(1934) The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika, 26, 404-413.
Hahm,G.J., and Meeker,W.Q.(1991) Statistical intervals: A guide for practitioners. Wiley.
廣津千尋『離散データ解析』教育出版　1982
蓑谷千凰彦『すぐに役立つ統計分布』東京図書 1998
竹内啓・藤野和建『２項分布とポアソン分布』東京大学出版会　1981
田中豊・垂水共之『Windows版統計解析ハンドブック基礎統計』共立出版 1997
東京大学教養部統計学教室編『基礎統計学１　統計学入門』東京大学出版会 1991
内田治『すぐわかるＥＸＣＥＬによる統計解析』東京図書 1996
Zar,J.H. 1996 Biostatistical analysis. 3rd ed. Prentice-Hall

堀　啓造(home page) e-mail hori@ec.kagawa-u.ac.jp

母比率の信頼区間の求め方 (２項分布)