からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【Python】3.4.4:二項分布の母比率の信頼区間の実装【統計検定2級のノート】

はじめに

 「統計検定2級」の独学時のまとめノートです。各種の統計手法について「数式・プログラム・図」を用いて解説します。
 『統計学基礎』に沿って学習を進めます。本の補助として読んでください。

 この記事では、二項分布の成功確率パラメータの信頼区間をPythonでスクラッチ実装して、計算過程を可視化します。

【前節の内容】

www.anarchive-beta.com

【この節の内容】

3.4.4 二項分布の母比率の信頼区間の実装

 二項分布(Binomial distribution)の母比率(population proportion)に関する信頼区間(confidence interval)を実装して、人工的に生成したデータを用いて、パラメータの区間推定を行います。
 二項分布については「二項分布の定義式 - からっぽのしょこ」、正規分布については「1次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

 利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
from scipy.stats import binom, norm
import matplotlib.pyplot as plt


区間推定の実装

 まずは、二項分布の母比率に関する信頼区間の計算における一連の処理を確認します。
 母分布(二項分布)を設定して、母分布に従う標本(トイデータ)を生成します。続いて、生成した標本を用いて、母比率の信頼区間の計算(パラメータの区間推定)を行います。

母集団の設定

 母集団のデータが従う分布(母分布・population distribution・二項分布)  x \sim \mathrm{Bin}(N, p_{\mathrm{pop}}) の母数(パラメータ)  N, p_{\mathrm{pop}} を設定します。
 二項分布のグラフ作成については「【Python】二項分布の作図 - からっぽのしょこ」を参照してください。

 母分布のパラメータ  N, p_{\mathrm{pop}} を設定します。

# サンプルサイズを指定
N = 20
print(N)

# 母比率を指定
p_pop = 0.6
print(p_pop)
20
0.6

 二項分布の試行回数(サンプルサイズ・正の整数)  N、成功確率(母比率・0から1の値)  0 \leq p_{\mathrm{pop}} \leq 1 を指定します。

  p_{\mathrm{pop}} が区間推定における真の値であり、この値が含まれる区間を求めるのがここでの目的(推定・推測)です。

 母分布の確率変数  x の作図範囲を設定します。

# x軸の範囲を設定
x_lower = 0
x_upper = N
x_min   = x_lower - 0.5
x_max   = x_upper + 0.5
print(x_min, x_max)

# x軸の値を作成
x_vec = np.arange(start=x_lower, stop=x_upper+1, step=1)
print(x_vec[:5])
-0.5 20.5
[0 1 2 3 4]

 この例では、指定したパラメータを使って、 x 軸の範囲を設定しています。

 母分布の確率を計算します。

# 母分布の確率を計算
pop_prob_vec = binom.pmf(k=x_vec, n=N, p=p_pop)
print(pop_prob_vec[:5])
[1.09951163e-08 3.29853488e-07 4.70041221e-06 4.23037099e-05
 2.69686150e-04]

  x 軸の値ごとに、二項分布に従う確率  \mathrm{Bin}(x \mid N, p_{\mathrm{pop}}) を計算します。
 二項分布の確率質量関数 sp.stats.binom.pmf() の確率変数の引数 k x、試行回数の引数 n N、成功確率の引数 p p_{\mathrm{pop}} を指定します。

 母分布のグラフを作成します。

・作図コード(クリックで展開)

# 凡例の表示位置を設定
legend_loc = 'upper left' if p_pop >= 0.5 else 'upper right'
# 母分布のラベルを作成
pop_param_lbl  = f'$N = {N}, '
pop_param_lbl += 'p_{pop} = '+f'{p_pop:.2f}$'

# 母分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.plot(
    x_vec, pop_prob_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
ax.scatter(
    x=x_vec, y=pop_prob_vec, 
    color='black', s=50, 
    zorder=10
) # 母分布
ax.axvline(
    x=N*p_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population proportion', 
    zorder=11
) # 母比率

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[N*p_pop], 
    labels=['$N p_{pop}$']
) # パラメータのラベル
ax.set_ylabel('$Bin(x \\mid n, p_{pop})$')
ax.set_title(pop_param_lbl, loc='left')
ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

plt.show()

母集団のデータが従う分布(二項分布)のグラフ

 母分布(二項分布)を曲線と点、母比率(成功確率パラメータ)に対応する位置(サンプルサイズ倍の値)を垂直線(赤色・破線)で示します。

 試行回数  N・成功確率  p_{\mathrm{pop}} (失敗確率  1-p_{\mathrm{pop}} )の二項分布  \mathrm{Bin}(x \mid N, p_{\mathrm{pop}}) について、真の値  p_{\mathrm{pop}} を区間で求めていきます。

標本の生成

 設定した母分布(二項分布)  x \sim \mathrm{Bin}(N, p_{\mathrm{pop}}) に従う標本(サンプル)  x を作成します。
 二項分布のサンプリングや集計については「【Python】二項分布の作図 - からっぽのしょこ」を参照してください。

 母分布から標本  x を生成します。

# 標本を生成
x_obs, = np.random.binomial(n=N, p=p_pop, size=1)
print(x_obs)
10

 二項分布に従う乱数(実現値)  x を生成します。
 二項分布の乱数生成関数 np.random.binomial() の試行回数の引数 n N、成功確率の引数 p p_{\mathrm{pop}}、サンプルサイズの引数 size 1 を指定します。

 標本  x の標本統計量  \hat{p} を計算します。

