はじめに

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

　この記事では、信頼区間をPythonでスクラッチ実装して、計算過程を可視化します。

【前節の内容】

【この節の内容】

はじめに
3.5.4 二項分布の母比率の差の区間推定の実装
参考文献
おわりに

3.5.4 二項分布の母比率の差の区間推定の実装

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

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

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

区間推定の実装

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

母集団の設定

　母集団のデータが従う分布(母分布・population distribution・二項分布) $x \sim \mathrm{Bin}(N_{\mathrm{pop}}, p_{\mathrm{pop}})$ の母数(パラメータ) $N_{\mathrm{pop}}, p_{\mathrm{pop}}$ を設定します。
　ただし、個々の母集団に注目する際は(表記がゴチャゴチャするのでpopを省略して)、1つ目の母集団のパラメータを $N_1, p_1$ 、2つ目の母集団のパラメータを $N_2, p_2$ で表します。
　二項分布のグラフ作成については「【Python】二項分布の作図 - からっぽのしょこ」を参照してください。

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

# サンプルサイズを指定
N_lt = [20, 15]
print(N_lt)

# 母比率を指定
p_pop_lt = [0.6, 0.25]
print(p_pop_lt)

[20, 15]
[0.6, 0.25]

　二項分布の試行回数(サンプルサイズ・正の整数) $N_{\mathrm{pop}}$ 、成功確率(母比率・0から1の値) $0 \leq p_{\mathrm{pop}} \leq 1$ を指定します。
　2つの母集団の値を指定して、パラメータごとのリストに格納します。

　 $p_1, p_2$ が区間推定における真の値であり、この値の差が含まれる区間を求めるのがここでの目的(推定・推測)です。

　母比率の差 $\delta_{\mathrm{pop}}$ を計算します。

# 母比率の差を計算
delta_pop = p_pop_lt[0] - p_pop_lt[1]
#delta_pop = np.subtract(*p_pop_lt)
print(delta_pop)

0.35

　母比率の差 $\delta_{\mathrm{pop}}$ を、次の式で計算します。

$\displaystyle \delta_{\mathrm{pop}} = p_1 - p_2$

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

# x軸の範囲を設定
x_lower = 0
x_upper = max(N_lt)
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$ 軸の範囲を設定しています。

　また作図用に、正規化した確率変数 $p$ の作図範囲を設定します。

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

-0.025 1.025

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

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

# 母分布の確率を計算
pop_prob_lt = [
    binom.pmf(
        k=x_vec, n=N_lt[idx], p=p_pop_lt[idx]
    ) for idx in range(2)
]
print(pop_prob_lt[0][:5])

[1.09951163e-08 3.29853488e-07 4.70041221e-06 4.23037099e-05
 2.69686150e-04]

　 $x$ 軸の値ごとに、二項分布に従う確率 $\mathrm{Bin}(x \mid N_{\mathrm{pop}}, p_{\mathrm{pop}})$ を計算します。
　二項分布の確率質量関数 sp.stats.binom.pmf() の確率変数の引数 k に $x$ 、試行回数の引数 n に $N_{\mathrm{pop}}$ 、成功確率の引数 p に $p_{\mathrm{pop}}$ を指定します。
　リスト内包表記を使って、母集団ごとに確率密度を計算してリストに格納します。

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

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

# 配色を設定
cmap = plt.get_cmap('tab10') # カラーマップを指定

# 母分布のラベルを作成
pop_param_lbl  = f'$N_1 = {N_lt[0]}, '
pop_param_lbl += f'p_1 = {p_pop_lt[0]:.2f}$\n'
pop_param_lbl += f'$N_2 = {N_lt[1]}, '
pop_param_lbl += f'p_2 = {p_pop_lt[1]:.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軸の設定用

for pop_idx in range(2):
    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color=cmap(pop_idx), s=50, 
        zorder=10
    ) # 母分布
    x = N_lt[pop_idx] * p_pop_lt[pop_idx] # 逆変換
    ax.axvline(
        x=x, 
        color='red', linewidth=1.0, linestyle='--', 
        label='population proportions' if pop_idx == 0 else None, 
        zorder=11
    ) # 母比率

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[N*p for N, p in zip(N_lt, p_pop_lt)], 
    labels=['$N_1 p_1$', '$N_2 p_2$']
) # パラメータのラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(pop_param_lbl, loc='left')
ax.legend(loc='upper right', 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()

　母集団ごとに色分けして、母分布(二項分布)を曲線と点、母比率(成功確率パラメータ)に対応する位置(サンプルサイズ倍の値)を垂直線(赤色・破線)で示します。
　2つの母集団を $k$ で表すとして、横軸は確率変数 $x_k$ に対応します。

　母集団ごとに、母分布のグラフを確認します。

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

# 確率軸の範囲を設定
u = 0.01
prob_max = np.max(pop_prob_lt)
prob_max = np.ceil(prob_max /u)*u # u単位で切り上げ

# 余白を設定
margin_ratio = 0.05

# 母分布を作図
fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)

