はじめに

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

　この記事では、正規分布の平均パラメータの差の信頼区間をPythonでスクラッチ実装して、計算過程を可視化します。

【前節の内容】

【この節の内容】

はじめに
3.5.1 正規分布の母平均の差の区間推定の実装：母分散が既知の場合
参考文献
おわりに

3.5.1 正規分布の母平均の差の区間推定の実装：母分散が既知の場合

　正規分布(Normal distribution・ガウス分布・Gaussian distribution)の母平均の差(the difference in population means)に関する信頼区間(confidence interval)を実装して、人工的に生成したデータを用いて、パラメータの区間推定を行います。この記事では、母分布の分散パラメータが既知(known variance parameter)の場合を扱います。
　正規分布については「1次元ガウス分布の定義式 - からっぽのしょこ」、母平均の信頼区間については「【Python】3.4.1：正規分布の母平均の信頼区間の実装：母分散が既知の場合【統計検定2級のノート】 - からっぽのしょこ」を参照してください。

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

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

区間推定の実装

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

母集団の設定

　母集団のデータが従う分布(母分布・population distribution・正規分布) $x \sim \mathcal{N}(\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ の母数(パラメータ) $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2$ を設定します。母分散 $\sigma_{\mathrm{pop}}^2$ は既知とします。
　ただし、個々の母集団に注目する際は(表記がゴチャゴチャするのでpopを省略して)、1つ目の母集団のパラメータを $\mu_1, \sigma_1^2$ 、2つ目の母集団のパラメータを $\mu_2, \sigma_2^2$ で表します。
　正規分布のグラフ作成については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

　母分布のパラメータ $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2$ を設定します。

# 母平均を指定
mu_pop_lt = [10.0, 4.0]
print(mu_pop_lt)

# 母標準偏差を指定
sigma_pop_lt = [3.0, 2.0]
#sigma_pop_lt = [np.sqrt(sgm2) for sgm2 in sigma2_pop_lt]
print(sigma_pop_lt)

# 母分散を計算
#sigma2_pop_lt = [9.0, 4.0]
sigma2_pop_lt = [sgm**2 for sgm in sigma_pop_lt]
print(sigma2_pop_lt)

[10.0, 4.0]
[3.0, 2.0]
[9.0, 4.0]

　正規分布の平均(母平均・実数値) $\mu_{\mathrm{pop}}$ 、標準偏差(母標準偏差・正の値) $\sigma_{\mathrm{pop}} \gt 0$ を指定して、分散(母分散・正の値) $\sigma_{\mathrm{pop}}^2$ を計算します。または、分散 $\sigma_{\mathrm{pop}}^2 \gt 0$ を指定して、標準偏差 $\sigma_{\mathrm{pop}} = \sqrt{\sigma_{\mathrm{pop}}^2}$ を計算します。
　2つの母集団の値を指定して、パラメータごとのリストに格納します。

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

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

# 母平均の差を計算
delta_pop = mu_pop_lt[0] - mu_pop_lt[1]
#delta_pop = np.subtract(*mu_pop_lt)
print(delta_pop)

6.0

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

$\displaystyle \delta_{\mathrm{pop}} = \mu_1 - \mu_2$

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

# x軸の範囲を設定
#x_min = -2.0
#x_max = 19.0
k = 3.0
u = 1.0
x_min = np.min([mu - k*sgm for mu, sgm in zip(mu_pop_lt, sigma_pop_lt)]) # 基準値を指定
x_max = np.max([mu + k*sgm for mu, sgm in zip(mu_pop_lt, sigma_pop_lt)]) # 基準値を指定
#x_min = np.hstack(x_lt).min() # 基準値を指定
#x_max = np.hstack(x_lt).max() # 基準値を指定
x_min = np.floor(x_min /u)*u # u単位で切り下げ
x_max = np.ceil(x_max /u)*u  # u単位で切り上げ
print(x_min, x_max)

# x軸の値を作成
x_vec = np.linspace(start=x_min, stop=x_max, num=1001)
print(x_vec[:5])

-2.0 19.0
[-2.    -1.979 -1.958 -1.937 -1.916]

　この例では、指定したパラメータ(または生成したサンプル)を使って、 $x$ 軸の範囲を設定しています。

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

# 母分布の確率密度を計算
pop_dens_lt = [
    norm.pdf(
        x=x_vec, loc=mu_pop_lt[idx], scale=sigma_pop_lt[idx]
    ) for idx in range(2)
]
print(pop_dens_lt[0][:5])

[4.46100753e-05 4.58756849e-05 4.71748889e-05 4.85085095e-05
 4.98773871e-05]

　 $x$ 軸の値ごとに、正規分布に従う確率密度 $\mathcal{N}(x \mid \mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ を計算します。
　正規分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x に $x$ 、平均の引数 loc に $\mu_{\mathrm{pop}}$ 、標準偏差の引数 scale に $\sigma_{\mathrm{pop}}$ を指定します。
　リスト内包表記を使って、母集団ごとに確率密度を計算してリストに格納します。

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

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

# 母分布のラベルを作成
pop_param_lbl  = '$\\mu_1 = '+f'{mu_pop_lt[0]:.2f}, '
pop_param_lbl += '\\sigma_1 = '+f'{sigma_pop_lt[0]:.2f}$\n'
pop_param_lbl += '$\\mu_2 = '+f'{mu_pop_lt[1]:.2f}, '
pop_param_lbl += '\\sigma_2 = '+f'{sigma_pop_lt[1]:.2f}$\n'
pop_param_lbl += '$\\delta_{pop} = \\mu_2 - \\mu_2 = '+f'{delta_pop:.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軸の設定用

for pop_idx in range(2):
    ax.plot(
        x_vec, pop_dens_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    ax.axvline(
        x=mu_pop_lt[pop_idx], 
        color='red', linewidth=1.0, linestyle='--', 
        label='population means' if pop_idx == 0 else None, 
        zorder=11
    ) # 母平均
ax.hlines(
    y=0.0, xmin=mu_pop_lt[0], xmax=mu_pop_lt[1], 
    color='red', linewidth=2.0, 
    label='population mean difference', 
    zorder=12
) # 母平均の差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =mu_pop_lt, 
    labels=['$\\mu_1$', '$\\mu_2$']
) # パラメータのラベル
ax.set_ylabel('$N(x \\mid \\mu_{pop}, \\sigma_{pop}^2)$')
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つの母集団に関して、母分布(正規分布)を曲線、母平均(平均パラメータ)を垂直線(赤色・破線)で示します。また、2つの母平均の差を線分(赤色)で示します。

　平均 $\mu_{\mathrm{pop}}$ ・分散 $\sigma_{\mathrm{pop}}^2$ (標準偏差 $\sigma_{\mathrm{pop}}$ )の正規分布 $\mathcal{N}(x \mid \mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ について、2つの母集団の真の値 $\mu_1, \mu_2$ の差 $\mu_1 - \mu_2$ を区間で求めていきます。

標本の生成

　設定した母分布(正規分布) $x_n \sim \mathcal{N}(\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ に従う標本(サンプル) $\mathbf{X}$ を作成します。
　ただし、1つ目の母集団の標本に関しては $N_1, x_{1,n}, \mathbf{X}_1, \bar{x}_1, \hat{\sigma}_1$ 、2つ目の母集団の標本に関しては $N_2, x_{2,n}, \mathbf{X}_2, \bar{x}_2, \hat{\sigma}_2$ のように添字を使って表します。
　正規分布のサンプリングや集計については「【Python】1次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。

　母分布から標本 $\mathbf{X}$ を生成します。

# サンプルサイズを指定
N_lt = [16, 20]

# 標本を生成
x_lt = [
    np.random.normal(
        loc=mu_pop_lt[idx], scale=sigma_pop_lt[idx], size=N_lt[idx]
    ) for idx in range(2)
]
print(x_lt[0][:5])

[ 7.0392477  12.60066478 12.59538839 14.40491739  8.34312259]

　サンプルサイズ $N$ を指定して、正規分布に従う乱数(実現値) $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成します。
　正規分布の乱数生成関数 np.random.normal() の平均の引数 loc に $\mu_{\mathrm{pop}}$ 、標準偏差の引数 scale に $\sigma_{\mathrm{pop}}$ 、サンプルサイズの引数 size に $N$ を指定します。
　リスト内包表記を使って、母集団ごとに標本を生成してリストに格納します。

　標本 $\mathbf{X}$ を集計します。

# 階級幅を指定
class_size = 1.0

# 階級数を計算
class_num = ((x_max - x_min) // class_size).astype(dtype='int') + 1
print(class_num)

# 階級幅を再設定:(割り切れない場合)
class_size = (x_max - x_min) / (class_num-1)
print(class_size)

# 階級値を作成
center_vec = np.arange(start=x_min, stop=x_max+class_size, step=class_size)
print(center_vec[:5])

# 境界値の範囲を設定
bound_min = x_min - 0.5*class_size
bound_max = x_max + 0.5*class_size

# 境界値を作成
bound_vec = np.arange(start=bound_min, stop=bound_max+0.5*class_size, step=class_size)
print(bound_vec[:5])

# 密度を集計
obs_dens_lt = [
    np.histogram(
        a=x_n, bins=class_num, range=(bound_min, bound_max), density=True
    )[0] for x_n in x_lt
]
print(obs_dens_lt[0][:5])

22
1.0
[-2. -1.  0.  1.  2.]
[-2.5 -1.5 -0.5  0.5  1.5]
[0.     0.     0.     0.     0.0625]

　階級幅 $w$ を指定して、階級数を設定します。
　 $x$ 軸の範囲で階級値、範囲を階級幅 $w$ 1つ分拡げて(上限・下限に $\pm \frac{w}{2}$ して)境界値を作成して、標本 x_n に含まれる要素数(度数 $N_x$ )をカウントして、密度 $\frac{N_x}{w N}$ を求めます。
　リスト内包表記を使って、母集団ごとに標本を集計してリストに格納します。

　標本 $\mathbf{X}$ の標本統計量 $\bar{x}, \hat{\sigma}^2$ を計算します。

# 標本平均を計算
#x_bar_lt = [np.sum(x_n) / N for N, x_n in zip(N_lt, x_lt)]
x_bar_lt = [np.mean(x_n) for x_n in x_lt]
print(x_bar_lt)

# 不偏分散を計算
#sigma2_hat_lt = [
#    np.sum((x_n - x_bar)**2) / (N-1) for N, x_n, x_bar in zip(N_lt, x_lt, x_bar_lt)
#]
sigma2_hat_lt = [np.var(x_n, ddof=1) for x_n in x_lt]
print(sigma2_hat_lt)

# 不偏標準偏差を計算
sigma_hat_lt = [np.sqrt(sgm2) for sgm2 in sigma2_hat_lt]
print(sigma_hat_lt)

[10.391483931726619, 3.3453256988544253]
[13.72375253755655, 3.116352133437969]
[3.704558345816212, 1.7653192723804862]

　標本 $\mathbf{X}$ を用いて、標本平均 $\bar{x}$ 、不偏分散 $\hat{\sigma}^2$ 、不偏標準偏差 $\hat{\sigma}$ を、次の式で計算します。

$\displaystyle \begin{aligned} \bar{x} &= \frac{1}{N} \sum_{n=1}^N x_n \\ \hat{\sigma}^2 &= \frac{1}{N-1} \sum_{n=1}^N (x_n - \bar{x})^2 \\ \hat{\sigma} &= \sqrt{\hat{\sigma}^2} \end{aligned}$

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

　標本平均の差 $d$ を計算します。

# 標本平均の差を計算
d_obs = x_bar_lt[0] - x_bar_lt[1]
#d_obs = np.subtract(*x_bar_lt)
print(d_obs)

7.046158232872193

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

$\displaystyle d = \bar{x}_1 - \bar{x}_2$

　標本 $\mathbf{X}$ のグラフを作成します。

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

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

# 密度軸の範囲を設定
u = 0.05
dens_max = max(np.max(pop_dens_lt), np.max(obs_dens_lt))
dens_max = np.ceil(dens_max /u)*u # u単位で切り上げ

# 度数軸の範囲を設定
freq_max = dens_max * class_size*np.max(N_lt)

# 余白を設定
margin_ratio = 0.05

# 標本のラベルを作成
obs_param_lbl  = f'$N_1 = {N_lt[0]}, '
obs_param_lbl += '\\bar{x}_1 = '+f'{x_bar_lt[0]:.2f}, '
obs_param_lbl += '\\hat{\\sigma}_1 = '+f'{sigma_hat_lt[0]:.2f}$\n'
obs_param_lbl += f'$N_2 = {N_lt[1]}, '
obs_param_lbl += '\\bar{x}_2 = '+f'{x_bar_lt[1]:.2f}, '
obs_param_lbl += '\\hat{\\sigma}_2 = '+f'{sigma_hat_lt[1]:.2f}$\n'
obs_param_lbl += '$d = \\bar{x}_2 - \\bar{x}_2 = '+f'{d_obs:.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軸の設定用
ax2y = ax.twinx() # 第2軸の設定用

for pop_idx in range(2):
    ax.bar(
        x=center_vec, height=obs_dens_lt[pop_idx], 
        width=class_size, 
        color=cmap(pop_idx), alpha=0.5, 
        label='sample' if pop_idx == 0 else None, 
        zorder=10
    ) # 標本
    ax.scatter(
        x=x_lt[pop_idx], y=np.zeros(N_lt[pop_idx]), 
        facecolor='none', edgecolor=cmap(pop_idx), s=50, clip_on=False, 
        label='sample values' if pop_idx == 0 else None, 
        zorder=11
    ) # 標本
    ax.axvline(
        x=x_bar_lt[pop_idx], 
        color='black', linewidth=1.0, linestyle='--', 
        label='sample means' if pop_idx == 0 else None, 
        zorder=12
    ) # 標本平均
ax.hlines(
    y=0.0, xmin=x_bar_lt[0], xmax=x_bar_lt[1], 
    color='black', linewidth=2.0, 
    label='sample mean difference', 
    zorder=13
) # 標本平均の差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =x_bar_lt, 
    labels=['$\\bar{x}_1$', '$\\bar{x}_2$']
) # 統計量のラベル
ax.set_ylabel('$\\frac{N_x}{w N}$')
ax2y.set_ylabel('$\\frac{N_x}{N}$')
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) # (目盛の共通化用)
ax.set_ylim(ymin=-margin_ratio*dens_max, ymax=(1.0+margin_ratio)*dens_max)   # (目盛の共通化用), 余白を追加
dens_vals = ax.get_yticks()                     # 密度を取得
freq_vals = dens_vals * class_size*np.max(N_lt) # 度数を計算
ax2y.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 度数軸目盛
ax2y.set_ylim(ymin=-margin_ratio*freq_max, ymax=(1.0+margin_ratio)*freq_max) # (目盛の共通化用), 余白を追加

plt.show()

　母集団ごとに色分けして、標本(正規分布の乱数)の実現値を点(くり抜き)、密度・度数を棒グラフ(塗りつぶし)、標本平均を垂直線(破線)で示します。また、2つの標本平均の差を線分で示します。

　サンプルサイズ $N_1, N_2$ が十分に大きいと、標本 $\mathbf{X}_1, \mathbf{X}_2$ のヒストグラムの形状がそれぞれの母分布 $\mathcal{N}(\mu_1, \sigma_1^2), \mathcal{N}(\mu_2, \sigma_2^2)$ に近付きます。

標本分布の計算

　標本平均の差 $d$ を確率変数とする分布(標本分布・sampling distribution・正規分布) $d \sim \mathcal{N}(\mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}^2)$ の分散パラメータ $\sigma_{\mathrm{smp}}^2$ を求めます。平均パラメータ $\mu_{\mathrm{smp}}$ は未知です。

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

# 標本分布の平均を計算
mu_smp = mu_pop_lt[0] - mu_pop_lt[1]
#mu_smp = np.subtract(*mu_pop_lt)
print(mu_smp)

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

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

6.0
0.7625
0.873212459828649

　母数 $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2$ とサンプルサイズ $N$ を用いて、標本分布の平均 $\mu_{\mathrm{smp}}$ 、分散 $\sigma_{\mathrm{smp}}^2$ 、標準偏差 $\sigma_{\mathrm{smp}}$ を、次の式で計算します。

$\displaystyle \begin{aligned} \mu_{\mathrm{smp}} &= \mu_1 - \mu_2 \\ &= \delta_{\mathrm{pop}} \\ \sigma_{\mathrm{smp}}^2 &= \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2} \\ \sigma_{\mathrm{smp}} &= \sqrt{\sigma_{\mathrm{smp}}^2} \\ &= \sqrt{ \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2} } \end{aligned}$

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

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