# 標本比率を計算
p_hat = x_obs / N
print(p_hat)
0.5

 標本  x を用いて、標本比率  \hat{p} を、次の式で計算します。

 \displaystyle
\hat{p}
    = \frac{x}{N}

 標本  x のグラフを作成します。

・作図コード(クリックで展開)

# 標本のラベルを作成
obs_param_lbl  = '$x_{obs} = '+f'{x_obs}, '
obs_param_lbl += '\\hat{p} = \\frac{x_{obs}}{N} = '+f'{p_hat:.2f}$'

# 標本を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.plot(
    x_vec, pop_prob_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
ax.scatter(
    x=x_vec, y=pop_prob_vec, 
    color='black', s=50, 
    zorder=10
) # 母分布
ax.scatter(
    x=x_obs, y=0.0, 
    color='black', alpha=0.5, s=100, 
    label='sample value', 
    zorder=11
) # 標本
ax.axvline(
    x=x_obs, 
    color='black', linewidth=1.0, linestyle='--', 
    zorder=12
) # 標本

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[x_obs], 
    labels=['$x_{obs}$']
) # 標本のラベル
ax.set_ylabel('$Bin(x \\mid n, p_{pop})$')
ax.set_title(obs_param_lbl, loc='left')
ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

plt.show()

標本(二項乱数)のグラフ

 標本(二項分布の乱数)の実現値を点と垂直線(破線)で示します。

  N は、二項分布のパラメータであり、標本の数には対応しません。

標本分布の計算

 標本比率  \hat{p} を確率変数とする分布(標本分布・sampling distribution・正規分布)  \hat{p} \sim \mathcal{N}(\mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}^2) の分散パラメータ  \sigma_{\mathrm{smp}}^2 を求めます。平均パラメータ  \mu_{\mathrm{smp}} は未知です。
 正規分布のグラフ作成については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

 標本分布のパラメータ  \mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}^2 を計算します。

# 標本分布の平均を計算
mu_smp = p_pop
print(mu_smp)

# 標本分布の分散を計算
sigma2_smp = p_hat * (1.0-p_hat) / N
print(sigma2_smp)

# 標本分布の標準偏差を計算
sigma_smp = np.sqrt(sigma2_smp)
print(sigma_smp)
0.6
0.0125
0.11180339887498948

 母数  N, p_{\mathrm{pop}} と標本  x (から求めた統計量)を用いて、標本分布の平均  \mu_{\mathrm{smp}}、分散  \sigma_{\mathrm{smp}}^2、標準偏差  \sigma_{\mathrm{smp}} を、次の式で計算します。

 \displaystyle
\begin{aligned}
\mu_{\mathrm{smp}}
   &= p_{\mathrm{pop}}
\\
\sigma_{\mathrm{smp}}^2
   &= \frac{
          p_{\mathrm{pop}}
          (1 - p_{\mathrm{pop}})
      }{
          N
      }
\\
   &\approx
      \frac{
          \hat{p}
          (1 - \hat{p})
      }{
          N
      }
\\
\sigma_{\mathrm{smp}}
   &= \sqrt{\sigma_{\mathrm{smp}}^2}
\\
   &= \sqrt{
          \frac{
              p_{\mathrm{pop}}
              (1 - p_{\mathrm{pop}})
          }{
              N
          }
      }
\\
   &\approx
      \sqrt{
          \frac{
              \hat{p}
              (1 - \hat{p})
          }{
              N
          }
      }
\end{aligned}

 母比率  p_{\mathrm{pop}} は未知なので、 \mu_{\mathrm{smp}} も未知の値です。そのため推定には使いませんが、作図(真の値の可視化)用に作成しておきます。

 標本分布の確率変数  \hat{p} の作図範囲を設定します。

# p軸の範囲を設定
p_min = x_min / N
p_max = x_max / N
print(p_min, p_max)

# p軸の値を作成
p_vec = np.linspace(start=p_min, stop=p_max, num=1001)
print(p_vec[:5])
-0.025 1.025
[-0.025   -0.02395 -0.0229  -0.02185 -0.0208 ]

 グラフ間で対応させるため、 x 軸の範囲を正規化して、 p 軸の範囲を設定しています。

 標本分布の確率密度を計算します。

# 標本分布の確率密度を計算
smp_dens_vec = norm.pdf(x=p_vec, loc=mu_smp, scale=sigma_smp)
print(smp_dens_vec[:5])
[5.84256805e-07 6.15722590e-07 6.48825770e-07 6.83648380e-07
 7.20276396e-07]

  p 軸の値ごとに、正規分布に従う確率密度  \mathcal{N}(\hat{p} \mid \mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}^2) を計算します。
 正規分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x p、平均の引数 loc \mu_{\mathrm{smp}}、標準偏差の引数 scale \sigma_{\mathrm{smp}} を指定します。

 標本分布のグラフを作成します。

・作図コード(クリックで展開)

# 標本分布のラベルを作成
smp_param_lbl  = '$\\hat{p}_{obs} = '+f'{p_hat:.2f}, '
smp_param_lbl += '\\mu_{smp} = p_{pop} = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} \\approx \\sqrt{\\frac{\\hat{p} (1-\\hat{p})}{N}} = '+f'{sigma_smp:.2f}$'

