母比率の区間推定

母比率の推定

インターネット通販での商品評価で、「10人中8人が『良い』と答えた」とか「3人中3人が『良い』と答えた」というのを見かけます。単純に「良いと答えた人の比率」を見れば、前者は80%で後者は100%ですが、後者は3人しか答えていないので100%という数字もあまり信頼出来ません。
そこで、仮にその商品を買った人全員に丁寧にお願いして、良かったかどうか尋ねた場合の「良いと答えた人の比率」を、ランダムに選ばれた少数の回答から推測する問題を考えます。実際には、凄く良かったと感激した人や、凄く不満で腹がたった人が回答する可能性が高いので、「ランダムに選ばれた」も怪しいのですが、そこまで考えると「どれくらい感激、立腹するとどれくらい回答するか？」「どのようにすればランダムに選ばれた回答を得られるか」まで調べなければなりません。今回の講義の目的は比率の推定方法ですので、そこはばっさり省いてランダムに選ばれた回答だと見做せる問題について考えます。

区間推定

例えば「10人中8人が『良い』と答えた」場合、商品を買った人全員の中で満足している人の比率は恐らく80%と推定出来ますが、もしかしたら81%かもしれませんし79%かもしれません。
でも「いや、1000人中満足したのは8人だけで、今回は偶々その8人選ばれた」とか、「1000人中不満なのは2人だけで、今回は偶々その2人選ばれた」と言われると、その可能性はゼロではありませんが、そのようなことを言い出すと「満足している人の比率は約0%から約100%の間」という全く意味のない推定しか得られませんので、そのような確率の低いことは無視します。
どれくらい確率が低いことを無視するか、というのが信頼水準で、例えば信頼水準95%と言うと「最も起こりそうな95%だけに注目して、残り5%は無視する」ということになります。
「いやいや、5%も無視してはいけない」という場合は信頼水準を例えば99%にして、1%しか無視しないことにします。
言い換えると、ある程度確率の低いことを無視する時にしか信頼区間(に限らず検定なども)は使えなくて、「どんなに可能性の低いことも考慮しなければならない」場合には、最初に戻って「いや、1000人中満足したのは8人だけで、今回は偶々その8人選ばれた」場合も考慮して「満足している人の比率は約0%から約100%の間」という結論になります。

信頼区間の求め方

説明のために、10人に質問したら8人が満足と答えた場合、そして信頼水準を95%とします。他の値を使うときは以下の数字を読み替えてください。
「商品を買った人全員の中で満足している人の比率」として様々な値を想定します。そして各々の比率の値に対して、95%の確率で最も起こりやすそうな範囲を求めます。
例えば比率として80%を想定した場合は

人数	確率
0	0%
1	0%
2	0%
3	0%
4	1%
5	3%
6	9%
7	20%
8	30%
9	27%
10	11%

になります。四捨五入のために足して100%にならないのは気にしないでください。実際に計算するときにはもっと有効桁を増やしますから。
この中で、確率の高い方から順に、合計が95%を超えるまで足していきます。この場合は6人の場合から10人の場合まで足すと97%になって95%を超えます。従って

• もし10人に聞いて6人(あるいはそれ以上の人)が満足と答えたら、「満足している人の比率はもしかしたら80%で、今回は偶々6人が満足と答えた」と考えて、80%は信頼区間の中に入ります。
  
• もし10人に聞いて5人(あるいはそれ以下の人)が満足と答えたら、「満足している人の比率はもしかしたら80%なのに今回は運悪く5人しか満足と答えなかっただけかもしれない。でもその確率は低いので無視する」と考えて、80%は信頼区間の中に入りません。

例として想定している調査結果は10人中8人ですので、80%は信頼区間の中に入ります。
同様に様々な比率を想定して、それが「10人中8人」という結果に対する信頼区間に入るかどうか調べることで、信頼区間がどこからどこまでなのかが分かります。

計算法