# 元の分布を描画
for pop_idx in range(2):
    ax   = axes[0, pop_idx]
    ax2x = ax.twiny() # 第2軸の設定用

    # ラベルを作成
    tmp_param_lbl  = f'$N_{pop_idx+1} = {N_lt[pop_idx]}, '
    tmp_param_lbl += f'p_{pop_idx+1} = {p_pop_lt[pop_idx]:.2f}$'

    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color=cmap(pop_idx), s=50, 
        zorder=10
    ) # 母分布
    x = N_lt[pop_idx] * p_pop_lt[pop_idx] # 逆変換
    ax.axvline(
        x=x, 
        color='red', linewidth=1.0, linestyle='--', 
        label='population proportions' if pop_idx == 0 else None, 
        zorder=11
    ) # 母比率

    ax.set_xlabel('$x$')
    ax.set_xticks(x_vec)
    ax2x.set_xticks(
        ticks =[N_lt[pop_idx] * p_pop_lt[pop_idx]], 
        labels=[f'$N_{pop_idx+1} p_{pop_idx+1}$']
    ) # パラメータのラベル
    ax.set_title(tmp_param_lbl, loc='left')
    ax.grid()
    ax.set_xlim(xmin=0.0, xmax=N_lt[pop_idx])   # (目盛の共通化用)
    ax2x.set_xlim(xmin=0.0, xmax=N_lt[pop_idx]) # (目盛の共通化用)
    ax.set_ylim(ymin=-margin_ratio*prob_max, ymax=(1.0+margin_ratio)*prob_max) # (目盛の共通化用)

axes[0, 0].legend(loc='upper left', prop={'size': 8})
axes[0, 0].set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
axes[0, 1].set_yticklabels([])

# 正規化した分布を描画
for pop_idx in range(2):
    ax   = axes[1, pop_idx]
    ax2x = ax.twiny() # 第2軸の設定用
    
    ax.plot(
        x_vec/N_lt[pop_idx], pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec/N_lt[pop_idx], y=pop_prob_lt[pop_idx], 
        color=cmap(pop_idx), s=50, 
        zorder=10
    ) # 母分布
    ax.axvline(
        x=p_pop_lt[pop_idx], 
        color='red', linewidth=1.0, linestyle='--', 
        zorder=11
    ) # 母比率

    ax.set_xlabel('$\\hat{p} = \\frac{x}{n}$')
    ax2x.set_xticks(
        ticks =[p_pop_lt[pop_idx]], 
        labels=[f'$p_{pop_idx}$']
    ) # パラメータのラベル
    ax.grid()
    ax.set_xlim(xmin=0.0, xmax=1.0)   # (目盛の共通化用)
    ax2x.set_xlim(xmin=0.0, xmax=1.0) # (目盛の共通化用)
    ax.set_ylim(ymin=-margin_ratio*prob_max, ymax=(1.0+margin_ratio)*prob_max) # (目盛の共通化用)

axes[1, 0].set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
axes[1, 1].set_yticklabels([])

plt.show()

　上段の図は、横軸をx軸としたグラフです。母分布の図(1つ前の図)を母集団ごとに分割したグラフと言えます。x軸に関して、(母分布の図のときとは異なり、)それぞれ $0$ から $N_k$ の範囲を表示しています。(この図では表示していませんが、)負の値は $p(x_k \lt 0) = 0$ であり、サンプルサイズ(試行回数)より大きい値も $p(x_k \gt N_k) = 0$ であり、確率変数 $x_k$ は $0$ から $N_k$ の整数をとります。
　下段の図は、横軸をp軸としたグラフです。p軸に関して、それぞれ $0$ から $1$ の範囲を表示しています。 $x_k = 0$ のとき $\frac{x_k}{N_k} = \frac{0}{N_k} = 0$ であり、 $x_k = N_k$ のとき $\frac{x_k}{N_k} = \frac{N_k}{N_k} = 1$ であり、正規化した確率変数 $\hat{p}_k = \frac{x_k}{N_k}$ は $0$ から $1$ の連続値をとります。
　正規化については「標本の生成」で扱います。

　横軸を正規化した母分布のグラフを作成します。

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

# 母分布を作図
fig, axes = plt.subplots(
    nrows=2, ncols=1, 
    figsize=(9, 12), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)

# 元分布を描画
ax   = axes[0]
ax2x = ax.twiny() # 第2軸の設定用