# 標本分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Normal distribution', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.plot(
    p_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sampling distribution', 
    zorder=10
) # 標本分布
ax.scatter(
    x=p_hat, y=0.0, 
    color='#00A968', s=100, 
    zorder=11
) # 標本比率
ax.axvline(
    x=p_hat, 
    color='black', linewidth=1.0, linestyle='--', 
    label='sample proportion\nsample value', 
    zorder=12
) # 標本比率

ax.set_xlabel('$\\hat{p} = \\frac{x}{n}$')
ax2x.set_xticks(
    ticks =[p_hat], 
    labels=['$\\hat{p}_{obs}$']
) # 標本のラベル
ax.set_ylabel('$N(\\hat{p} \\mid \\mu_{smp}, \\sigma_{smp}^2)$')
ax.set_title(smp_param_lbl, loc='left')
ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)

plt.show()

標本比率が従う分布(正規分布)のグラフ

 標本分布(正規分布)を曲線、標本比率(の標本)を点(青緑色)と垂直線(破線)で示します。

 標本比率  \hat{p} の理論分布は、中心極限定理により、平均  p_{\mathrm{pop}}・分散  \frac{p_{\mathrm{pop}} (1-p_{\mathrm{pop}})}{N} (標準偏差  \sqrt{\frac{p_{\mathrm{pop}} (1-p_{\mathrm{pop}})}{N}} )の正規分布  \mathcal{N}(\hat{p} \mid p_{\mathrm{pop}}, \frac{p_{\mathrm{pop}} (1-p_{\mathrm{pop}})}{N}) になります。ただし、母比率  p_{\mathrm{pop}} が未知なので、変わりに標本比率  \hat{p} を用いて、分散  \frac{\hat{p} (1-\hat{p})}{N} (標準偏差  \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} )として近似的に計算します。平均も本来は分からない値です(真の値の可視化にのみ使っています)。
 サンプルサイズ  N が大きいほど、標本分布の分散  \sigma_{\mathrm{smp}}^2 \approx \frac{\hat{p} (1-\hat{p})}{N} (標準偏差  \sigma_{\mathrm{smp}} \approx \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} )が小さくなります。標本分布の平均  \mu_{\mathrm{smp}} = p_{\mathrm{pop}} は、 N の影響を受けません。

 サンプルサイズ  N が十分に大きいと、標本比率  \hat{p} が母比率  p_{\mathrm{pop}} に近付きます。

標準化統計量の分布の計算

 標本比率  \hat{p} を標準化した  z を確率変数とする分布(標準化統計量の分布・distribution of a standardized statistic・標準正規分布)  z \sim \mathcal{N}(0, 1) における中央領域  [z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}}] = z_{1-\frac{\alpha}{2}} \leq z \leq z_{\frac{\alpha}{2}} の下限・上限  z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}} を設定します。

 標本比率(の標本)  \hat{p} を標準化します。

# 標本統計量を標準化
z_obs = (p_hat - mu_smp) / sigma_smp
print(z_obs)
-0.8944271909999157

 標本平均  \bar{x} を用いて、標準化した変数  z を、次の式で計算します。

 \displaystyle
\begin{aligned}
z  &= \frac{
           \hat{p} - \mu_{\mathrm{smp}}
       }{
           \sigma_{\mathrm{smp}}
       }
\\
   &= (\hat{p} - p_{\mathrm{pop}})
      \sqrt{
          \frac{
              N
          }{
              \hat{p}
              (1 - \hat{p})
          }
      }
\end{aligned}

 母比率  p_{\mathrm{pop}} は未知なので、 z_{\mathrm{obs}} も未知の値です。そのため推定には使いませんが、作図(真の値の可視化)用に作成しておきます。

 有意水準  \alpha を設定します。

# 信頼係数を指定
gamma = 0.95

# 有意水準を計算
alpha = 1.0 - gamma
print(alpha)
0.050000000000000044

 信頼係数(0から1の値)  \gamma = 1-\alpha を指定して、有意水準(0から1の値)  \alpha = 1-\gamma を計算します。

 中央領域の範囲  [z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}}] を計算します。

# 中央領域の範囲を計算
cr_bound_upper = norm.ppf(q=1.0-0.5*alpha, loc=0.0, scale=1.0)
cr_bound_lower = -cr_bound_upper
print(cr_bound_lower, cr_bound_upper)
-1.959963984540054 1.959963984540054

 連続型の確率分布に関して、確率密度関数を  f(x) として、累積分布関数は

 \displaystyle
F(x)
    = P(X \leq x)
    = \int_{-\infty}^x
          f(u)
      \mathrm{d} u

で定義され、分位点関数はその逆関数  F^{-1}(x) で定義されます。
 パーセンタイル(百分位数)を  p、クオンタイル(分位数)を  q として、累積分布関数と分位点関数は

 \displaystyle
\begin{aligned}
p  &= F(q)
\\
q  &= F^{-1}(p)
\end{aligned}

の関係です。

 標準正規分布の分位点関数  F^{-1}(x) を用いて、上側確率が  p (下側確率が  1-p )となる分位点を  z_p で表すことにします(下側確率を基準にする場合は  z_{1-p} と表記する方が自然な表現だと思われます)。

 \displaystyle
z_{p}
    = F^{-1}(1-p)

  100 (1-\alpha) \% 中央領域の境界値  z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}} を、次の式で計算します。

 \displaystyle
\begin{aligned}
z_{1-\frac{\alpha}{2}}
   &= F^{-1} \Bigl(
          \frac{\alpha}{2}
      \Bigr)