# d軸の範囲を設定
d_min = x_min - mu_pop_lt[1]
d_max = x_max - mu_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])

-6.0 15.0
[-6.    -5.979 -5.958 -5.937 -5.916]

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

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

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

[4.47713015e-42 6.22880732e-42 8.66081765e-42 1.20354328e-41
 1.67152684e-41]

　「母集団の設定」のときと同様にして、正規分布に従う確率密度 $\mathcal{N}(d \mid \mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}^2)$ を計算します。

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

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

# 標本分布のラベルを作成
smp_param_lbl  = '$d_{obs} = '+f'{d_obs:.2f}, '
smp_param_lbl += '\\mu_{smp} = \\mu_1 - \\mu_2 = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} = \\sqrt{\\frac{\\sigma_1}{N_1} + \\frac{\\sigma_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 value', 
    zorder=12
) # 標本平均の差

ax.set_xlabel('$d = \\bar{x}_1 - \\bar{x}_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 right', 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$ の理論分布は、正規分布の性質により、平均 $\mu_1 - \mu_2$ ・分散 $\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}$ (標準偏差 $\sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}}$ )の正規分布 $\mathcal{N}(d \mid \mu_1 - \mu_2, \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2})$ になります。ただし、母平均 $\mu_1, \mu_2$ が未知なので、標本分布の平均は、本来は分からない値です(真の値の可視化にのみ使っています)。
　サンプルサイズ $N_1, N_2$ が大きいほど、標本分布の分散 $\sigma_{\mathrm{smp}}^2 = \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}$ (標準偏差 $\sigma_{\mathrm{smp}} = \sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}}$ )が小さくなります。標本分布の平均 $\mu_{\mathrm{smp}} = \mu_1 - \mu_2$ は、 $N_1, N_2$ の影響を受けません。