# ラベルを作成
tmp_param_lbl  = f'$N_1 = {N_lt[0]}, '
tmp_param_lbl += f'p_1 = {p_pop_lt[0]:.2f}$\n'
tmp_param_lbl += f'$N_2 = {N_lt[1]}, '
tmp_param_lbl += f'p_2 = {p_pop_lt[1]:.2f}$'

for pop_idx in range(2):
    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color=cmap(pop_idx), s=50, 
        zorder=10
    ) # 母分布
    x = N_lt[pop_idx] * p_pop_lt[pop_idx] # 逆変換
    ax.axvline(
        x=x, 
        color='red', linewidth=1.0, linestyle='--', 
        zorder=11
    ) # 母比率

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[N*p for N, p in zip(N_lt, p_pop_lt)], 
    labels=['$N_1 p_1$', '$N_2 p_2$']
) # パラメータのラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(tmp_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

# 正規化した分布を描画
ax   = axes[1]
ax2x = ax.twiny() # 第2軸の設定用

# ラベルを作成
tmp_param_lbl = '$\\delta_{pop} = p_1 - p_2 = '+f'{delta_pop:.2f}$'

for pop_idx in range(2):
    ax.plot(
        x_vec/N_lt[pop_idx], pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec/N_lt[pop_idx], y=pop_prob_lt[pop_idx], 
        color=cmap(pop_idx), s=50, 
        zorder=10
    ) # 母分布
    ax.axvline(
        x=p_pop_lt[pop_idx], 
        color='red', linewidth=1.0, linestyle='--', 
        label='population proportions' if pop_idx == 0 else None, 
        zorder=11
    ) # 母比率
ax.hlines(
    y=0.0, xmin=p_pop_lt[0], xmax=p_pop_lt[1], 
    color='red', linewidth=2.0, 
    label='population proportion difference', 
    zorder=12
) # 母比率の差

ax.set_xlabel('$\\hat{p} = \\frac{x}{p}$')
ax2x.set_xticks(
    ticks =p_pop_lt, 
    labels=['$p_1$', '$p_2$']
) # パラメータのラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(tmp_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)

# ラベルの装飾を調整(軸の対応用)
p_vals = axes[1].get_xticks()   # p軸の値を取得
x_vals  = np.max(N_lt) * p_vals # x軸の値を計算
axes[0].set_xticks(ticks=x_vals)
axes[0].set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

plt.show()

　上図は、横軸がx軸の母分布の図(2つ前の図)です。
　下図は、横軸をx軸からp軸に変換(正規化)した図です。2つの母集団に関して、母比率(成功確率パラメータ)を垂直線(赤色・破線)で示します。また、2つの母比率の差を線分(赤色)で示します。
　2つの母集団を $k$ で表すとして、横軸は正規化した確率変数 $\hat{p}_k = \frac{x_k}{N_k}$ に対応します。

　横軸(確率変数)について、確率変数の最大値(試行回数パラメータ)で割ることで、0から1の値に正規化しています。正規化により、2つの母分布の横方向の位置関係が変わっています。
　2つの母分布を0から1(0から最大値)の範囲に変換することで、2つの母比率の幅(符号付きの幅)が母比率の差に対応します。

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

標本の生成

　設定した母分布(二項分布) $x \sim \mathrm{Bin}(N_{\mathrm{pop}}, p_{\mathrm{pop}})$ に従う標本(サンプル) $x$ を作成します。
　ただし、1つ目の母集団の標本に関しては $x_1, \hat{p}_1$ 、2つ目の母集団の標本に関しては $x_2, \hat{p}_2$ のように添字を使って表します。
　二項分布のサンプリングや集計については「【Python】二項分布の作図 - からっぽのしょこ」を参照してください。

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

# 標本を生成
x_lt = [
    np.random.binomial(
        n=N_lt[idx], p=p_pop_lt[idx], size=1
    )[0] for idx in range(2)
]
print(x_lt)

[13, 4]

　二項分布に従う乱数(実現値) $x$ を生成します。
　二項分布の乱数生成関数 np.random.binomial() の試行回数の引数 n に $N_{\mathrm{pop}}$ 、成功確率の引数 p に $p_{\mathrm{pop}}$ 、サンプルサイズの引数 size に $1$ を指定します。
　リスト内包表記を使って、母集団ごとに標本を生成してリストに格納します。

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

# 標本比率を計算
p_hat_lt = [x / N for N, x in zip(N_lt, x_lt)]
print(p_hat_lt)

[0.65, 0.26666666666666666]

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

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

　2つの母集団の値を計算して、リストに格納します。

　標本比率の差 $d$ を計算します。