\\
z_{\frac{\alpha}{2}}
   &= F^{-1} \Bigl(
          1-\frac{\alpha}{2}
      \Bigr)
\end{aligned}

 ただし、標準正規分布の場合、 z_{1-\frac{\alpha}{2}} = -z_{\frac{\alpha}{2}} です。

 正規分布の分位点関数 sp.stats.norm.ppd() のパーセンタイルの引数 q 1-\frac{\alpha}{2}、平均の引数 loc 0、標準偏差の引数 scale 1 を指定します。

 標準化統計量の分布の確率変数  z の作図範囲を設定します。

# z軸の範囲を設定
z_min = (p_min - p_hat) / sigma_smp
z_max = (p_max - p_hat) / sigma_smp
print(z_min, z_max)

# z軸の値を作成
z_vec = np.linspace(start=z_min, stop=z_max, num=1001)
print(z_vec[:5])
-4.695742752749559 4.695742752749558
[-4.69574275 -4.68635127 -4.67695978 -4.6675683  -4.65817681]

 グラフ間の対応させるため(真の値  \mu_{\mathrm{smp}} ではなく観測値  \hat{p}_{\mathrm{obs}} を使って)、 p 軸の範囲を標準化して、 z 軸の範囲を設定しています。

 標準化統計量の分布の確率密度を計算します。

# 標準化分布の確率密度を計算
std_dens_vec = norm.pdf(x=z_vec, loc=0.0, scale=1.0)
print(std_dens_vec[:5])
[6.49850418e-06 6.79120182e-06 7.09645685e-06 7.41477866e-06
 7.74669596e-06]

 「標本分布の計算」のときと同様にして、正規分布に従う確率密度  \mathcal{N}(z \mid 0, 1) を計算します。

 標準化統計量の分布のグラフを作成します。

・作図コード(クリックで展開)

# 中央領域を計算
cr_z_vec    = np.linspace(start=cr_bound_lower, stop=cr_bound_upper, num=501)
cr_dens_vec = norm.pdf(x=cr_z_vec, loc=0.0, scale=1.0)

# 外側領域を計算
tail_z_vec    = np.hstack([
    np.linspace(start=z_min, stop=cr_bound_lower, num=251), 
    np.nan, # (塗りつぶしの分割用)
    np.linspace(start=cr_bound_upper, stop=z_max, num=251)
])
tail_dens_vec = norm.pdf(x=tail_z_vec, loc=0.0, scale=1.0)
# 標準化分布のラベルを作成
std_param_lbl  = '$z_{obs} = '+f'{z_obs:.2f}, '
std_param_lbl += '\\mu_{std} = '+f'{0:.2f}, '
std_param_lbl += '\\sigma_{std} = '+f'{1:.2f}$\n'
std_param_lbl += f'$\\alpha = {alpha:.2f}, '
std_param_lbl += 'z_{1-\\frac{\\alpha}{2}} = '+f'{cr_bound_lower:.2f}, '
std_param_lbl += 'z_{\\frac{\\alpha}{2}} = '+f'{cr_bound_upper:.2f}$'

# 標準化分布を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('standard Normal distribution', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.plot(
    z_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='standardized statistic distribution', 
    zorder=10
) # 標準化分布
ax.fill_between(
    x=cr_z_vec, y1=np.zeros_like(cr_z_vec), y2=cr_dens_vec, 
    facecolor='purple', alpha=0.5, 
    label='central region', 
    zorder=11
) # 中央領域
ax.fill_between(
    x=tail_z_vec, y1=np.zeros_like(tail_z_vec), y2=tail_dens_vec, 
    facecolor='blue', alpha=0.5, 
    label='tail regions', 
    zorder=12
) # 外側領域
ax.hlines(
    y=0.0, xmin=cr_bound_lower, xmax=cr_bound_upper, 
    color='purple', linewidth=2.0, 
    zorder=13
) # 中央領域
ax.axvline(
    x=z_obs, 
    color='black', linewidth=1.0, linestyle='--', 
    zorder=14
) # 標準化変数
for idx, z in enumerate([cr_bound_lower, cr_bound_upper]):
    ax.axvline(
        x=z, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='central bounds' if idx == 0 else None, 
        zorder=15
    ) # 中央領域の境界値

ax.set_xlabel('$z = \\frac{\\hat{p} - \\mu_{smp}}{\\sigma_{smp}}$')
ax2x.set_xticks(
    ticks =[z_obs, cr_bound_lower, cr_bound_upper], 
    labels=['$z_{obs}$', '$z_{1-\\frac{\\alpha}{2}}$', '$z_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$N(z \\mid 0, 1)$')
ax.set_title(std_param_lbl, loc='left')
ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=z_min, xmax=z_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=z_min, xmax=z_max) # (目盛の共通化用)

plt.show()

標準化した標本比率が従う分布(正規分布)のグラフ

 標準化統計量の分布(標準正規分布)を曲線、標準化した標本比率(の標本)を垂直線(破線)、境界値を垂直線(鎖線)、中央領域(境界値の内側)を線分と塗りつぶし(紫色)、外側の領域(境界値の外側)を塗りつぶし(青色)で示します。

 標準化変数  z の理論分布は、標準化により、平均  0・分散  1 (標準偏差  1 )の正規分布(標準正規分布)  \mathcal{N}(z \mid 0, 1) になります。
 標準正規分布は、中心が0で左右対称な形状なので、 z_{1-\frac{\alpha}{2}} = -z_{\frac{\alpha}{2}} であり、 |z_{1-\frac{\alpha}{2}}| = z_{\frac{\alpha}{2}} です。

 有意水準  \alpha が小さい(0に近い)ほど、信頼係数  1-\alpha が大きくなり(1に近付き)ます。
  1-\alpha は、中央領域の範囲(面積)

 \displaystyle