　サンプルサイズ $N_1, N_2$ が十分に大きいと、標本平均の差 $d = \bar{x}_1 - \bar{x}_2$ が母平均の差 $\delta_{\mathrm{pop}} = \mu_1 - \mu_2$ に近付きます。

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

　標本平均の差 $d$ を標準化した $z$ を確率変数とする分布(標準化統計量の分布・標準正規分布) $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)

1.198056923142715

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

$\displaystyle \begin{aligned} z &= \frac{ d - \mu_{\mathrm{smp}} }{ \sigma_{\mathrm{smp}} } \\ &= (d - \delta_{\mathrm{pop}}) \sqrt{ \frac{N_1}{\sigma_1^2} + \frac{N_2}{\sigma_2^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])

-14.940416946675551 9.108713094506912
[-14.94041695 -14.91636782 -14.89231869 -14.86826956 -14.84422043]

　グラフ間で対応させるため(真の値 $\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])

[1.34945628e-49 1.93231082e-49 2.76531096e-49 3.95512110e-49
 5.65359139e-49]

　「母集団の設定」のときと同様にして、正規分布に従う確率密度 $\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 right', 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}}$ です。

　標本 $\mathbf{X}_1, \mathbf{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)

5.334693260756412 8.757623204987974

　信頼区間の境界値 $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 Mean 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 right', 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{\sigma_1}{N_1} + \frac{\sigma_2}{N_2}}$ の範囲です。 $z_{\frac{\alpha}{2}} \sqrt{\frac{\sigma_1}{N_1} + \frac{\sigma_2}{N_2}}$ は、確率変数 $d$ の標準偏差 $\sigma_{\mathrm{smp}} = \sqrt{\frac{\sigma_1}{N_1} + \frac{\sigma_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'\\mu_1 = {mu_pop_lt[0]:.2f}, '
pop_param_lbl += f'\\sigma_1 = {sigma_pop_lt[0]:.2f}, '
pop_param_lbl += '\\bar{x}_1 = '+f'{x_bar_lt[0]:.2f}, '
pop_param_lbl += '\\hat{\\sigma}_1 = '+f'{sigma_hat_lt[0]:.2f}$\n'
pop_param_lbl += f'$N_2 = {N_lt[1]}, '
pop_param_lbl += f'\\mu_2 = {mu_pop_lt[1]:.2f}, '
pop_param_lbl += f'\\sigma_2 = {sigma_pop_lt[1]:.2f}, '
pop_param_lbl += '\\bar{x}_2 = '+f'{x_bar_lt[1]:.2f}, '
pop_param_lbl += '\\hat{\\sigma}_2 = '+f'{sigma_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} = \\mu_1 - \\mu_2 = '+f'{mu_smp:.2f}, '
smp_param_lbl += '\\sigma_{smp} = \\sqrt{\\frac{\\sigma_1}{N_1} + \\frac{\\sigma_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 Mean Difference Confidence Interval (known variance)', fontsize=20)

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

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