実際にはこんな計算はやっていられないので、正規分布で近似して解析的に扱えるようにしてから逆算しています。How Not To Sort By Average Rating
でも回答者数が少ないときは正規近似が良くないですし、最近のパソコンは性能が上がっているので、近似に頼ることなく信頼区間の定義通りに計算するプログラムをExcel VBAを使って書いてみましょう。

まずエクセルのA1セルから以下のように入力します。

試行回数	10
発生回数	8
信頼水準	0.95
最低桁数	2
推定値
信頼下限
信頼上限

そして以下のプログラムを貼りつけて、ConfidenceIntervalを実行すると、推定値、信頼限界がB5～B7セルに書かれます。
受講生以外で、VBAを使ったことがない人はこの講義の第一回を見てください。
試行回数が多く、発生回数の比率が50%付近の場合は、正規近似とほぼ同じ結果になります。でも10回しか試行しなかったとか、1万回中1回しか発生しなかった、という場合は正規近似が良くないので、近似を使わずに直接求めるこのプログラムを使ってください。

Sub ConfidenceInterval()
' (c) 2009 Kaoru Fueda
Dim n As Long, x As Long, level As Double, decimals As Integer, mean As Double, sup As Double, inf As Double
n = Cells(1, 2) ' 標本数
If n < 1 Then MsgBox ("少なくとも1回は試行してください")
x = Cells(2, 2) '成功数
If x < 0 Or x > n Then MsgBox ("発生回数は0回以上かつ試行回数以下でしょう")
level = Cells(3, 2) 'level 信頼水準
If level < 0 Or level > 1 Then MsgBox ("信頼水準は確率ですから0から1の間です")
decimals = Cells(4, 2)
If decimals < 1 Then decimals = 1

decimals = decimals - 1
Do
　　decimals = decimals + 1
　　mean = x / n
　　sup = Application.Min(1, Round(mean + Application.NormInv(1 - (1 - level) / 2, 0, 1) * Sqr(mean * (1 - mean) / n), decimals + 1))
　　inf = Application.Max(0, Round(mean - Application.NormInv(1 - (1 - level) / 2, 0, 1) * Sqr(mean * (1 - mean) / n), decimals + 1))

　　Cells(5, 2) = Round(mean, decimals)

　　If shinraikukan(n, inf, level, 1) >= x Then
　　　　Do
　　　　　　inf = Application.Max(0, inf - 1 / (10 ^ (decimals + 1)))
　　　　Loop While inf > 0 And shinraikukan(n, inf, level, 1) >= x
　　　　If shinraikukan(n, inf, level, 1) < x Then
　　　　　　inf = inf + 1 / (10 ^ (decimals + 1))
　　　　End If
　　Else
　　　　Do
　　　　　　inf = inf + 1 / (10 ^ (decimals + 1))
　　　　Loop Until shinraikukan(n, inf, level, 1) >= x
　　End If
　　Cells(6, 2) = Round(inf, decimals)

　　If shinraikukan(n, sup, level, 0) <= x Then
　　　　Do
　　　　　　sup = Application.Min(1, sup + 1 / (10 ^ (decimals + 1)))
　　　　Loop While sup < 1 And shinraikukan(n, sup, level, 0) <= x
　　　　If shinraikukan(n, sup, level, 0) > x Then
　　　　　　sup = sup - 1 / (10 ^ (decimals + 1))
　　　　End If
　　Else
　　　　Do
　　　　　　sup = sup - 1 / (10 ^ (decimals + 1))
　　　　Loop Until shinraikukan(n, sup, level, 0) <= x
　　End If
　　Cells(7, 2) = Round(sup, decimals)
Loop Until ((x = n Or x = 0) And Cells(7, 2) > Cells(6, 2)) Or (Cells(5, 2) > Cells(6, 2) And Cells(7, 2) > Cells(5, 2))
End Sub