\begin{aligned}
P(z_{1-\frac{\alpha}{2}} \leq z \leq z_{\frac{\alpha}{2}})
   &= \int_{z_{1-\frac{\alpha}{2}}}^{z_{\frac{\alpha}{2}}}
          \mathcal{N}(z \mid 0, 1)
      \mathrm{d} z
\\
   &= 1 - \alpha
\end{aligned}

に対応します。
 また、 \alpha は、外側の領域の範囲(面積)

 \displaystyle
\begin{aligned}
P(z \lt z_{1-\frac{\alpha}{2}})+ P(z \gt z_{\frac{\alpha}{2}}) 
   &= 1 - P(z_{1-\frac{\alpha}{2}} \leq z \leq z_{\frac{\alpha}{2}})
\\
   &= \alpha
\end{aligned}

に対応し、左右の片側の面積は

 \displaystyle
\begin{aligned}
P(z \lt z_{1-\frac{\alpha}{2}}) 
   &= \int_{-\infty}^{z_{1-\frac{\alpha}{2}}}
          \mathcal{N}(z \mid 0, 1)
      \mathrm{d} z
\\
   &= \frac{\alpha}{2}
\\
P(z \gt z_{\frac{\alpha}{2}}) 
   &= \int_{z_{\frac{\alpha}{2}}}^{\infty}
          \mathcal{N}(z \mid 0, 1)
      \mathrm{d} z
\\
   &= \frac{\alpha}{2}
\end{aligned}

です。

 標本  x_{\mathrm{obs}} から求めた標本比率  \hat{p}_{\mathrm{obs}} を標準化した確率変数(標本)  z_{\mathrm{obs}} が中央領域(紫色の範囲)  [z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}}] に含まれる確率  P(z_{1-\frac{\alpha}{2}} \leq z_{\mathrm{obs}} \leq z_{\frac{\alpha}{2}}) が信頼係数  1-\alpha に対応します。言い換えると、外側の領域(青色の範囲)  [-\infty, z_{1-\frac{\alpha}{2}}], [z_{\frac{\alpha}{2}}, \infty] に含まれる(中央領域に含まれない)確率  P(z_{\mathrm{obs}} \lt z_{1-\frac{\alpha}{2}}) + P(z_{\mathrm{obs}} \gt z_{\frac{\alpha}{2}}) が有意水準  \alpha に対応します。

 \displaystyle
\begin{aligned}
P(|z_{\mathrm{obs}}| \leq z_{\frac{\alpha}{2}})
   &= P(z_{1-\frac{\alpha}{2}} \leq z_{\mathrm{obs}} \leq z_{\frac{\alpha}{2}})
\\
   &= 1 - \alpha
\\
P(|z_{\mathrm{obs}}| \gt z_{\frac{\alpha}{2}})
   &= P(z_{\mathrm{obs}} \lt z_{1-\frac{\alpha}{2}})
      + P(z_{\mathrm{obs}} \gt z_{\frac{\alpha}{2}})
\\
   &= \alpha
\end{aligned}

 母比率  p_{\mathrm{pop}} が未知なので、標準化後の値  z_{\mathrm{obs}} は、本来は分からない値です(真の値の可視化にのみ使っています)。しかし、母数に関わらず、 z_{\mathrm{obs}} が従う分布はパラメータが固定(標準正規分布)になるので、 z_{\mathrm{obs}} が含まれる確率は分かります。

 標本比率  \hat{p} や標準化した変数  z は確率変数(標本として得られる実現値)なので、それぞれ対応する分布(理論分布)に従って確率的に値が決まるとみなせます。一方、次に求める信頼区間は、標本として得られた実現値(与えられた値)  \hat{p}_{\mathrm{obs}} z_{\mathrm{obs}} に基づいて確定的に値が決まり(一意に定まり)ます。

信頼区間の計算

 母比率  p_{\mathrm{pop}} の信頼区間  [L, U] = L \leq p \leq U の下限・上限  L, U を求めます。

 信頼区間の範囲  [L, U] を計算します。

# 信頼区間の範囲を計算
ci_bound_lower = p_hat + cr_bound_lower * sigma_smp
ci_bound_upper = p_hat + cr_bound_upper * sigma_smp
print(ci_bound_lower, ci_bound_upper)
0.2808693648558546 0.7191306351441453

 信頼区間の境界値  L, U を、次の式で計算します。

 \displaystyle
\begin{aligned}
L  &= \hat{p}
      + z_{1-\frac{\alpha}{2}}
        \sigma_{\mathrm{smp}}
\\
   &= \hat{p}
      - z_{\frac{\alpha}{2}}
        \sigma_{\mathrm{smp}}
\\
U  &= \hat{p}
      + z_{\frac{\alpha}{2}}
        \sigma_{\mathrm{smp}}
\\
   &= \hat{p}
      - z_{1-\frac{\alpha}{2}}
        \sigma_{\mathrm{smp}}
\end{aligned}

 信頼区間のグラフを作成します。