for pop_idx in range(2):
    ax.bar(
        x=center_vec, height=obs_dens_lt[pop_idx], 
        width=class_size, 
        color=cmap(pop_idx), alpha=0.5, 
        label='sample' if pop_idx == 0 else None, 
        zorder=8
    ) # 標本
    ax.scatter(
        x=x_lt[pop_idx], y=np.zeros(N_lt[pop_idx]), 
        facecolor='none', edgecolor=cmap(pop_idx), s=50, clip_on=False, 
        zorder=9
    ) # 標本
for pop_idx in range(2):
    ax.plot(
        x_vec, pop_dens_lt[pop_idx], 
        color='black', linewidth=1.0, 
        label='population distribution' if pop_idx == 0 else None, 
        zorder=10
    ) # 母分布
    for idx, mu in enumerate([mu_pop_lt[pop_idx], x_bar_lt[pop_idx]]):
        ax.axvline(
            x=mu, 
            color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
            label=['population means', 'sample means'][idx] if pop_idx == 0 else None, 
            zorder=[20, 21][idx]
        ) # 母・標本平均
for delta in [ci_bound_lower, ci_bound_upper]:
    x = delta + mu_pop_lt[1] # 並行移動
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値
ax.hlines(
    y=0.0, xmin=mu_pop_lt[0], xmax=mu_pop_lt[1], 
    color='red', linewidth=2.0, 
    zorder=30
) # 母平均の差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =mu_pop_lt+x_bar_lt+[ci_bound_lower+mu_pop_lt[1], ci_bound_upper+mu_pop_lt[1]], 
    labels=['$\\mu_1$', '$\\mu_2$', '$\\bar{x}_1$', '$\\bar{x}_2$', '$L + \\mu_2$', '$U + \\mu_2$']
) # 信頼区間のラベル
ax.set_ylabel('$N(x \\mid \\mu_{pop}, \\sigma_{pop}^2)$')
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) # (目盛の共通化用)
ax.set_ylim(ymin=-margin_ratio*dens_max, ymax=(1.0+margin_ratio)*dens_max) # 余白を追加