Function shinraikukan(n As Long, p As Double, level As Double, flag As Integer) As Long
' (c) 2009 Kaoru Fueda
'二項分布B(n,p)で確率level以上の最も短い区間の上端、下端を返す。
Dim upper As Long, lower As Long
'n 標本数
If n <= 0 Then n = 1
'p 母比率
If p < 0 Then p = 0
If p > 1 Then p = 1
'level 信頼水準
If level < 0 Then level = 0
If level > 1 Then level = 1
'一応変な値は訂正していますが、呼び出す側でチェックしてください。

upper = Round(n * p, 0)
lower = Round(n * p, 0)

'lower=0かつupper=nになったら、計算誤差によりlevelに達しなくても終了
Do While BinomialDistribution(upper, n, p, 1) - BinomialDistribution(lower - 1, n, p, 1) < level And (lower > 0 Or upper < n)
　　If lower = 0 Then
　　　　upper = upper + 1
　　Else
　　　　If upper = n Then
　　　　　　lower = lower - 1
　　　　Else
　　　　　　If BinomialDistribution(lower - 1, n, p, 0) > BinomialDistribution(upper + 1, n, p, 0) Then
　　　　　　　　lower = lower - 1
　　　　　　Else
　　　　　　　　If BinomialDistribution(lower - 1, n, p, 0) < BinomialDistribution(upper + 1, n, p, 0) Then
　　　　　　　　　　upper = upper + 1
　　　　　　　　Else
　　　　　　　　　　lower = lower - 1
　　　　　　　　　　upper = upper + 1
　　　　　　　　End If
　　　　　　End If
　　　　End If
　　End If
Loop
If flag = 0 Then
　　shinraikukan = lower
Else
　　shinraikukan = upper
End If
End Function

Function BinomialDistribution(x As Long, n As Long, p As Double, flag As Integer) As Double
' (c) 2009 Kaoru Fueda
'組み込み関数BinomDistは0≦x≦nを満たさなければエラーになるので、その点を定義通りに拡張
'BinomDistで計算出来無いほど数が多い場合は正規近似の近似精度は高いので正規近似を用いる
If flag = 0 Then
　　If x < 0 Or x > n Then
　　　　BinomialDistribution = 0
　　Else
　　　　temp = Application.BinomDist(x, n, p, 0)
　　　　If IsNumeric(temp) Then
　　　　　　BinomialDistribution = temp
　　　　Else
　　　　　　temp = Application.BinomDist(n - x, n, 1 - p, 0)
　　　　　　If IsNumeric(temp) Then
　　　　　　　　BinomialDistribution = temp
　　　　　　Else
　　　　　　　　BinomialDistribution = Application.NormDist(x, n * p, Sqr(n * p * (1 - p)), 0)
　　　　　　End If
　　　　End If
　　End If
Else
　　If x < 0 Then
　　　　BinomialDistribution = 0
　　Else
　　　　If x >= n Then
　　　　　　BinomialDistribution = 1
　　　　Else
　　　　　　temp = Application.BinomDist(x, n, p, 1)
　　　　　　If IsNumeric(temp) Then
　　　　　　　　BinomialDistribution = temp
　　　　　　Else
　　　　　　　　temp = Application.BinomDist(n - x - 1, n, 1 - p, 1)
　　　　　　　　If IsNumeric(temp) Then
　　　　　　　　　　BinomialDistribution = 1 - temp
　　　　　　　　Else
　　　　　　　　　　BinomialDistribution = Application.NormDist(x + 0.5, n * p, Sqr(n * p * (1 - p)), 1)
　　　　　　　　End If
　　　　　　End If
　　　　End If
　　End If
End If
End Function

計算例

10人中8人満足なら、満足している人の比率の95%信頼区間は45%～96%
100人中80人なら71%～87%
1000人中800人なら77%～82%になります。
80%に関して左右対称ではありませんね。比率50%付近はバラツキ(標準偏差)が大きく、0%,100%付近では小さいので、信頼区間は50%側が広く0%, 100%側が狭くなります。

練習

これまでswanを1万回観測したら1万回すべて白いswanだった。swanの中の「黒いswan」の比率を信頼水準99%で推定しなさい。
注意:点推定すれば比率は0と推定されます。