・作図コード(クリックで展開)

# 信頼区間のラベルを作成
ci_param_lbl  = f'$\\alpha = {alpha:.2f}, '
ci_param_lbl += 'L = \\hat{p} + z_{1-\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_lower:.2f}, '
ci_param_lbl += 'U = \\hat{p} + z_{\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_upper:.2f}$'

# 信頼区間を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Population Proportion Confidence Interval', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=p_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population proportion', 
    zorder=10
) # 母比率
ax.hlines(
    y=0.0, xmin=ci_bound_lower, xmax=ci_bound_upper, 
    color='purple', linewidth=2.0, 
    label='confidence interval', 
    zorder=11
) # 信頼区間
for idx, p in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.axvline(
        x=p, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='confidence bounds' if idx == 0 else None, 
        zorder=12
    ) # 信頼区間の境界値
    
ax.set_xlabel('$p$')
ax2x.set_xticks(
    ticks =[p_pop, ci_bound_lower, ci_bound_upper], 
    labels=['$p_{pop}$', '$L$', '$U$']
) # 信頼区間のラベル
ax.set_yticks(ticks=[0.0], labels=['']) # (目盛の非表示化用)
ax.set_title(ci_param_lbl, loc='left')
ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)

plt.show()

母比率の信頼区間のグラフ

 母比率を垂直線(赤色・破線)、信頼区間の下限・上限を垂直線(鎖線)、信頼区間を線分(紫色)で示します。

 母比率  p_{\mathrm{pop}} の信頼区間は、標本比率  \hat{p} を中心に  \pm z_{\frac{\alpha}{2}} \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} の範囲です。 z_{\frac{\alpha}{2}} \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} は、確率変数  \hat{p} の標準偏差  \sigma_{\mathrm{smp}} \approx \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} z_{\frac{\alpha}{2}} 個分のサイズを表します。
 サンプルサイズ  N が大きいほど、標本分布の標準偏差  \sigma_{\mathrm{smp}} が小さくなり、信頼区間の範囲(幅)  U-L = 2 z_{\frac{\alpha}{2}} \sigma_{\mathrm{smp}} が小さくなります。

 以上で、二項分布の母比率の信頼区間の計算を実装しました。

スポンサードリンク

母集団と信頼区間の対応関係

 次は、母集団やパラメータ・統計量と信頼区間の対応関係を図で確認します。

 母比率の信頼区間の計算のグラフを作成します。

・作図コード(クリックで展開)

# ラベルを作成
pop_param_lbl  = f'$N = {N}, '
pop_param_lbl += 'p_{pop} = '+f'{p_pop:.2f}$'
smp_param_lbl  = '$\\hat{p} = '+f'{p_hat:.2f}$\n'
smp_param_lbl += '$\\mu_{smp} = p_{pop} = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} \\approx \\sqrt{\\frac{\\hat{p} (1-\\hat{p})}{N}} = '+f'{sigma_smp:.2f}$\n'
smp_param_lbl += '$L = \\hat{p} - z_{\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_lower:.2f}, '
smp_param_lbl += 'U = \\hat{p} + z_{\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_upper:.2f}$'
std_param_lbl  = f'$z = {z_obs:.2f}$\n'
std_param_lbl += f'$\\alpha = {alpha:.2f}, '
std_param_lbl += 'z_{1-\\frac{\\alpha}{2}} = '+f'{cr_bound_lower:.2f}, '
std_param_lbl += 'z_{\\frac{\\alpha}{2}} = '+f'{cr_bound_upper:.2f}$'

