母比率の信頼区間
二項分布$B(n, p)$に従う確率変数$X$の母比率$p$の$100(1-\alpha)%$信頼区間は、標本比率を$\hat{p} = X/n$、標準正規分布の上側$100\alpha%$点を$z_\alpha$とすると、
\Biggl[\hat{p}-z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}, \hspace{8pt}\hat{p}+z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\Biggr]
となります。これは教科書でよく見る信頼区間の作り方で、wald信頼区間と呼ばれています。
この信頼区間はベルヌーイ分布$B(1, p)$従う独立な確率変数$X_1, X_2, \ldots, X_n$の標本平均$\bar{X}$が中心極限定理により正規分布$N(p, p(1-p)/n)$に分布収束することから導かれます。したがって、上記の公式はサンプルサイズが十分に大きい時に成り立つ式であることがわかります。しかし、現実にはあまりサンプルサイズが大きくない場合にも母比率の信頼区間を用いて何かしらの評価したいことがあると思います。そこで、そこでサンプルサイズをいろいろ変えて信頼区間の性能を評価してみます。
今回の記事では信頼区間の性能を評価するのに被覆確率 (coverage probability)を使用します。被覆確率とは信頼区間が真の母比率$p$を含む確率のことです。そこで、二項分布$B(n, p)$からたくさんの乱数を発生させ、たくさんの信頼区間を作ります。その中で真の母比率$p$を含んでいる信頼区間の数を数えて、その割合を被覆確率の推定値とします。サンプルサイズ$N$と母比率$p$をいろいろ変えながらそれぞれ$10000$個の$95$%wald信頼区間を作り、被覆確率を計算してみます。
上のグラフで赤い線は被覆確率が$0.95$の位置を示しています。$95$%信頼区間を作っているので被覆確率は$0.95$になって欲しいです。しかし、上のグラフを見ると特にサンプルサイズが小さい時は被覆確率が$0.95$を大きく下回る場合があることがわかります。つまり、サンプルサイズが小さい時にwald信頼区間を用いると「$95$%信頼区間を作りたいのに実は$60$%信頼区間を作ってしまっている」みたいな状況になってしまうことがあり得ます。また、上のグラフを見るとサンプルサイズがある程度大きくても母比率が$0$や$1$に近いと被覆確率が小さくなってしまうことがわかります。wald信頼区間は教科書にもよく出てくるし、多くのソフトウェアで標準的に出力されたりしますが、（サンプルサイズが小さい時は特に）それほど性質の良いものではないことがわかります。
実は二項分布のような離散的な確率分布では真のパラメータの値によって信頼区間の被覆確率が変化してしまうことが知られています。真のパラメータの値が分からないから信頼区間で推定しようとしているのにこれでは困ります。真のパラメータがどんな値であっても被覆確率が目的の信頼水準を下回らないような区間が作れないなぁ、という問題の解決策（の一つ）がClopper-Pearsonの信頼区間です。
Clopper-Pearsonの信頼区間
Clopper-Pearson信頼区間の構成方法
Clopper-Pearsonの信頼区間は二項分布に基づく正確な検定によって構成される信頼区間です。正確な検定とは正規分布への近似を利用しないという意味です。次のような具体例を用いてClopper-Pearsonの信頼区間について説明します。成功確率$\theta$のベルヌーイ試行を$5$回行ったとすると成功回数$X$は二項分布$B(5, \theta)$に従います。仮に成功回数の実現値が$2$回だったとして$\theta$の95%Clopper-Pearson信頼区間を構成してみることにします。$95$%信頼区間なので、$\alpha=0.05$です。まずは信頼区間の上限を調べるために以下のような片側検定を考えます。有意水準は$\alpha/2 =0.05/2= 0.025$とします。
H_0: \theta = \theta_0, \hspace{10pt} H_1: \theta < \theta_0
ここで、$\theta_0$がいくつの時に帰無仮説が棄却されるか調べてみます。例えば、$\theta_0=0.2$とした時の$p$値は、試行回数5の二項分布の確率質量関数を$P(x; \theta_0)$とすると、
P(0; \hspace{2pt}0.2) + P(1; \hspace{2pt}0.2) + P(2;\hspace{2pt} 0.2) = 0.942
となるので帰無仮説は棄却されません。一方で$\theta_0=0.9$とした時の$p$値は
P(0;\hspace{2pt} 0.9) + P(1;\hspace{2pt} 0.9) + P(2;\hspace{2pt} 0.9) = 0.008
となるので帰無仮説は棄却されます。では、$\theta_0$の値をいろいろ変えた時の$p$値の変化をプロットして確認してみることにします。
上のグラフを見るとおよそ$\theta_0=0.85$より大きくなると$p$値が$0.025$を下回り、帰無仮説が棄却されることがわかります。このように考えると「帰無仮説が棄却されない$\theta_0$の上限」$\fallingdotseq$「得られたデータと矛盾しない$\theta_0$の上限」が分かります。このようにして算出した$0.85$を$95$%信頼区間の上限とします。
続いて信頼区間の下限を算出する時は以下のような先ほどとは逆向きの片側検定を考えます。有意水準は同様に$0.05/2 = 0.025$とします。
H_0: \theta = \theta_0, \hspace{10pt} H_1: \theta > \theta_0
すると、$p$値は以下のように計算することができます。
P(2; \theta_0) + P(3;\theta_0) + P(4;\theta_0) + P(5;\theta_0)
では、先ほどと同様に$\theta_0$の値をいろいろ変えた時の$p$値の変化をプロットして確認してみることにします。
上のグラフを見ると大体$\theta_0=0.05$以下になると$p$値が$0.025$を下回り、帰無仮説が棄却されることがわかります。したがって、データと矛盾しない$\theta_0$の下限である$0.05$を$95$%信頼区間の下限とします。このように算出した上限と下限を用いて「試行回数$5$回、成功回数$2$回」というデータが得られた時の母比率$\theta$の$95$%信頼区間を$[0.05, 0.85]$とするのがClopper-Pearsonの信頼区間です。
ここまでの信頼区間の導出方法を一般化してClopper-Pearsonの信頼区間を定義しておきましょう。ベルヌーイ試行$B(1,\theta)$を$n$回行い、成功回数が$x$回であったとします。このとき母比率$\theta$の$100(1-\alpha)%$Clopper-Pearson信頼区間は、
\sum_{i=x}^{n}{}_nC_i\theta^i(1-\theta)^{n-i} = \frac{\alpha}{2}
の解$\theta=\theta_L$を信頼区間の下限、
\sum_{i=0}^{x}{}_nC_i\theta^i(1-\theta)^{n-i} = \frac{\alpha}{2}
の解$\theta=\theta_U$を信頼区間の上限として、$[\theta_L, \theta_U]$となります。
このように信頼区間を構成することで母比率がどんな値であっても被覆確率が目的の信頼水準を下回ることを避けることができます。Clopper-Pearsonの信頼区間は二項分布に基づく正確な検定によって構成されるので「Clopper-Pearsonの正確信頼区間」と呼ばれることがあります。「正確」とは正規分布への近似をしていないという意味なので、得られる信頼区間の正確さとは関係ないようです。実際に$95$%信頼区間を作っても被覆確率がピッタリ$0.95$になるわけではありません。
被覆確率
サンプルサイズ$N$と母比率$p$をいろいろ変えながら$95$%Clopper-Pearson信頼区間の被覆確率を見てみましょう。
上のグラフで赤い線は被覆確率が$0.95$の位置を示しています。Clopper-Pearsonの信頼区間はサンプルサイズが小さい時でも被覆確率が$0.95$を下回っていないことがわかります。サンプルサイズが大きくなると、ときどき$0.95$を下回っていることがありますがシミュレーションで生成した信頼区間の数が$10000$個と少ないことが原因だと思われます。一方で、サンプルサイズが小さい時はClopper-Pearson信頼区間の被覆確率は$0.95$よりも大きめになり、信頼区間が必要以上に広くなってしまうことが知られています。
ベータ分布による表現
Clopper-Pearsonの信頼区間はベータ分布を用いて表現することが可能です。二項分布$B(n, \theta)$とベータ分布には以下のような関係が成立することが知られています。
\sum_{i=x}^{n}{}_nC_i\theta^i(1-\theta)^{n-i} = \int_{0}^{\theta}f(t;x, n-x+1)dt
ここで、$f(t;\alpha,\beta)$はベータ分布$Beta(\alpha,\beta)$の確率密度関数です。したがって、右辺はベータ分布の累積分布関数になっています。Clopper-Pearson信頼区間の定義とこの関係を利用するとClopper-Pearson信頼区間の下限は
\sum_{i=x}^{n}{}_nC_i\theta^i(1-\theta)^{n-i} = \int_{0}^{\theta}f(t;x, n-x+1)dt= \frac{\alpha}{2}
を満たす$\theta = \theta_L$となります。また、信頼区間の上限は
\sum_{i=0}^{x}{}_nC_i\theta^i(1-\theta)^{n-i} = 1-\sum_{i=x+1}^{n}{}_nC_i\theta^i(1-\theta)^{n-i}= \frac{\alpha}{2}
\int_{0}^{\theta}f(t;x+1, n-x)dt=1-\frac{\alpha}{2}
を満たす$\theta = \theta_U$となります。そこで、ベータ分布$Beta(\alpha,\beta)$の$100q$パーセント点を$Beta(q;\alpha,\beta)$とするとClopper-Pearson信頼区間は
\biggl[Beta\Bigl(\frac{\alpha}{2}; x, n-x+1\Bigr), \hspace{5pt}Beta\Bigl(1-\frac{\alpha}{2}; x+1, n-x\Bigr)\biggr]
のように表現することができます。そのため、Clopper-Pearsonの信頼区間を「ベータ分布に基づく信頼区間」と表現することがあります。この式を用いると、$5$回のベルヌーイ試行の成功回数が$2$回だった時の成功確率の$95$%信頼区間を以下のように計算できます。
from scipy.stats import beta
alpha = 0.05 # 95%信頼区間
n = 5 # 試行回数
x = 2 # 成功回数
lower_limit = beta.ppf(alpha/2, x, n-x+1)
upper_limit = beta.ppf(1-alpha/2, x+1, n-x)
print(f'95% Confidence interval: [{lower_limit:.5f}, {upper_limit:.5f}]')
> 95% Confidence interval: [0.05274, 0.85337]
F分布による表現
Clopper-Pearsonの信頼区間はF分布を用いて表現することも可能です。二項分布$B(n, \theta)$とF分布には$x = 1,\ldots,n-1$に対して、
\sum_{i=0}^{x}{}_nC_i\theta^i(1-\theta)^{n-i} = 1-\int_{0}^{m}f(t;\phi_1, \phi_2)dt
という関係があることが知られています。ここで、$f(t;\phi_1, \phi_2)$はF分布$F(\phi_1, \phi_2)$の確率密度関数、$\phi_1=2(x+1)$、$\phi_2=2(n-x)$、$m=\frac{\phi_2\theta}{\phi_1(1-\theta)}$となっています。右辺第2項はF分布の累積分布関数になっていることがわかります。そこで、Clopper-Pearson信頼区間の上限は定義とこの関係をそのまま用いて、
\sum_{i=0}^{x}{}_nC_i\theta^i(1-\theta)^{n-i} = 1-\int_{0}^{m}f(t;\phi_1, \phi_2)dt = \frac{\alpha}{2}
\int_{0}^{m}f(t;\phi_1, \phi_2)dt = 1-\frac{\alpha}{2}
となるようなF分布の$100(1-\alpha/2)$パーセント点$m$を求めて、これを$\theta$について解くことで求めることができます。また、Clopper-Pearson信頼区間の下限は
\sum_{i=x}^{n}{}_nC_i\theta^i(1-\theta)^{n-i} =1-\sum_{i=0}^{x-1}{}_nC_i\theta^i(1-\theta)^{n-i}
=1-\Bigl\{1-\int_{0}^{m}f(t;\phi_1, \phi_2)dt\Bigr\}=\int_{0}^{m}f(t;\phi_1, \phi_2)dt=\frac{\alpha}{2}
の関係を用いてF分布の$100(\alpha/2)$パーセント点$m$を$\theta$について解くことで求めることができます。したがって、Clopper-Pearsonの信頼区間はF分布$F(\phi_1, \phi_2)$の$100q$パーセント点を$F(q; \phi_1, \phi_2)$とすると、
\Biggl[\Biggl\{1+\frac{n-x+1}{xF\bigl[\frac{\alpha}{2};2x,2(n-x+1)\bigr]}\Biggr\}^{-1},\hspace{5pt}
\Biggl\{1+\frac{n-x}{(x+1)F\bigl[1-\frac{\alpha}{2};2(x+1),2(n-x)\bigr]}\Biggr\}^{-1} \Biggr]
のようにF分布を用いて表現することができます。そのため、Clopper-Pearsonの信頼区間を「F分布に基づく信頼区間」と表現することがあります。この式を用いると、$5$回のベルヌーイ試行の成功回数が$2$回だった時の成功確率の$95$%信頼区間を以下のように計算できます。
from scipy.stats import f
alpha = 0.05 # 95%信頼区間
n = 5 # 試行回数
x = 2 # 成功回数
lower_limit = (1 + ((n-x+1)/(x*f.ppf(alpha/2, 2*x, 2*n-2*x+2))))**-1
upper_limit = (1 + ((n-x)/((x+1)*f.ppf(1-alpha/2, 2*x+2, 2*n-2*x))))**-1
print(f'95% Confidence interval: [{lower_limit:.5f}, {upper_limit:.5f}]')
> 95% Confidence interval: [0.05274, 0.85337]
終わりに
今回は最もよく用いられる（と思われる）Clopper-Pearsonの信頼区間について調べました。母比率の信頼区間の構成方法はかなりたくさんあって、Clopper-Pearsonよりも良い性質の方法もあったりします。
参考文献
- Måns Thulin. "The cost of using exact confidence intervals for a binomial proportion." Electron. J. Statist. 8 (1) 817 - 840, 2014. https://doi.org/10.1214/14-EJS909
- Pastor, Dominique. (2005). On the coverage probability of the Clopper-Pearson confidence interval.
- https://www.math.kyoto-u.ac.jp/~ichiro/lectures/2012st.pdf
記事内のグラフ作成に用いたプログラム：