# 標本比率の差を計算
d_obs = p_hat_lt[0] - p_hat_lt[1]
#d_obs = np.subtract(*p_hat_lt)
print(d_obs)

0.38333333333333336

　標本比率の差 $d$ を、次の式で計算します。

$\displaystyle d = \hat{p}_1 - \hat{p}_2$

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

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

# 標本のラベルを作成
obs_param_lbl   = f'$x_1 = {x_lt[0]}, '
obs_param_lbl  += '\\hat{p}_1 = \\frac{x_1}{N_1} = '+f'{p_hat_lt[0]:.2f}$\n'
obs_param_lbl  += f'$x_2 = {x_lt[1]}, '
obs_param_lbl  += '\\hat{p}_2 = \\frac{x_2}{N_2} = '+f'{p_hat_lt[1]:.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軸の設定用

for pop_idx in range(2):
    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        label='sample values' if pop_idx == 0 else None, 
        zorder=11
    ) # 標本
    ax.axvline(
        x=x_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        zorder=12
    ) # 標本

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =x_lt, 
    labels=['$x_1$', '$x_2$']
) # 標本のラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(obs_param_lbl, loc='left')
ax.legend(loc='upper right', 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()

　母集団ごとに色分けして、標本(二項分布の乱数)の実現値を点と垂直線(破線)で示します。
　2つの母集団を $k$ で表すとして、横軸は確率変数 $x_k$ に対応します。

　母集団ごとに、標本のグラフを確認します。

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

# 余白を設定
margin_ratio = 0.05

# 標本を作図
fig, axes = plt.subplots(
    nrows=2, ncols=2, 
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)

# 元の分布を描画
for pop_idx in range(2):
    ax   = axes[0, pop_idx]
    ax2x = ax.twiny() # 第2軸の設定用

    # ラベルを作成
    tmp_param_lbl  = f'$x_{pop_idx+1} = {x_lt[pop_idx]}, '
    tmp_param_lbl += f'\\hat{{p}}_{pop_idx+1} = \\frac{{x_{pop_idx+1}}}{{N_{pop_idx+1}}} = {p_hat_lt[pop_idx]:.2f}$'

    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        label='sample value' if pop_idx == 0 else None, 
        zorder=11
    ) # 標本
    ax.axvline(
        x=x_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        zorder=12
    ) # 標本

    ax.set_xlabel('$x$')
    ax.set_xticks(x_vec)
    ax2x.set_xticks(
        ticks =[x_lt[pop_idx]], 
        labels=[f'$x_{pop_idx+1}$']
    ) # 標本のラベル
    ax.set_title(tmp_param_lbl, loc='left')
    ax.grid()
    ax.set_xlim(xmin=0.0, xmax=N_lt[pop_idx])   # (目盛の共通化用)
    ax2x.set_xlim(xmin=0.0, xmax=N_lt[pop_idx]) # (目盛の共通化用)
    ax.set_ylim(ymin=-margin_ratio*prob_max, ymax=(1.0+margin_ratio)*prob_max) # (目盛の共通化用)

axes[0, 0].legend(loc='upper left', prop={'size': 8})
axes[0, 0].set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
axes[0, 1].set_yticklabels([])

# 正規化した分布を描画
for pop_idx in range(2):
    ax   = axes[1, pop_idx]
    ax2x = ax.twiny() # 第2軸の設定用

    ax.plot(
        x_vec/N_lt[pop_idx], pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec/N_lt[pop_idx], y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=p_hat_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        zorder=11
    ) # 標本比率
    ax.axvline(
        x=p_hat_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        zorder=12
    ) # 標本比率ｈ

    ax.set_xlabel('$\\hat{p} = \\frac{x}{n}$')
    ax2x.set_xticks(
        ticks =[p_hat_lt[pop_idx]], 
        labels=[f'$\\hat{{p}}_{pop_idx+1}$']
    ) # 標本のラベル
    ax.grid()
    ax.set_xlim(xmin=0.0, xmax=1.0)   # (目盛の共通化用)
    ax2x.set_xlim(xmin=0.0, xmax=1.0) # (目盛の共通化用)
    ax.set_ylim(ymin=-margin_ratio*prob_max, ymax=(1.0+margin_ratio)*prob_max) # (目盛の共通化用)

axes[1, 0].set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
axes[1, 1].set_yticklabels([])

plt.show()

　上段の図は、 $0$ から $N_k$ の整数値のx軸を横軸として、確率変数 $x_k$ を示したグラフです。標本の図(1つ前の図)を母集団ごとに分割したグラフと言えます。
　下段の図は、 $0$ から $1$ の連続値のp軸を横軸として、正規化した確率変数(標本比率) $\hat{p}_k = \frac{x_k}{N_k}$ を示したグラフです。

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

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