# 母比率の信頼区間の計算を作図
fig, axes = plt.subplots(
    nrows=3, ncols=1, 
    figsize=(9, 12), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Population Proportion Confidence Interval', fontsize=20)

##### 母分布の作図 -----

# 母分布を描画
ax   = axes[0]
ax2x = ax.twiny()

ax.scatter(
    x=x_obs, y=0.0, 
    color='black', alpha=0.5, s=100, clip_on=False, 
    label='sample value', 
    zorder=9
) # 標本
ax.plot(
    x_vec, pop_prob_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
ax.scatter(
    x=x_vec, y=pop_prob_vec, 
    color='black', s=30, 
    zorder=10
) # 母分布
for idx, x in enumerate([N*p_pop, x_obs]):
    ax.axvline(
        x=x, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母比率・標本
for p in [ci_bound_lower, ci_bound_upper]:
    x = N * p # 逆変換
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[N*p_pop, x_obs, N*ci_bound_lower, N*ci_bound_upper], 
    labels=['$N p_{pop}$', '$x_{obs}$', '$N L$', '$N U$']
) # 信頼区間のラベル
ax.set_ylabel('$Bin(x \\mid n, p_{pop})$')
ax.set_title(pop_param_lbl, loc='left')
#ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

# ラベルの装飾を調整(表示順の変更用)
order = [1, 0] # 表示順を指定
handles, labels = ax.get_legend_handles_labels() # 凡例情報を取得
ax.legend(
    handles=[handles[i] for i in order], 
    labels =[labels[i] for i in order], 
    loc=legend_loc, prop={'size': 8}
)

##### 標本分布の作図 -----

# 標本分布を描画
ax   = axes[1]
ax2x = ax.twiny()

ax.scatter(
    x=p_hat, y=0.0, 
    color='#00A968', s=100, clip_on=False, 
    label='sample', 
    zorder=9
) # 標本比率
ax.plot(
    p_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sampling distribution', 
    zorder=10
) # 標本分布
for idx, p in enumerate([p_pop, p_hat]):
    ax.axvline(
        x=p, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        label=['population proportion', 'saple proportion'][idx], 
        zorder=[20, 21][idx]
    ) # 母比率・標本比率
for idx, p in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.axvline(
        x=p, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='confidence bounds' if idx == 0 else None, 
        zorder=22
    ) # 信頼区間の境界値
ax.hlines(
    y=0.0, xmin=ci_bound_lower, xmax=ci_bound_upper, 
    color='purple', linewidth=2.0, 
    label='confidence interval', 
    zorder=30
) # 信頼区間

ax.set_xlabel('$\\hat{p} = \\frac{x}{n}$')
ax2x.set_xticks(
    ticks =[p_pop, p_hat, ci_bound_lower, ci_bound_upper], 
    labels=[
        '$p_{pop}$', 
        '$\\hat{p}_{obs}$', 
        '$\\hat{p}_{obs} - z_{\\frac{\\alpha}{2}} \\sigma_{smp}$', 
        '$\\hat{p}_{obs} + z_{\\frac{\\alpha}{2}} \\sigma_{smp}$'
    ]
) # 信頼区間のラベル
ax.set_ylabel('$N(\\hat{p} \\mid \\mu_{smp}, \\sigma_{smp}^2)$')
ax.set_title(smp_param_lbl, loc='left')
#ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)

# ラベルの装飾を調整(重なる対策用)
halignments = ['center', 'center', 'right', 'left'] # 表示位置を指定
labels      = ax2x.get_xticklabels() # 軸情報を取得
for label, ha in zip(labels, halignments):
    label.set_ha(ha)

# ラベルの装飾を調整(表示順の変更用)
order = [1, 2, 3, 4, 5, 0] # 表示順を指定
handles, labels = ax.get_legend_handles_labels() # 凡例情報を取得
ax.legend(
    handles=[handles[i] for i in order], 
    labels =[labels[i] for i in order], 
    loc=legend_loc, prop={'size': 8}
)

##### 標準化分布の作図 -----

# 標準化分布を描画
ax   = axes[2]
ax2x = ax.twiny()

ax.fill_between(
    x=cr_z_vec, y1=np.zeros_like(cr_z_vec), y2=cr_dens_vec, 
    facecolor='purple', alpha=0.5, 
    label='central region', 
    zorder=9
) # 中央領域
ax.plot(
    z_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='standardized statistic distribution', 
    zorder=10
) # 標準化分布
for idx, p in enumerate([p_pop, p_hat]):
    z = (p - mu_smp) / sigma_smp # 標準化
    ax.axvline(
        x=z, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母比率・標本比率
for idx, z in enumerate([cr_bound_lower, cr_bound_upper]):
    ax.axvline(
        x=z, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='central bounds' if idx == 0 else None, 
        zorder=22
    ) # 中央領域の境界値
ax.hlines(
    y=0.0, xmin=cr_bound_lower, xmax=cr_bound_upper, 
    color='purple', linewidth=2.0, 
    zorder=30
) # 中央領域

ax.set_xlabel('$z = \\frac{\\hat{p} - \\mu_{smp}}{\\sigma_{smp}}$')
ax2x.set_xticks(
    ticks =[z_obs, cr_bound_lower, cr_bound_upper], 
    labels=['$z_{obs}$', '$z_{1-\\frac{\\alpha}{2}}$', '$z_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$N(z \\mid 0, 1)$')
ax.set_title(std_param_lbl, loc='left')
#ax.legend(loc=legend_loc, prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=z_min, xmax=z_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=z_min, xmax=z_max) # (目盛の共通化用)

# ラベルの装飾を調整(表示順の変更用)
order = [1, 2, 0] # 表示順を指定
handles, labels = ax.get_legend_handles_labels() # 凡例情報を取得
ax.legend(
    handles=[handles[i] for i in order], 
    labels =[labels[i] for i in order], 
    loc=legend_loc, prop={'size': 8}
)

plt.show()

二項分布の母比率の信頼区間の計算の図

 標本分布(中図)と標準化統計量の分布(下図)とで、それぞれの標準偏差で横軸のスケールを調整して、信頼区間の境界値(鎖線)と中央領域の境界値(鎖線)を対応させています(そのため形状の違いが分かりにくくなっています)。また、母比率(赤色・破線)と標準化した標本比率(黒色・破線)の位置が対応するように、横軸をスライドしています(そのため母比率と標準化した母比率や、標本比率と標準化した標本比率の位置が対応していません)。

 標準正規分布における  100 (1-\alpha) \% 中央領域(中央  1-\alpha 領域)の境界値  z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}} と、母比率  p_{\mathrm{pop}} 100 (1-\alpha) \% 信頼区間の境界値  L, U が、対応するのを確認できます。
 母比率  p_{\mathrm{pop}} と標本比率  \hat{p} の位置関係や、 \hat{p} の標準偏差  \sigma_{\mathrm{smp}} \approx \sqrt{\frac{\hat{p} (1-\hat{p})}{N}} の大きさ、上側確率が  \frac{\alpha}{2} (下側確率が  1-\frac{\alpha}{2} )となる分位点  z_{\frac{\alpha}{2}} の大きさなどが、母比率  p_{\mathrm{pop}} と信頼区間  [L, U] の位置関係に影響するのが分かります。

 以上で、二項分布の母比率の信頼区間の計算における3つの分布の対応関係を確認しました。

スポンサードリンク

シミュレーション設定の影響

 最後は、二項分布と信頼区間の対応関係について、設定を変更したときの信頼区間の変化をアニメーションで確認します。
 作図コードについては「Toukei-Kentei/code/confidence_interval/estimate_proportion.py at main · anemptyarchive/Toukei-Kentei · GitHub」を参照してください。

サンプルサイズと信頼区間の関係

 サンプルサイズ  N を変化させたときの信頼区間の変化をアニメーションにします。

 サンプルサイズ  N が大きくなるのに応じて、母比率  p_{\mathrm{pop}} の信頼区間の範囲が小さくなっていくのを確認できます。

試行回数と信頼区間の関係

 サンプリングと推定を複数回行った各結果をアニメーションにします。

 試行回数  N・成功確率  p_{\mathrm{pop}}^2 の二項分布  x_i \sim \mathrm{Bin}(N, p_{\mathrm{pop}}) から  I 回(複数回)サンプリングを行い、各試行の標本を用いて(試行ごとに)信頼区間を推定することを考えます。

 左図は、 i 回目(各試行)の推定の様子を表したグラフです。フレームごとに推定結果を切り替えています。
 右図は、 I 回(全試行)の推定の結果(  I 個の信頼区間)を並べたグラフです。フレームごとに推定結果を追加しています。

  i 回目のサンプリングで得られた  1 個の実現値を標本  x_i とします。 i 回目の標本  x_i から求めた標本比率を  \hat{p}_i、確率変数  \hat{p}_i を標準化した確率変数を  z_i、信頼区間を  [L_i, U_i] で表します。
  \hat{p}_i が従う分布(理論分布)は、中心極限定理により、平均  p_{\mathrm{pop}}・分散  \frac{p_{\mathrm{pop}} (1-p_{\mathrm{pop}})}{N} の正規分布  \hat{p}_i \sim \mathcal{N}(p_{\mathrm{pop}}, \frac{p_{\mathrm{pop}} (1-p_{\mathrm{pop}})}{N}) になるため、サンプルサイズ  N が一定であれば試行ごとに変わりません。 \hat{p}_i は確率変数なので、試行ごとに(標本によって)実現値は変わります。
 さらに、 z_i が従う分布(理論分布)も、標準化により、平均  0・分散  1 の正規分布  z_i \sim \mathcal{N}(0, 1) になるため、試行ごとに変わりません。 \hat{p}_i の値によってグラフの表示範囲や表示位置が変わっていますが、 z 軸に対する確率密度(分布の形状)や中央領域は固定です。 z_i も確率変数なので、試行ごとに(標本によって)実現値は変わります。

 左図では、パラメータ  N, p_{\mathrm{pop}}^2 を固定した二項分布から標本(実現値)を100個生成(100回サンプリング)して、有意水準  \alpha を固定して各試行の標本を用いて区間推定を行っています。右図では、100個の信頼区間(100回分の推定結果)を並べています。
  \alpha = 0.05 の場合、95%信頼区間なので、95個程度の区間(95回程度の推定結果)で  p_{\mathrm{pop}} を捉えられている。言い換えると、5個程度の区間(5回程度の推定結果)で  p_{\mathrm{pop}} を捉えられていないのを確認できます。

 母比率  p_{\mathrm{pop}} 100 (1-\alpha) \% 信頼区間とは、同じ母集団から得られた多数の標本に対してそれぞれ区間推定を行ったとき、 100 (1-\alpha) \% の頻度(割合)で区間に真の値  p_{\mathrm{pop}} が含まれるように設計された手順によって得られる区間です。
 同じ手順(設定)で繰り返し試行したとき、真の値が区間に含まれない頻度(割合)をどの程度許容するかは、 \alpha によって決まります。パラメータなどの設定を固定した場合、同じ標本からは同じ区間が求まり(一意に定まり)ます。
 1回の試行において、真の値が区間に含まれる確率が  1-\alpha と解釈することはできません。

 以上で、二項分布の母比率の信頼区間の計算を確認しました。

 この記事では、二項分布の母比率に関する信頼区間を扱いました。次の記事では、2つの正規分布の母平均に関する信頼区間を扱います。

参考文献

おわりに

 Nがサンプルサイズ(二項分布の乱数の数)ではなく二項分布の試行回数パラメータであり、二項分布のサンプル(乱数)は1つのみを扱う、という点がこれまでの内容と異なりますね。母集団をベルヌーイ分布として扱うのであれば、Nはサンプルサイズ(ベルヌーイ分布の乱数の数)と言えるのかもしれませんが。

 この記事で3.4節(1つの母集団に関する信頼区間)の内容は完了です。今回の内容が一番シンプルだったような気もします。ただし、この記事では「中心極限定理により」の一言で済ませている部分に深入りすると、途端に重い内容になる気もします。

 2026年5月6日は、いぎなり東北産の桜ひなのさんの22歳のお誕生日です。

 スタプラ恵比寿方面から九州方面を巡った後にひなもんに引き寄せられて東北方面に辿り着きました。昨年、リリイベで一見産しまして、単独ライブにも行けました。今年は、合同フェスですがまたライブに行く予定です。楽しみ♪
 新参オタクですが、メジャーデビューしたこれからの活動も楽しみです。おめでとうございます!

【次節の内容】

www.anarchive-beta.com