# ラベルの装飾を調整(重なる対策用)
rotations = [0, 0, 0, 0, 45, 45] # 表示角度を指定
labels    = ax2x.get_xticklabels() # 軸情報を取得
for label, r in zip(labels, rotations):
    label.set_rotation(r)

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

# 標本分布を描画
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 mean 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 = \\mu_1 - \\mu_2, d = \\bar{x}_1 - \\bar{x}_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 right', 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 right', 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 right', prop={'size': 8}
)

plt.show()

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

　標準正規分布における $100 (1-\alpha) \%$ 中央領域(中央 $1-\alpha$ 領域)の境界値 $z_{1-\frac{\alpha}{2}}, z_{\frac{\alpha}{2}}$ と、母平均の差 $\delta_{\mathrm{pop}} = \mu_1 - \mu_2$ の $100 (1-\alpha) \%$ 信頼区間の境界値 $L, U$ が、対応するのを確認できます。
　2つの母平均 $\mu_1, \mu_2$ と標本平均 $\bar{x}_1, \bar{x}_2$ の位置関係や、標本平均の差 $d = \bar{x}_1 - \bar{x}_2$ の標準偏差 $\sigma_{\mathrm{smp}} = \sqrt{\frac{\sigma_1}{N_1} + \frac{\sigma_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_mean_difference_known_variance.py at main · anemptyarchive/Toukei-Kentei · GitHub」を参照してください。

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

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

　サンプルサイズ $N_1, N_2$ が大きくなるのに応じて、標本分布の標準偏差パラメータ $\sigma_{\mathrm{smp}}$ が小さくなり、母平均の差 $\delta_{\mathrm{pop}} = \mu_1 - \mu_2$ の信頼区間の範囲が小さくなっていくのを確認できます。

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

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

　2つの母集団に関して、平均 $\mu_{\mathrm{pop}}$ ・分散 $\sigma_{\mathrm{pop}}^2$ の正規分布 $x_1, \cdots, x_N \sim \mathcal{N}(\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ から $I$ 回(複数回)サンプリングを行い、各試行の標本を用いて(試行ごとに)信頼区間を推定することを考えます。

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

　 $i$ 回目のサンプリングで得られた $N_1, N_2$ 個の実現値をまとめて標本 $\mathbf{X}_{i,1} = \{x_{i,1,1}, \cdots, x_{i,1,N_1}\}, \mathbf{X}_{i,2} = \{x_{i,2,1}, \cdots, x_{i,2,N_2}\}$ とします。 $i$ 回目の標本 $\mathbf{X}_{i,1}, \mathbf{X}_{i,2}$ から求めた標本平均を $\bar{x}_{i,1}, \bar{x}_{i,2}$ 、標本平均の差を $d_i$ 、確率変数 $d_i$ を標準化した確率変数を $z_i$ 、信頼区間を $[L_i, U_i]$ で表します。
　 $d_i$ が従う分布(理論分布)は、正規分布の性質により、平均 $\mu_1 - \mu_2$ ・分散 $\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}$ の正規分布 $d_i \sim \mathcal{N}(\mu_1 - \mu_2, \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2})$ になるため、サンプルサイズ $N_1, N_2$ が一定であれば試行ごとに変わりません。 $d_i$ は確率変数なので、試行ごとに(標本によって)実現値は変わります。
　さらに、 $z_i$ が従う分布(理論分布)も、標準化により、平均 $0$ ・分散 $1$ の正規分布 $z_i \sim \mathcal{N}(0, 1)$ になるため、試行ごとに変わりません。 $d_i$ の値によってグラフの表示範囲や表示位置が変わっていますが、 $z$ 軸に対する確率密度(分布の形状)や中央領域は固定です。 $z_i$ も確率変数なので、試行ごとに(標本によって)実現値は変わります。

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

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

　以上で、母分散が既知の正規分布の母平均の差の信頼区間の計算を確認しました。

　この記事では、母分散が既知の場合の正規分布の母平均の差に関する信頼区間を扱いました。次の記事では、母分散が未知の場合を扱います。

参考文献

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

東京図書

Amazon

おわりに

　3.5節(2つの母集団に関する信頼区間)に入りました。
　2つの正規分布(母集団)を扱うので一見複雑になったように感じましたが、区間推定の対象となる母平均の差や標本平均の差が1つの正規分布に従うので、推定の流れは3.4節(1つの母集団に関する信頼区間)と同じであり、それほど難しくなったわけではないですね。パラメータなどの計算式も複雑になりますが、標本分布のパラメータを個別において(定義して)やれば、計算の意図は同じなので同様の式で表現でき、恐くなくなります。

　2026年5月7日は、モーニング娘。の元メンバーの佐藤優樹さんのお誕生日です♪

　おめでとうございます！

　たぶん毎年言っているのですが、この方がいなければ私はアイドルにハマっていませんでした。そして、アイドルにハマっていなければこれほどブログを続けられなかったでしょう。(このブログは技術ブログに偽装したオタクブログなので。)
　というわけで、皆さん1曲聴いてください。よろしくお願いします。

【次節の内容】

からっぽのしょこ

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

【Python】3.5.1：正規分布の母平均の差の区間推定の実装：母分散が既知の場合【統計検定2級のノート】

はじめに

3.5.1 正規分布の母平均の差の区間推定の実装：母分散が既知の場合

区間推定の実装

母集団の設定

標本の生成

標本分布の計算

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

信頼区間の計算

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

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

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

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

参考文献

おわりに