# 標本を作図
fig, axes = plt.subplots(
    nrows=2, ncols=1, 
    figsize=(9, 12), dpi=100, facecolor='white', 
    constrained_layout=True
)
fig.suptitle('Binomial distribution', fontsize=20)

# 標本を描画
ax   = axes[0]
ax2x = ax.twiny() # 第2軸の設定用

# ラベルを作成
tmp_param_lbl  = f'$x_1 = {x_lt[0]}, '
tmp_param_lbl += '\\hat{p}_1 = \\frac{x_1}{N_1} = '+f'{p_hat_lt[0]:.2f}$\n'
tmp_param_lbl += f'$x_2 = {x_lt[1]}, '
tmp_param_lbl += '\\hat{p}_2 = \\frac{x_2}{N_2} = '+f'{p_hat_lt[1]:.2f}$'

for pop_idx in range(2):
    ax.plot(
        x_vec, pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec, y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        label='sample values' if pop_idx == 0 else None, 
        zorder=11
    ) # 標本
    ax.axvline(
        x=x_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        zorder=12
    ) # 標本

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =x_lt, 
    labels=['$x_1$', '$x_2$']
) # 標本のラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(tmp_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

# 標本比率を描画
ax   = axes[1]
ax2x = ax.twiny() # 第2軸の設定用

# ラベルを作成
tmp_param_lbl = '$d = \\hat{p}_2 - \\hat{p}_2 = '+f'{d_obs:.2f}$'

for pop_idx in range(2):
    ax.plot(
        x_vec/N_lt[pop_idx], pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec/N_lt[pop_idx], y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=p_hat_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        zorder=11
    ) # 標本比率
    ax.axvline(
        x=p_hat_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        label='sample proportions' if pop_idx == 0 else None, 
        zorder=12
    ) # 標本比率
ax.hlines(
    y=0.0, xmin=p_hat_lt[0], xmax=p_hat_lt[1], 
    color='black', linewidth=2.0, 
    label='sample proportion difference', 
    zorder=13
) # 標本比率の差

ax.set_xlabel('$\\hat{p} = \\frac{x}{p}$')
ax2x.set_xticks(
    ticks =p_hat_lt, 
    labels=['$\\hat{p}_1$', '$\\hat{p}_2$']
) # 標本のラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(tmp_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)

# ラベルの装飾を調整(軸の対応用)
p_vals = axes[1].get_xticks()   # p軸の値を取得
x_vals  = np.max(N_lt) * p_vals # x軸の値を計算
axes[0].set_xticks(ticks=x_vals)
axes[0].set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

plt.show()

　上図は、横軸がx軸の標本の図(2つ前の図)です。
　下図は、横軸をx軸からp軸に変換(正規化)した図です。2つの母集団に関して、標本比率(正規化した標本)を点と垂直線(破線)で示します。また、2つの標本比率の差を線分で示します。
　2つの母集団を $k$ で表すとして、横軸は正規化した確率変数 $\hat{p}_k = \frac{x_k}{N_k}$ に対応します。

　「母集団の設定」のときと同様に、正規化により、2つの標本の横方向の位置関係が変わっています。
　2つの(母分布と)標本を0から1の範囲に変換することで、2つの標本比率の幅(符号付きの幅)が標本比率の差に対応します。

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

標本分布の計算

　標本比率の差 $d$ を確率変数とする分布(標本分布・sampling distribution・正規分布) $d \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_lt[0] - p_pop_lt[1]
#mu_smp = np.subtract(*p_pop_lt)
print(mu_smp)

# 標本分布の分散を計算
sigma2_smp  = p_hat_lt[0] * (1.0-p_hat_lt[0]) / N_lt[0]
sigma2_smp += p_hat_lt[1] * (1.0-p_hat_lt[1]) / N_lt[1]
#sigma2_smp = np.sum([p * (1.0-p) / N for N, p in zip(N_lt, p_hat_lt)])
print(sigma2_smp)

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

0.35
0.024412037037037038
0.1562435183840822

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

$\displaystyle \begin{aligned} \mu_{\mathrm{smp}} &= p_1 - p_2 \\ &= \delta_{\mathrm{pop}} \\ \sigma_{\mathrm{smp}}^2 &= \frac{p_1 (1-p_1)}{N_1} + \frac{p_2 (1-p_2)}{N_2} \\ &\approx \frac{\hat{p}_1 (1-\hat{p}_1)}{N_1} + \frac{\hat{p}_2 (1-\hat{p}_2)}{N_2} \\ \sigma_{\mathrm{smp}} &= \sqrt{\sigma_{\mathrm{smp}}^2} \\ &= \sqrt{ \frac{\hat{p}_1 (1-\hat{p}_1)}{N_1} + \frac{\hat{p}_2 (1-\hat{p}_2)}{N_2} } \\ &\approx \sqrt{ \frac{\hat{p}_1 (1-\hat{p}_1)}{N_1} + \frac{\hat{p}_2 (1-\hat{p}_2)}{N_2} } \end{aligned}$

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

