二項分布の信頼区間の計算方法の解説と精度比較
はじめに
今回は二項分布の信頼区間の計算方法について解説します。その計算方法は、主として厳密に信頼区間を計算する Clopper-Pearson 法と二項分布を正規分布に近似して計算する方法が２種類、計３通りの方法が存在します。今回はこれらの計算方法とPythonを用いたシミュレーションにより信頼区間の精度（真値が信頼区間にどの程度含まれるか）の比較結果について説明します。
概要
各手法の計算式
二項分布 $${B(n,p)}$$ に従う確率変数から得られた観測値を $${x_{0}}$$ とします。成功確率 $${p}$$ は未知とします。
このとき、$${p}$$ の$${100(1-\alpha)}$$ % 信頼区間の下限、上限を順に $${p_{LL},p_{UL}}$$ とすると、各手法における $${p_{LL},p_{UL}}$$ は以下の式を満たします。
1. Clopper-Pearson 法の信頼区間
二項分布の式を用いた場合
$$
\begin{aligned}
\sum_{k=x_{0}}^{n}{}_{n}C_{k} p_{LL}^{k}(1-p_{LL})^{n-k}
&=\alpha/2 \\
\sum_{k=0}^{x_{0}}{}_{n}C_{k} p_{UL}^{k}(1-p_{UL})^{n-k}
&=\alpha/2 \\
\end{aligned}
$$
ベータ分布の式を用いた場合
$$
\begin{aligned}
\int_{0}^{p_{LL}} \dfrac{1}{\Beta(x_{0},n-x_{0}+1)} p^{x_{0}-1}(1-p)^{n-x_{0}}dp =F_{\Beta(x_{0},n-x_{0}+1)}(p_{LL}) &=\alpha/2 \\
\int_{0}^{p_{UL}} \dfrac{1}{\Beta(x_{0}+1,n-x_{0})} p^{x_{0}}(1-p)^{n-x_{0}-1}dp =F_{\Beta(x_{0}+1,n-x_{0})}(p_{UL}) &=1-\alpha/2 \\
\end{aligned}
$$
$${F_{\Beta(a,b)}(p)}$$ は $${Beta(a,b)}$$ に従う分布の累積分布関数を表します。
詳細は "1. Clopper-Pearson 法 による信頼区間の計算" をご覧ください。
2. 正規近似による方法１
$$
\left\{\begin{aligned}
p_{LL}&=\dfrac{n(2x_{0}+z_{\alpha/2}^{2})-\sqrt{D}}{2n(n+z_{\alpha/2}^{2})} \\
p_{UL}&=\dfrac{n(2x_{0}+z_{\alpha/2}^{2})+\sqrt{D}}{2n(n+z_{\alpha/2}^{2})} \\
D&=nz_{\alpha/2}^{2}(nz_{\alpha/2}^{2}+4x_{0}(n-x_{0}))
\end{aligned} \right.
$$
$${z_{\alpha}}$$ は標準正規分布の上側 $${100\alpha}$$ % 点を表します。
詳細は "2. 正規近似による二項分布の信頼区間の計算方法１" をご覧ください。
3. 正規近似による方法２
$$
\left\{\begin{aligned}
p_{LL}&=\dfrac{x_{0}}{n}-\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
p_{UL}&=\dfrac{x_{0}}{n}+\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
\end{aligned} \right.
$$
詳細は "3. 正規近似による二項分布の信頼区間の計算方法２" をご覧ください
各方法の信頼区間の精度の比較結果
簡易シミュレーションにより、”二項分布に従う乱数から上記３方法による信頼区間の計算”を１０万回繰り返した結果、信頼区間に真値が含まれる確率が最も高いのはClopper-Pearson 法、真値が含まれる確率が最も信頼係数に近いのは正規近似による方法１という結果になりました。
詳細は "4. Python シミュレーションによる検証" をご覧ください
1. Clopper-Pearson 法 による信頼区間の計算
Clopper-Pearson 法は一般的な分布の信頼区間の計算方法に基づいたものです。よく信頼区間の計算例として取り上げられる”分散が既知の正規分布より得られた観測値から平均の信頼区間を求める” 場合の計算も同様の方法を用いています。そこで、Clopper-Pearsonの方法を紹介する前に、まず一般的な分布（未知パラメータが１つ）の場合の信頼区間の計算方法について概説し、次に上記の正規分布の平均の信頼区間の例がその計算方法に基づいていることを説明してから、Clopper-Pearson 法について説明します。
1.1. 未知パラメータが１つの確率変数の信頼区間の計算方法
一般に未知パラメータが $${\theta}$$ の１つのみで構成される確率変数 $${X}$$ について確率関数を $${f_{X}(x;\theta)}$$ 、累積分布関数を $${F_{X}(x;\theta)}$$ とし、確率変数 $${X}$$ から観測値 $${x_{0}}$$ が得られたとします。このとき、$${\theta}$$ の $${100(1-\alpha)}$$ % 信頼区間の下限 $${\theta_{LL}}$$ 、上限 $${\theta_{UL}}$$ については多くの場合、以下の式で表されます。
$${x_{0}}$$ における累積分布関数の値 $${F_{X}(x_{0};\theta)=P(X \leq x_{0}|\theta)}$$ が $${\theta}$$ について単調増加のとき、
$$
\left\{\begin{aligned}
P(X\leq x_{0}|\theta_{LL})&=\alpha/2\\
P(X\geq x_{0}|\theta_{UL})&=\alpha/2\\
\end{aligned} \hspace{25pt} (1.1a) \right.
$$
$${x_{0}}$$ における累積分布関数の値 $${F_{X}(x_{0};\theta)=P(X \leq x_{0}|\theta)}$$ が $${\theta}$$ について単調減少のとき、
$$
\left\{\begin{aligned}
P(X\geq x_{0}|\theta_{LL})&=\alpha/2\\
P(X\leq x_{0}|\theta_{UL})&=\alpha/2\\
\end{aligned} \hspace{25pt} (1.1b) \right.
$$
$${P(X\leq x_{0}|\theta)}$$ はパラメータが $${\theta}$$ の場合に確率変数 $${X}$$ から得られた観測値が $${x_{0}}$$ 以下となる確率を意味します。
図で表すと以下のようになります。
1.2. 分散が既知の場合の正規分布の平均の信頼区間の計算
本項ではよく例として挙げられる正規分布に従う分布において分散が既知の場合の平均の信頼区間の計算が上式 $${(1.1b)}$$ による計算方法と同じであることを説明します。
$${X}$$ を平均が未知で $${\mu}$$ 、分散が既知で1 の正規分布 $${N(\mu,1)}$$ に従う確率変数とします。
この確率変数 $${X}$$ から独立に $${n}$$ 個の観測値 $${x_{01},\ ...\ ,x_{0n}}$$ が得られたとします。
このとき、$${n}$$ 個の観測値から平均 $${\mu}$$ の信頼区間を求めるために、$${n}$$ 個の観測値の平均を表す確率変数を $${\bar{X}}$$ として、 $${\dfrac{\bar{X} -\mu }{1/\sqrt{n}}}$$ が標準正規分布 $${N(0,1)}$$ に従うことを利用し、以下のように $${\mu}$$ の $${100(1-\alpha)}$$ % 信頼区間を求めます。
$${z_{\alpha}}$$ を標準正規分布の上側 $${100\alpha}$$ % 点とし、$${\bar{x}_{0}=\dfrac{x_{01}+\cdots+x_{0n}}{n}}$$ とすると、以下の式が導かれます。
$$
\begin{aligned}
-z_{\alpha /2} \leq \sqrt{n}\left( \bar{x}_{0}-\mu\right) \leq z_{\alpha /2}
\end{aligned}
$$
上式を利用して、 $${\mu}$$ の $${100(1-\alpha)}$$ % 信頼区間は
$$
\begin{aligned}
\bar{x}_{0} -\dfrac{z_{\alpha /2}}{\sqrt{n}} \leq \mu \leq \bar{x}_{0} +\dfrac{z_{\alpha /2}}{\sqrt{n}}
\end{aligned}
$$
となります。信頼区間の下限、上限をそれぞれ $${\mu_{LL},\ \mu_{UL}}$$ とすると、以下のように表すことができます。
$$
\begin{aligned}
\mu_{LL}&= \bar{x}_{0} -\dfrac{z_{\alpha /2}}{\sqrt{n}} \\
\mu_{UL}&= \bar{x}_{0} +\dfrac{z_{\alpha /2}}{\sqrt{n}}
\end{aligned}
$$
まず、信頼区間の下限 $${\mu_{LL}}$$ について変形すると、
$$
\begin{aligned}
\sqrt{n}(\mu_{LL}-\bar{x})&= -z_{\alpha /2} \\
\int_{-\sqrt{n}(\mu_{LL}-\bar{x}_{0})}^{\infty}\dfrac{1}{\sqrt{2\pi}}\exp{\left[-\dfrac{z^{2}}{2}\right]}dz
&=\alpha/2\\
\int_{-\infty}^{\sqrt{n}(\mu_{LL}-\bar{x}_{0})}\dfrac{1}{\sqrt{2\pi}}\exp{\left[-\dfrac{z^{2}}{2}\right]}dz
&= \alpha/2\\
\end{aligned}
$$
と表すことができます。次に、以下のように $${z \to \bar{x}}$$ の変数変換します。
$$
\left\{\begin{aligned}
z &= \sqrt{n}(\mu_{LL}-\bar{x}) \\
dz &= -\sqrt{n} \bar{x} \\
\bar{x} &\ :\infty \to \bar{x}_{0}
\end{aligned}\right.
$$
$$
\begin{aligned}
\int_{-\infty}^{\sqrt{n}(\mu_{LL}-\bar{x}_{0})}\dfrac{1}{\sqrt{2\pi}}\exp{\left[-\dfrac{z^{2}}{2}\right]}dz
&=- \int_{\infty}^{\bar{x}_{0}}\dfrac{\sqrt{n}}{\sqrt{2\pi }}\exp{\left[-\dfrac{n(\bar{x}-\mu_{LL})^{2}}{2}\right]}d\bar{x}\\
&=\int_{\bar{x}_{0}}^{\infty}\dfrac{1}{\sqrt{2\pi \cdot \dfrac{1}{n}}}\exp{\left[-\dfrac{(\bar{x}-\mu_{LL})^{2}}{2\cdot \dfrac{1}{n}} \right]}d\bar{x} \\
\end{aligned}
$$
よって、以下の式が導かれます。
$$
\begin{aligned}
\int_{\bar{x}_{0}}^{\infty}\dfrac{1}{\sqrt{2\pi \cdot \dfrac{1}{n}}}\exp{\left[-\dfrac{(\bar{x}-\mu_{LL})^{2}}{2\cdot \dfrac{1}{n}} \right]}d\bar{x}=P(\bar{X} \geq \bar{x}_{0}|\mu_{LL}) = \alpha/2\\
\end{aligned}
$$
続いて、信頼区間の上限 $${\mu_{UL}}$$ について変形すると、
$$
\begin{aligned}
\sqrt{n}(\mu_{UL}-\bar{x})&= z_{\alpha /2} \\
\int_{\sqrt{n}(\mu_{UL}-\bar{x}_{0})}^{\infty}\dfrac{1}{\sqrt{2\pi}}\exp{\left[-\dfrac{z^{2}}{2}\right]}dz
&= \alpha /2\\
\end{aligned}
$$
と表すことができます。$${\mu_{LL}}$$ の場合と同様に $${z \to \bar{x}}$$ の変数変換します。
$$
\left\{\begin{aligned}
z &= \sqrt{n}(\mu_{UL}-\bar{x}) \\
dz &= -\sqrt{n} \bar{x} \\
\bar{x} &\ :\bar{x}_{0} \to -\infty
\end{aligned}\right.
$$
$$
\begin{aligned}
\int_{\sqrt{n}(\mu_{UL}-\bar{x}_{0})}^{\infty}\dfrac{1}{\sqrt{2\pi}}\exp{\left[-\dfrac{z^{2}}{2}\right]}dz
&= -\int_{\bar{x}_{0}}^{-\infty}\dfrac{\sqrt{n}}{\sqrt{2\pi}}\exp{\left[-\dfrac{n(\bar{x}-\mu_{UL})^{2}}{2} \right]}d\bar{x}\\
&= \int_{-\infty}^{\bar{x}_{0}}\dfrac{1}{\sqrt{2\pi \cdot \dfrac{1}{n}}}\exp{\left[-\dfrac{(\bar{x}-\mu_{UL})^{2}}{2\cdot \dfrac{1}{n}} \right]}d\bar{x}\\
\end{aligned}
$$
よって、以下の式が導かれます。
$$
\begin{aligned}
\int_{-\infty}^{\bar{x}_{0}}\dfrac{1}{\sqrt{2\pi \cdot \dfrac{1}{n}}}\exp{\left[-\dfrac{(\bar{x}-\mu_{UL})^{2}}{2\cdot \dfrac{1}{n}} \right]}d\bar{x}=P(\bar{X} \leq \bar{x}_{0}|\mu_{UL})=\alpha/2\\
\end{aligned}
$$
以上より、$${\mu_{LL},\mu_{UL}}$$ について以下の式が得られます。
$$
\left\{\begin{aligned}
P(\bar{X} \geq \bar{x}_{0}|\mu_{LL})&=\alpha/2\\
P(\bar{X} \leq \bar{x}_{0}|\mu_{UL})&=\alpha/2\\
\end{aligned}\right.
$$
縦軸に $${P(\bar{X} \leq \bar{x}_{0}|\mu)}$$ 、横軸に $${\mu}$$ をとったグラフと信頼区間は以下の図のようになります。累積分布関数に等しい $${P(\bar{X} \leq \bar{x}_{0}|\mu)}$$ は $${\mu}$$ に対して単調減少となります。
以上より、分散が既知の正規分布の観測値から平均の信頼区間を求める際には、式 $${(1.1b)}$$ と同様の計算をしていることが分かります。
1.3. Clopper-Pearson 法による二項分布の信頼区間の計算
$${X}$$ を二項分布 $${B(n,p)}$$ に従う確率変数とし、得られた観測値を $${x_{0}}$$ とします。成功確率 $${p}$$ は未知とします。
式 $${(1.1b)}$$ に従って、$${p}$$ の$${100(1-\alpha)}$$ % 信頼区間を求めます。信頼区間の下限、上限を順に $${p_{LL},p_{UL}}$$ とすると、以下のように表せます。
$$
\left\{ \begin{aligned}
\sum_{k=x_{0}}^{n}{}_{n}C_{k} p_{LL}^{k}(1-p_{LL})^{n-k}
&=\alpha/2 \\
\sum_{k=0}^{x_{0}}{}_{n}C_{k} p_{UL}^{k}(1-p_{UL})^{n-k}
&=\alpha/2
\end{aligned}\hspace{25pt} (1.3.1) \right.
$$
この式は Clopper-Pearson の式と同じです。よって Clopper-Pearson 法は 1.2 のように ”分散が既知の正規分布より得られた観測値から平均 $${\mu}$$ の信頼区間を求める” 場合と同様の計算をしていることが分かります。
上式 $${(1.3.1)}$$ は 二項分布の累積分布関数ですが、ベータ関数の累積分布関数を用いて計算することもできます。その場合は以下のような式になります。
$$
\left\{\begin{aligned}
\int_{0}^{p_{LL}} \dfrac{1}{\Beta(x_{0},n-x_{0}+1)} p^{x_{0}-1}(1-p)^{n-x_{0}}dp =F_{\Beta(x_{0},n-x_{0}+1)}(p_{LL}) &=\alpha/2 \\
\int_{0}^{p_{UL}} \dfrac{1}{\Beta(x_{0}+1,n-x_{0})} p^{x_{0}}(1-p)^{n-x_{0}-1}dp =F_{\Beta(x_{0}+1,n-x_{0})}(p_{UL}) &=1-\alpha/2 \\
\end{aligned} \hspace{25pt} (1.3.2) \right.
$$
$${F_{\Beta(a,b)}(p)}$$ は $${Beta(a,b)}$$ に従う分布の累積分布関数を表します。
上式は二項分布の上側確率とベータ分布の累積分布関数に以下の関係があることから求めることができます。
$$
\begin{aligned}
F_{\Beta(x,n-x+1)}(p) &=\sum_{k=x}^{n} {}_{n}C_{k}p^{k}(1-p)^{n-k} \hspace{25pt} (1.3.3)
\end{aligned}
$$
以下、証明です。部分積分を複数回繰り返すことで導出できます。
$$
\begin{aligned}
F_{\Beta(x,n-x+1)}(p)
&=\int_{0}^{p} \dfrac{1}{\Beta(x,n-x+1)} p^{x-1}(1-p)^{n-x}dp \\
&=\int_{0}^{p} \dfrac{n!}{(n-x)!(x-1)!} p^{x-1}(1-p)^{n-x}dp \\
&=\int_{0}^{p} \dfrac{n!}{(n-x)!x!} (p^{x})'(1-p)^{n-x}dp \\
&= {}_{n}C_{x}p^{x}(1-p)^{n-x}+\int_{0}^{p} \dfrac{n!}{(n-x)!x!}(n-x)p^{x}(1-p)^{n-x-1}dp \\
&= {}_{n}C_{x}p^{x}(1-p)^{n-x}+\int_{0}^{p} \dfrac{n!}{(n-x-1)!(x+1)!}(p^{x+1})'(1-p)^{n-x-1}dp \\
&= \cdots \\
&=\sum_{k=x}^{n} {}_{n}C_{k}p^{k}(1-p)^{n-k}
\end{aligned}
$$
2. 正規近似による二項分布の信頼区間の計算方法１
二項分布の正規近似による信頼区間の計算方法１について説明します。
$${X}$$ を二項分布 $${B(n,p)}$$ に従う確率変数として、$${X}$$ を正規分布 $${N \left(np,\sqrt{np(1-p)} \right)}$$ と近似することにより求めます。
得られた観測値を $${x_{0}}$$ 、$${z_{\alpha}}$$ を標準正規分布の上側 $${100\alpha}$$ % 点 として、 $${p}$$ の $${100(1-\alpha)}$$ % 信頼区間を求めると、
$$
\begin{aligned}
-z_{\alpha/2} \leq \dfrac{x_{0}-np}{\sqrt{np(1-p)}}& \leq z_{\alpha/2} \\
\dfrac{(x_{0}-np)^{2}}{np(1-p)}& \leq z^{2}_{\alpha/2} \\
x_{0}^{2} -2npx_{0} + n^{2}p^{2} \leq &z^{2}_{\alpha/2}np - z^{2}_{\alpha/2}np^{2} \\
n(n+z^{2}_{\alpha/2})p^{2}-n(2x_{0}+&z_{\alpha/2}^{2})p+x_{0}^{2} \leq 0\\
\end{aligned}
$$
以上より、２次方程式の解より $${p_{LL},p_{UL}}$$ は以下のようになります。
$$
\left\{\begin{aligned}
p_{LL}&=\dfrac{n(2x_{0}+z_{\alpha/2}^{2})-\sqrt{D}}{2n(n+z_{\alpha/2}^{2})} \\
p_{UL}&=\dfrac{n(2x_{0}+z_{\alpha/2}^{2})+\sqrt{D}}{2n(n+z_{\alpha/2}^{2})} \\
D&=n^{2}(2x_{0}+z_{\alpha/2}^{2})^{2}-4nx_{0}^{2}(n+z_{\alpha/2}^{2}) \\
&=nz_{\alpha/2}^{2}(nz_{\alpha/2}^{2}+4x_{0}(n-x_{0}))
\end{aligned} \hspace{15pt}(2.1)\right.
$$
3. 正規近似による二項分布の信頼区間の計算方法２
二項分布の正規近似による信頼区間の計算方法２について説明します。
$${X}$$ を二項分布 $${B(n,p)}$$ に従う確率変数、得られた観測値を $${x_{0}}$$ として、$${X}$$ を正規分布 $${N \left(np,\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \right)}$$ と近似することにより求めます。
正規近似による計算方法１において正規近似した分散の $${p}$$ を $${\dfrac{x_{0}}{n}}$$ に置き換えて計算を簡略化したものです。
$${z_{\alpha}}$$ を標準正規分布の上側 $${100\alpha}$$ % 点 として、 $${p}$$ の $${100(1-\alpha)}$$ % 信頼区間を求めると、
$$
\begin{aligned}
-z_{\alpha/2} \leq \dfrac{x_{0}-np}{\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)}}& \leq z_{\alpha/2} \\
-z_{\alpha/2}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \leq x_{0}&-np \leq z_{\alpha/2}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
\dfrac{x_{0}}{n}-\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \leq &p \leq \dfrac{x_{0}}{n}+\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
\end{aligned}
$$
以上より、$${p_{LL},p_{UL}}$$ は以下のようになります。
$$
\left\{\begin{aligned}
p_{LL}&=\dfrac{x_{0}}{n}-\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
p_{UL}&=\dfrac{x_{0}}{n}+\dfrac{z_{\alpha/2}}{n}\sqrt{x_{0}\left(1-\dfrac{x_{0}}{n}\right)} \\
\end{aligned} \hspace{15pt} (3.1)\right.
$$
4. Python シミュレーションによる検証
前置き
本章では、pythonを使用したシミュレーションにより、以下の３方法による信頼区間の精度、すなわち信頼区間内にどの程度真値が含まれるかを比較します。
方法１：Clopper-Pearson 法を用いた信頼区間の計算方法：下限、上限 $${p_{LL1},p_{UL1}}$$
方法２：2章の正規近似１による信頼区間の計算方法：下限、上限 $${p_{LL2},p_{UL2}}$$
方法３：3章の正規近似２による信頼区間の計算方法：下限、上限 $${p_{LL3},p_{UL3}}$$
比較方法は以下の手順で実施します。
二項分布に従う乱数を生成
得られた乱数（観測値）から、各方法において未知パラメータ $${p}$$ の信頼区間を計算
上記の試行を100000回繰り返し、計算した信頼区間について真値が含まれる確率（どの程度真値が含まれるか）により各方法を評価
結果のグラフ
$${n}$$ を $${5}$$ から $${10000}$$ まで変化させたときの信頼区間に真値が含まれる確率を以下に示します。
条件：$${n=5,\ 10,\ 20,\ 50,\ ...\ ,10000}$$ 、$${p=0.5}$$ 、信頼係数 $${=0.95}$$ 、生成乱数の数（反復回数） $${N=100000}$$
まとめ
以上の結果から、Clopper-Pearson 法を用いた方が、信頼区間が真値に含まれる確率が一番高くなります。真値に含まれる確率が信頼係数に近いという観点では、正規近似による方法１の方が望ましいという結果になります。
5. Python ソースコード
以下、簡易シミュレーションで使用したPythonのソースコードです。
numpy, pandas, scipy は標準ライブラリにありませんので、インストールが必要になります。
バージョン情報
python : 3.12.0
numpy : 1.26.2
pandas : 2.1.3
scipy : 1.13.0
import numpy as np
import scipy
import pandas as pd
import os
#Clopper Pearson 法による信頼区間の計算
def CALC_CI_CLOPPER_PEARSON(n,x,q_ul,q_ll):
p_ul = scipy.stats.beta.ppf(q_ul, x+1, n-x)
p_ll = scipy.stats.beta.ppf(q_ll, x, n-x+1)
return p_ul,p_ll
#正規近似による信頼区間の計算１
def CALC_CI_NORM_APPROX1(n,x,p_ul,p_ll):
a = n * (n + p_ul **2)
b = -n* (2*x+p_ul**2)
D = n * p_ul**2 * (n *p_ul**2+ 4*x*(n-x))
p_ul_bin = (-b + np.sqrt(D)) / (2*a)
p_ll_bin = (-b - np.sqrt(D)) / (2*a)
return p_ul_bin,p_ll_bin
#正規近似による信頼区間の計算２
def CALC_CI_NORM_APPROX2(n,x,p_ul,p_ll):
a = x/n
b = p_ul/n*np.sqrt(x*(1-a))
p_ul_bin = a+b
p_ll_bin = a-b
return p_ul_bin, p_ll_bin
#メイン関数
def main():
#### 定数の設定 ####
n = 10000 # 試行回数
p = 0.5 # 成功確率
size = 100000 # 生成乱数の数
q = 0.95 #信頼係数
####################
# 二項分布に従う乱数の生成
x_arr = np.random.binomial(n, p, size)
# 累積確率の計算
q_ul = (1 + q)/2 #上限
q_ll = (1 - q)/2 #下限
# 上限のパーセント点を計算（正規近似）
p_ul_norm = scipy.stats.norm.ppf(q_ul, loc=0, scale=1)
# 下限のパーセント点を計算（正規近似）
p_ll_norm = scipy.stats.norm.ppf(q_ll, loc=0, scale=1)
# 計算結果格納用の配列を作成
p_ul1_arr =[] #Clopper-Pearson 信頼区間上限用
p_ll1_arr =[] #Clopper-Pearson 信頼区間下限用
withinbounds_count1 = 0 #信頼領域内判定カウント
p_ul2_arr =[] #正規近似1 信頼区間上限用
p_ll2_arr =[] #正規近似1 信頼区間下限用
withinbounds_count2 = 0 #信頼領域内判定カウント
p_ul3_arr =[] #正規近似2 信頼区間上限用
p_ll3_arr =[] #正規近似2 信頼区間下限用
withinbounds_count3 = 0 #信頼領域内判定カウント
# 生成した乱数から信頼区間の計算
for x in x_arr:
# Clopper-Pearson 信頼区間の計算
p_ul1,p_ll1=CALC_CI_CLOPPER_PEARSON(n,x,q_ul,q_ll)
#領域内判定とカウント
if p_ll1 <= p <= p_ul1:
withinbounds_count1 +=1
#計算結果を追加
p_ll1_arr.append(p_ll1)
p_ul1_arr.append(p_ul1)
# 正規近似1 信頼区間の計算
p_ul2,p_ll2 = CALC_CI_NORM_APPROX1(n,x,p_ul_norm,p_ll_norm)
#領域内判定とカウント
if p_ll2 <= p <= p_ul2:
withinbounds_count2 +=1
#計算結果を追加
p_ll2_arr.append(p_ll2)
p_ul2_arr.append(p_ul2)
# 正規近似2 信頼区間の計算
p_ul3,p_ll3 = CALC_CI_NORM_APPROX2(n,x,p_ul_norm,p_ll_norm)
#領域内判定とカウント
if p_ll3 <= p <= p_ul3:
withinbounds_count3 +=1
#計算結果を追加
p_ll3_arr.append(p_ll3)
p_ul3_arr.append(p_ul3)
#リスト⇒二次元配列の変換 （列接続のため
p_ll1_arr = np.array(p_ll1_arr).reshape(size,1)
p_ul1_arr = np.array(p_ul1_arr).reshape(size,1)
p_ll2_arr = np.array(p_ll2_arr).reshape(size,1)
p_ul2_arr = np.array(p_ul2_arr).reshape(size,1)
p_ll3_arr = np.array(p_ll3_arr).reshape(size,1)
p_ul3_arr = np.array(p_ul3_arr).reshape(size,1)
x_arr = np.array(x_arr).reshape(size,1)
#データの接続、データ列の作成
outdata = np.concatenate((x_arr, p_ll1_arr,p_ul1_arr, p_ll2_arr, p_ul2_arr, p_ll3_arr, p_ul3_arr), axis=1)
# ファイル出力
## 表題作成
column1 = ['観測値',
'Clopper-Pearson：信頼区間下限', 'Clopper-Pearson：信頼区間上限',
'正規近似1：信頼区間下限', '正規近似1：信頼区間上限',
'正規近似2：信頼区間下限', '正規近似2：信頼区間上限'
]
pd_out = pd.DataFrame(outdata, columns= column1)
# 名前の設定
outname0 = f'二項分布の信頼区間計算_n={n}_p={p}_信頼係数={str(q).replace(".","p")}_試行回数={size}'
outname = outname0 + '.csv'
#同名ファイルが存在する場合の名前変更
filenum = 1
while( os.path.exists(outname) ):
outname = outname0 + f'_{filenum}.csv'
filenum += 1
# csvデータ出力 #shift-jis でないとエクセルで文字化け
pd.DataFrame(pd_out).to_csv(outname,header=True, encoding='shift-jis')
print(f'ファイルを保存しました。ファイル名:{outname}')
#print(x_arr)
#print(f'{results}')
print(f'----真値が信頼領域に含まれるかの判定-----\n'\
+f'Clopper-Pearson：{withinbounds_count1}回\n'\
+f' 正規近似1：{withinbounds_count2}回\n'\
+f' 正規近似2：{withinbounds_count3}回\n'\
+f'-----------------------------------------'\
)
if __name__ == "__main__":
while(1):
main()
input('ENTERキーで最初に戻ります。')