　標本分布の確率変数 $d$ の作図範囲を設定します。

# d軸の範囲を設定
d_min = p_min - p_pop_lt[1]
d_max = p_max - p_pop_lt[1]
print(d_min, d_max)

# d軸の値を作成
d_vec = np.linspace(start=d_min, stop=d_max, num=1001)
print(d_vec[:5])

-0.275 0.7749999999999999
[-0.275   -0.27395 -0.2729  -0.27185 -0.2708 ]

　グラフ間で対応させるため、 $p$ 軸の範囲を平行移動して、 $d$ 軸の範囲を設定しています。

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

# 標本分布の確率密度を計算
smp_dens_vec = norm.pdf(x=d_vec, loc=mu_smp, scale=sigma_smp)
print(smp_dens_vec[:5])

[0.00085598 0.00087928 0.00090318 0.00092768 0.00095281]

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

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

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

# 標本分布のラベルを作成
smp_param_lbl  = '$d_{obs} = '+f'{d_obs:.2f}, '
smp_param_lbl += '\\mu_{smp} = p_1 - p_2 = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} \\approx \\sqrt{\\frac{\\hat{p}_1 (1-\\hat{p}_1)}{N_1} + \\frac{\\hat{p}_2 (1-\\hat{p}_2)}{N_2}} = '+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(
    d_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sampling distribution', 
    zorder=10
) # 標本分布
ax.scatter(
    x=d_obs, y=0.0, 
    color='#00A968', s=100, 
    zorder=11
) # 標本比率の差
ax.axvline(
    x=d_obs, 
    color='black', linewidth=1.0, linestyle='--', 
    label='sample proportion\nsample value', 
    zorder=12
) # 標本比率の差

ax.set_xlabel('$d = \\hat{p}_1 - \\hat{p}_2$')
ax2x.set_xticks(
    ticks =[d_obs], 
    labels=['$d_{obs}$']
) # 標本のラベル
ax.set_ylabel('$N(d \\mid \\mu_{smp}, \\sigma_{smp}^2)$')
ax.set_title(smp_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=d_min, xmax=d_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=d_min, xmax=d_max) # (目盛の共通化用)

plt.show()

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

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

　サンプルサイズ $N_1, N_2$ が十分に大きいと、標本比率の差 $d = \hat{p}_1 - \hat{p}_2$ が母比率の差 $\delta_{\mathrm{pop}} = p_1 - p_2$ に近付きます。

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

　標本比率の差 $d$ を標準化した $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}}$ を設定します。

　標本比率の差(の標本) $d$ を標準化します。

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

0.21334218326671603

　標本比率の差 $d$ を用いて、標準化した変数 $z$ を、次の式で計算します。

$\displaystyle \begin{aligned} z &= \frac{ d - \mu_{\mathrm{smp}} }{ \sigma_{\mathrm{smp}} } \\ &= (d - \delta_{\mathrm{pop}}) \sqrt{ \frac{N_1}{\hat{p}_1 (1 - \hat{p}_1)} + \frac{N_2}{\hat{p}_2 (1 - \hat{p}_2)} } \end{aligned}$

　母比率の差 $\delta_{\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

　正規分布の分位点関数 sp.stats.norm.ppd() のパーセンタイルの引数 q に $1-\frac{\alpha}{2}$ 、平均の引数 loc に $0$ 、標準偏差の引数 scale に $1$ を指定します。
　パーセンタイル(百分位数)やクオンタイル(分位数)と分位点関数の関係については「母比率の信頼区間の実装」を参照してください。

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

# z軸の範囲を設定
z_min = (d_min - d_obs) / sigma_smp
z_max = (d_max - d_obs) / 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.213508119517636 2.506770653383909
[-4.21350812 -4.20678784 -4.20006756 -4.19334728 -4.186627  ]

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

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

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

[5.56869979e-05 5.72850714e-05 5.89263443e-05 6.06119038e-05
 6.23428623e-05]

　「標本分布の計算」のときと同様にして、正規分布に従う確率密度 $\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{d - \\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='upper left', 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}}$ です。

　標本 $x_1, x_2$ から求めた標本比率の差 $d_{\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}$

　母比率の差 $\delta_{\mathrm{pop}}$ が未知なので、標準化後の値 $z_{\mathrm{obs}}$ は、本来は分からない値です(真の値の可視化にのみ使っています)。しかし、母数に関わらず、 $z_{\mathrm{obs}}$ が従う分布はパラメータが固定(標準正規分布)になるので、 $z_{\mathrm{obs}}$ が含まれる確率は分かります。
　有意水準と中央領域や外側領域の関係については「母平均の信頼区間の実装」を参照してください。

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

信頼区間の計算

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

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

# 信頼区間の範囲を計算
ci_bound_lower = d_obs - cr_bound_upper * sigma_smp
ci_bound_upper = d_obs + cr_bound_upper * sigma_smp
print(ci_bound_lower, ci_bound_upper)

0.07710166448271039 0.6895650021839563

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

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

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

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

# 信頼区間のラベルを作成
ci_param_lbl  = f'$\\alpha = {alpha:.2f}, '
ci_param_lbl += 'L = d + z_{1-\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_lower:.2f}, '
ci_param_lbl += 'U = d + 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 Difference Confidence Interval', fontsize=20)
ax2x = ax.twiny() # 第2軸の設定用

ax.axvline(
    x=delta_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population mean difference', 
    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, delta in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.axvline(
        x=delta, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='confidence bounds' if idx == 0 else None, 
        zorder=12
    ) # 信頼区間の境界値
    
ax.set_xlabel('$\\delta = \\mu_1 - \\mu_2$')
ax2x.set_xticks(
    ticks =[delta_pop, ci_bound_lower, ci_bound_upper], 
    labels=['$\\delta_{pop}$', '$L$', '$U$']
) # 信頼区間のラベル
ax.set_yticks(ticks=[0.0], labels=['']) # (目盛の非表示化用)
ax.set_title(ci_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=d_min, xmax=d_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=d_min, xmax=d_max) # (目盛の共通化用)

plt.show()

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

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

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

スポンサードリンク

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

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

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

# 余白を設定
margin_ratio = 0.05

# ラベルを作成
pop_param_lbl  = f'$N_1 = {N_lt[0]}, '
pop_param_lbl += f'p_1 = {p_pop_lt[0]:.2f}, '
pop_param_lbl += '\\hat{p}_1 = '+f'{p_hat_lt[0]:.2f}$\n'
pop_param_lbl += f'$N_2 = {N_lt[1]}, '
pop_param_lbl += f'p_2 = {p_pop_lt[1]:.2f}, '
pop_param_lbl += '\\hat{p}_2 = '+f'{p_hat_lt[1]:.2f}$'
smp_param_lbl  = '$\\delta_{pop} = '+f'{delta_pop:.2f}, '
smp_param_lbl += f'd = {d_obs:.2f}$\n'
smp_param_lbl += '$\\mu_{smp} = p_1 - p_2 = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} \\approx \\sqrt{\\frac{\\hat{p}_1 (1-\\hat{p}_1)}{N_1} + \\frac{\\hat{p}_2 (1-\\hat{p}_2)}{N_2}} = '+f'{sigma_smp:.2f}$\n'
smp_param_lbl += '$L = d - z_{\\frac{\\alpha}{2}} \\sigma_{smp} = '+f'{ci_bound_lower:.2f}, '
smp_param_lbl += 'U = d + 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 Difference Confidence Interval', fontsize=20)

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

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

for pop_idx in range(2):
    ax.plot(
        x_vec/N_lt[pop_idx], pop_prob_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=x_vec/N_lt[pop_idx], y=pop_prob_lt[pop_idx], 
        color='black', s=50, 
        zorder=10
    ) # 母分布
    ax.scatter(
        x=p_hat_lt[pop_idx], y=0.0, 
        color=cmap(pop_idx), s=100, 
        label='sample' if pop_idx == 0 else None, 
        zorder=15
    ) # 標本比率
    for idx, p, in enumerate([p_pop_lt[pop_idx], p_hat_lt[pop_idx]]):
        ax.axvline(
            x=p, 
            color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
            label=['population proportions', 'sample proportions'][idx] if pop_idx == 0 else None, 
            zorder=[20, 21][idx]
        ) # 母・標本比率
for delta in [ci_bound_lower, ci_bound_upper]:
    x = delta + p_pop_lt[1] # 並行移動
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値
ax.hlines(
    y=0.0, xmin=p_pop_lt[0], xmax=p_pop_lt[1], 
    color='red', linewidth=2.0, 
    zorder=30
) # 母比率の差

ax.set_xlabel('$p, \\hat{p} = \\frac{x}{n}$')
ax2x.set_xticks(
    ticks =p_pop_lt+[p+1e-10 for p in p_hat_lt]+[ci_bound_lower+p_pop_lt[1], ci_bound_upper+p_pop_lt[1]], 
    labels=['$p_1$', '$p_2$', '$\\hat{p}_1$', '$\\hat{p}_2$', '$L + p_2$', '$U + p_2$']
) # パラメータのラベル
ax.set_ylabel('$Bin(x \\mid n_{pop}, p_{pop})$')
ax.set_title(pop_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=p_min, xmax=p_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=p_min, xmax=p_max) # (目盛の共通化用)
'''
# ラベルの装飾を調整(重なる対策用)
rotations = [0, 0, 0, 0, 30, 30] # 表示角度を指定
labels    = ax2x.get_xticklabels() # 軸情報を取得
for label, r in zip(labels, rotations):
    label.set_rotation(r)
'''
# ラベルの装飾を調整(表示順の変更用)
order = [0, 2, 3, 1] # 表示順を指定
handles, labels = ax.get_legend_handles_labels() # 凡例情報を取得
ax.legend(
    handles=[handles[i] for i in order], 
    labels =[labels[i] for i in order], 
    loc='upper left', prop={'size': 8}
)

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

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

ax.scatter(
    x=d_obs, y=0.0, 
    color='#00A968', s=100, clip_on=False, 
    zorder=9
) # 標本比率の差
ax.plot(
    d_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sampling distribution', 
    zorder=10
) # 標本分布
for idx, delta in enumerate([delta_pop, 0.0, d_obs]):
    ax.axvline(
        x=delta, 
        color=['red', 'red', 'black'][idx], linewidth=1.0, linestyle='--', 
        label=['population proportion difference', None, 'saple mean difference'][idx], 
        zorder=[20, 20, 21][idx]
    ) # 母・標本比率の差
for idx, delta in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.axvline(
        x=delta, 
        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, clip_on=False, 
    label='confidence interval', 
    zorder=30
) # 信頼区間

ax.set_xlabel('$\\delta = p_1 - p_2, d = \\hat{p}_1 - \\hat{p}_2$')
ax2x.set_xticks(
    ticks =[delta_pop, d_obs, ci_bound_lower, ci_bound_upper], 
    labels=[
        '$\\delta_{pop}$', 
        '$d_{obs}$', 
        '$d_{obs} - z_{\\frac{\\alpha}{2}} \\sigma_{smp}$', 
        '$d_{obs} + z_{\\frac{\\alpha}{2}} \\sigma_{smp}$'
    ]
) # 信頼区間のラベル
ax.set_ylabel('$N(d \\mid \\mu_{smp}, \\sigma_{smp}^2)$')
ax.set_title(smp_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=d_min, xmax=d_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=d_min, xmax=d_max) # (目盛の共通化用)
'''
# ラベルの装飾を調整(重なる対策用)
rotations   = [0, 0, 30, 30] # 表示角度を指定
halignments = [
    'center', 'center', # (重ならない場合)
    #'left' if delta_pop >= d_obs else 'right', # (重なる場合)
    #'right' if delta_pop >= d_obs else 'left', # (重なる場合)
    'right', 'left'
] # 表示位置を指定
labels      = ax2x.get_xticklabels() # 軸情報を取得
for label, r, ha in zip(labels, rotations, halignments):
    label.set_rotation(r)
    #label.set_ha(ha)
'''
##### 標準化分布の作図 -----

# 標準化分布を描画
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, delta in enumerate([delta_pop, d_obs]):
    z = (delta - 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{d - \\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='upper left', 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='upper left', prop={'size': 8}
)

plt.show()

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

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

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

スポンサードリンク

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

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

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

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

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

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

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

参考文献

改訂版　日本統計学会公式認定　統計検定２級対応　統計学基礎

東京図書

Amazon

おわりに

　2026年5月29日は、私立恵比寿中学の仲村悠菜さんの19歳のお誕生日です。

　まず声が可愛いです。声以外も可愛いですし、歌も上手いですし、魅力はたくさんありますが、可愛い声で喋っているところが私は一番好きです。おめでとうございます。

　そしてたった今、AMEFURASSHIの元メンバーの小島はなさんのBEYOOOOONDSの加入が発表されました！！！！つまり5月29日が加入日ってことになるんですよね？
　めでたいめでたい！！！！！兼オタとして歓喜しています。
　ビヨ曲もさることながら、「Magic of Love」や「イニミニマニモ」を歌う未来があるかもしれない……嗚呼♫

【次節の内容】

つづく

からっぽのしょこ

はてなブログの仕様変更の影響で、数式の表示が崩壊しています。読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】3.5.4：二項分布の母比率の差の区間推定の実装【統計検定2級のノート】

はじめに

3.5.4 二項分布の母比率の差の区間推定の実装

区間推定の実装

母集団の設定

標本の生成

標本分布の計算

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

信頼区間の計算

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

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

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

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

参考文献

おわりに