はじめに

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

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

【前節の内容】

いつか書く

【この節の内容】

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

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

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

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

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

区間推定の実装

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

母集団の設定

　母集団のデータが従う分布(母分布・正規分布) $x_n \sim \mathcal{N}(\mu, \sigma^2)$ の母数(パラメータ) $\mu, \sigma^2$ を設定します。 $\sigma^2$ は既知とします。
　正規分布については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

　母分布のパラメータ $\mu, \sigma^2$ を設定します。

# 母平均を指定
mu_pop = 4.0
print(mu_pop)

# 母標準偏差を指定
sigma_pop = 2.0
#sigma_pop = np.sqrt(sigma2_pop)
print(sigma_pop)

# 母分散を計算
#sigma2_pop = 4.0
sigma2_pop = sigma_pop**2
print(sigma2_pop)

4.0
2.0
4.0

　正規分布の平均(母平均・実数) $\mu$ 、標準偏差(母標準偏差・正の値) $\sigma \gt 0$ を指定して、分散(母分散・正の値) $\sigma^2$ を計算します。または、分散 $\sigma^2 \gt 0$ を指定して、標準偏差 $\sigma = \sqrt{\sigma^2}$ を計算します。
　 $\mu$ が区間推定における真の値であり、この値が含まれる区間を求めるのがここでの目的(推測)です。

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

# x軸の範囲を設定
#x_min = -2.0
#x_max = 10.0
k = 3.0
u = 1.0
x_size  = sigma_pop # 基準値を指定
x_size *= k # 定数倍
#x_size  = np.max(np.abs(x_n - mu_pop)) # 基準値を指定
x_size  = np.ceil(x_size /u)*u # u単位で切り上げ
x_min   = mu_pop - x_size
x_max   = mu_pop + x_size
print(x_min, x_max)

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

-2.0 10.0
[-2.    -1.988 -1.976 -1.964 -1.952]

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

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

# 母分布の確率密度を計算
pop_dens_vec = norm.pdf(x=x_vec, loc=mu_pop, scale=sigma_pop)
print(pop_dens_vec[:5])

[0.00221592 0.00225613 0.00229699 0.0023385  0.00238067]

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

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

# 母分布のラベルを作成
pop_param_lbl = f'$\\mu = {mu_pop:.2f}, \\sigma = {sigma_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軸の設定用

ax.plot(
    x_vec, pop_dens_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
ax.axvline(
    x=mu_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population mean', 
    zorder=11
) # 母平均

ax.set_xlabel('$x$')
ax2x.set_xticks(ticks=[mu_pop], labels=['$\\mu$']) # パラメータのラベル
ax.set_ylabel('$N(x \\mid \\mu, \\sigma^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()

　母分布(正規分布)を曲線(黒色)、母平均(平均パラメータ)を垂直線(赤色・破線)で示します。

　平均 $\mu$ ・分散 $\sigma^2$ (標準偏差 $\sigma$ )の正規分布 $\mathcal{N}(x \mid \mu, \sigma^2)$ について、真の値 $\mu$ を区間で求めていきます。

標本の生成

　設定した母分布(正規分布) $x_n \sim \mathcal{N}(\mu, \sigma^2)$ に従う標本(サンプル) $\mathbf{X}$ を作成します。
　正規分布のサンプリングや集計については「【Python】1次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。

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

# サンプルサイズを指定
N = 16

# 標本を生成
x_n = np.random.normal(loc=mu_pop, scale=sigma_pop, size=N)
print(x_n[:5])

[6.66317301 5.43055795 0.90919942 3.9832323  5.24267195]

　サンプルサイズ $N$ を指定して、正規分布に従う乱数(実現値) $\mathbf{X} = \{x_1, \cdots, x_N\}$ を生成します。
　正規分布の乱数生成関数 np.random.normal() の平均の引数 loc に $\mu$ 、標準偏差の引数 scale に $\sigma$ 、サンプルサイズの引数 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_vec, bound_vec = np.histogram(
    a=x_n, bins=class_num, range=(bound_min, bound_max), density=True
)
print(obs_dens_vec[:5])

13
1.0
[-2. -1.  0.  1.  2.]
[-2.5 -1.5 -0.5  0.5  1.5]
[0.     0.     0.     0.0625 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 = np.mean(x_n)
#x_bar = np.sum(x_n) / N
print(x_bar)

# 不偏分散を計算
sigma2_hat = np.var(x_n, ddof=1)
#sigma2_hat = np.sum((x_n - x_bar)**2) / (N-1)
print(sigma2_hat)

# 不偏標準偏差を計算
sigma_hat = np.sqrt(sigma2_hat)
print(sigma_hat)

4.3713903494493325
2.4033731783004253
1.5502816448311658

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

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

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

# 標本のラベルを作成
obs_param_lbl  = f'$N = {N}, '
obs_param_lbl += f'\\bar{{x}} = {x_bar:.2f}, \\hat{{\\sigma}} = {sigma_hat:.2f}$'

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

# 標本を作図
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軸の設定用

ax.bar(
    x=center_vec, height=obs_dens_vec, width=class_size, 
    color='#00A968', alpha=0.5, 
    label='sample', 
    zorder=10
) # 標本
ax.scatter(
    x=x_n, y=np.tile(0.0, reps=N), 
    color='black', alpha=0.33, s=50, clip_on=False, 
    label='sample values', 
    zorder=11
) # 標本
ax.axvline(
    x=x_bar, 
    color='black', linewidth=1.0, linestyle='--', 
    label='sample mean', 
    zorder=12
) # 標本平均

ax.set_xlabel('$x$')
ax2x.set_xticks(ticks=[x_bar], labels=['$\\bar{x}$']) # 統計量のラベル
ax.set_ylabel('$\\frac{N_x}{w N}$')
ax2y.set_ylabel('$\\frac{N_x}{N}$')
dens_vals = ax.get_yticks()          # 密度
freq_vals = dens_vals * class_size*N # 度数
ax2y.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 度数軸目盛
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=0.0, ymax=dens_max)                # (目盛の共通化用)
ax2y.set_ylim(ymin=0.0, ymax=dens_max*class_size*N) # (目盛の共通化用)

plt.show()

　標本(正規分布の乱数) $\mathbf{X}$ の実現値を点(黒色)、密度・度数を棒グラフ(青緑色)、標本平均を垂直線(黒色・破線)で示します。

　サンプルサイズ $N$ が十分に大きいと、標本 $\mathbf{X}$ のヒストグラムの形状が母分布 $\mathcal{N}(\mu, \sigma^2)$ に近付きます。

標本分布の計算

　標本平均 $\bar{x}$ を確率変数とする分布(標本分布・正規分布) $\bar{x} \sim \mathcal{N}(\mu, \frac{\sigma^2}{N})$ のパラメータ $\mu, \frac{\sigma^2}{N}$ を求めます。

　標本分布のパラメータ $\mu, \frac{\sigma^2}{N}$ を計算します。

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

# 標本分布の分散を計算
sigma2_smp = sigma2_pop / N
#sigma2_smp = sigma_smp**2
print(sigma2_smp)

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

4.0
0.25
0.5

　標本 $\mathbf{X}$ を用いて、標本分布の平均 $\mu_{\mathrm{smp}}$ 、分散 $\sigma_{\mathrm{smp}}^2$ 、標準偏差 $\sigma_{\mathrm{smp}}$ を、次の式で計算します。

$\displaystyle \begin{aligned} \mu_{\mathrm{smp}} &= \mu_{\mathrm{pop}} \\ \sigma_{\mathrm{smp}}^2 &= \frac{\sigma_{\mathrm{pop}}^2}{N} \\ \sigma_{\mathrm{smp}} &= \sqrt{\sigma_{\mathrm{smp}}^2} \\ &= \frac{\sigma_{\mathrm{pop}}}{\sqrt{N}} \end{aligned}$

　母分布のパラメータを $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}$ 、標本分布のパラメータを $\mu_{\mathrm{smp}}, \sigma_{\mathrm{smp}}$ で表しています(ここでのみの表記です)。

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

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

[4.29276747e-32 5.72386126e-32 7.62764863e-32 1.01587929e-31
 1.35220762e-31]

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

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

# 標本分布のラベルを作成
smp_param_lbl = f'$\\bar{{x}} = {x_bar:.2f}, '
smp_param_lbl += f'\\mu = {mu_smp:.2f}, \\frac{{\\sigma}}{{\\sqrt{{N}}}} = {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(
    x_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sample distribution', 
    zorder=10
) # 標本分布
ax.scatter(
    x=x_bar, y=0.0, 
    color='#00A968', s=100, clip_on=False, 
    zorder=11
) # 標本平均
ax.axvline(
    x=x_bar, 
    color='black', linewidth=1.0, linestyle='--', 
    label='sample mean', 
    zorder=12
) # 標本平均

ax.set_xlabel('$\\bar{x} = \\frac{1}{n} \\sum_{i=1}^n x_i$')
ax2x.set_xticks(ticks=[x_bar], labels=['$\\bar{x}$']) # 標本のラベル
ax.set_ylabel('$N(\\bar{x} \\mid \\mu, \\frac{\\sigma^2}{n})$')
ax.set_title(smp_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()

　標本分布(正規分布)を曲線(黒色)、標本平均を点(青緑色)と垂直線(黒色・破線)で示します。

　正規分布の再生性により、平均 $\mu$ ・分散 $\frac{\sigma^2}{N}$ (標準偏差 $\frac{\sigma}{\sqrt{N}}$ )の正規分布 $\mathcal{N}(\bar{x} \mid \mu, \frac{\sigma^2}{N})$ になります。
　サンプルサイズ $N$ が大きいほど、標本分布の分散 $\frac{\sigma^2}{N}$ (標準偏差 $\frac{\sigma}{\sqrt{N}}$ )が小さくなります。標本分布の平均 $\mu$ は、 $N$ の影響を受けません。

　サンプルサイズ $N$ が十分に大きいと、標本平均 $\bar{x}$ が母平均 $\mu$ に近付きます。

標準化分布の計算

　標本平均 $\bar{x}$ を標準化した $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}}$ を設定します。

　有意水準 $\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 = (x_min - x_bar) * np.sqrt(N) / sigma_pop
z_max = (x_max - x_bar) * np.sqrt(N) / sigma_pop
print(z_min, z_max)

# z軸の値を作成
z_vec = np.linspace(start=z_min, stop=z_max, num=1001)
print(z_vec[:5])

-12.742780698898665 11.257219301101335
[-12.7427807 -12.7187807 -12.6947807 -12.6707807 -12.6467807]

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

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

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

[2.19217641e-36 2.97556394e-36 4.03657402e-36 5.47275997e-36
 7.41565860e-36]

　「母集団の設定」のときと同様にして、正規分布に従う確率密度 $\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  = f'$\\alpha = {alpha:.2f}, '
std_param_lbl += f'z_{{1-\\frac{{\\alpha}}{{2}}}} = {cr_bound_lower:.2f}, '
std_param_lbl += f'z_{{\\frac{{\\alpha}}{{2}}}} = {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='standard 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=11
) # 外側領域
ax.hlines(
    y=0.0, xmin=cr_bound_lower, xmax=cr_bound_upper, 
    color='purple', linewidth=2.0, 
    zorder=12
) # 中央領域
for z, lbl in zip([cr_bound_lower, cr_bound_upper], ['central bounds', None]):
    ax.axvline(
        x=z, 
        color='black', linewidth=1.0, linestyle='-.', 
        label=lbl, 
        zorder=13
    ) # 中央領域の境界値

ax.set_xlabel('$z = (\\bar{x} - \\mu) \\frac{\\sqrt{n}}{\\sigma}$')
ax2x.set_xticks(
    ticks =[cr_bound_lower, cr_bound_upper], 
    labels=['$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()

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

　標準化により、平均 $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$ が小さいほど、信頼係数 $1-\alpha$ が大きくなります。
　 $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}}) &= \frac{\alpha}{2} + \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}$

です。

信頼区間の計算

　母平均 $\mu$ の信頼区間 $[L, U] = L \leq \mu \leq U$ の下限・上限 $L, U$ を求めます。

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

# 信頼区間の範囲を計算
ci_bound_lower = x_bar - cr_bound_upper * sigma_pop / np.sqrt(N)
ci_bound_upper = x_bar + cr_bound_upper * sigma_pop / np.sqrt(N)
print(ci_bound_lower, ci_bound_upper)

3.3914083571793054 5.35137234171936

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

$\displaystyle \begin{aligned} L &= \bar{x} + z_{1-\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}} \\ &= \bar{x} - z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}} \\ U &= \bar{x} + z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}} \\ &= \bar{x} - z_{1-\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}} \end{aligned}$

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

# 信頼区間のラベルを作成
ci_param_lbl  = f'$\\alpha = {alpha:.2f}, '
ci_param_lbl += f'L = {ci_bound_lower:.2f}, U = {ci_bound_upper:.2f}$'
ci_def_lbl    = '$L = \\bar{x} + z_{1-\\frac{\\alpha}{2}} \\frac{\\sigma}{\\sqrt{N}}$\n'
ci_def_lbl   += '$U = \\bar{x} + z_{\\frac{\\alpha}{2}} \\frac{\\sigma}{\\sqrt{N}}$'

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

ax.axvline(
    x=mu_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population mean', 
    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 x, lbl in zip([ci_bound_lower, ci_bound_upper], ['confidence bounds', None]):
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        label=lbl, 
        zorder=12
    ) # 信頼区間の境界値
    
ax.text(
    x=0.015, y=0.97, 
    s=ci_def_lbl, transform=ax.transAxes, ha='left', va='top', 
    bbox=dict(facecolor='white', alpha=0.8, edgecolor='black', linewidth=0.5), 
    size = 10, 
    zorder=100
) # 計算式のラベル
ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[mu_pop, ci_bound_lower, ci_bound_upper], 
    labels=['$\\mu$', '$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=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

plt.show()

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

　信頼区間は、標本平均 $\bar{x}$ を中心に $\pm z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}}$ の範囲です。 $z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}}$ は、確率変数 $\bar{x}$ の標準偏差 $\frac{\sigma}{\sqrt{N}}$ の $z_{\frac{\alpha}{2}}$ 個分のサイズを表します。
　サンプルサイズ $N$ が大きいほど、 $\frac{\sigma}{\sqrt{N}}$ が小さくなり、信頼区間の範囲(幅) $U-L = 2 z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{N}}$ が小さくなります。

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

スポンサードリンク

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

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

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

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

# ラベルを作成
ci_param_lbl   = f'$N = {N}, '
ci_param_lbl  += f'\\mu = {mu_pop:.2f}, \\sigma = {sigma_pop:.2f}$\n'
ci_param_lbl  += f'$L = {ci_bound_lower:.2f}, U = {ci_bound_upper:.2f}$'
smp_param_lbl  = f'$\\bar{{x}} = {x_bar:.2f}, \\hat{{\\sigma}} = {sigma_hat:.2f}$\n'
smp_param_lbl += f'$\\mu = {mu_smp:.2f}, \\frac{{\\sigma}}{{\\sqrt{{N}}}} = {sigma_smp:.2f}$'
std_param_lbl  = f'$\\alpha = {alpha:.2f}, '
std_param_lbl += f'z_{{1-\\frac{{\\alpha}}{{2}}}} = {cr_bound_lower:.2f}, '
std_param_lbl += f'z_{{\\frac{{\\alpha}}{{2}}}} = {cr_bound_upper:.2f}$'

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

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

ax.bar(
    x=center_vec, height=obs_dens_vec, width=class_size, 
    color='#00A968', alpha=0.5, 
    label='sample', 
    zorder=8
) # 標本
ax.scatter(
    x=x_n, y=np.tile(0.0, reps=N), 
    color='black', alpha=0.33, s=25, clip_on=False, 
    zorder=9
) # 標本
ax.plot(
    x_vec, pop_dens_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
ax.axvline(
    x=mu_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population mean', 
    zorder=20
) # 母平均
ax.axvline(
    x=x_bar, 
    color='black', linewidth=1.0, linestyle='--', 
    zorder=21
) # 標本平均
for x, lbl in zip([ci_bound_lower, ci_bound_upper], ['confidence bounds', None]):
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        label=lbl, 
        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('$x$')
ax2x.set_xticks(
    ticks =[mu_pop, ci_bound_lower, ci_bound_upper], 
    labels=['$\\mu$', '$L$', '$U$']
) # 信頼区間のラベル
ax.set_ylabel('$N(x \\mid \\mu, \\sigma^2)$')
ax.set_title(ci_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) # (目盛の共通化用)
margin_ratio = 0.05 # 余白を設定
ax.set_ylim(ymin=-margin_ratio*dens_max, ymax=(1.0+margin_ratio)*dens_max) # (余白の表示用)

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

ax.scatter(
    x=x_bar, y=0.0, 
    color='#00A968', s=100, clip_on=False, 
    zorder=9
) # 標本平均
ax.plot(
    x_vec, smp_dens_vec, 
    color='black', linewidth=1.0, 
    label='sample distribution', 
    zorder=10
) # 標本分布
ax.axvline(
    x=mu_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=20
) # 母平均
ax.axvline(
    x=x_bar, 
    color='black', linewidth=1.0, linestyle='--', 
    label='sample mean', 
    zorder=21
) # 標本平均
for x in [ci_bound_lower, ci_bound_upper]:
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値

ax.set_xlabel('$\\bar{x} = \\frac{1}{n} \\sum_{i=1}^n x_i$')
ax2x.set_xticks(
    ticks =[ci_bound_lower, x_bar, ci_bound_upper], 
    labels=[
        '$\\bar{x} - z_{\\frac{\\alpha}{2}} \\frac{\\sigma}{\\sqrt{N}}$', 
        '$\\bar{x}$', 
        '$\\bar{x} + z_{\\frac{\\alpha}{2}} \\frac{\\sigma}{\\sqrt{N}}$'
    ]
) # 信頼区間のラベル
ax.set_ylabel('$N(\\bar{x} \\mid \\mu, \\frac{\\sigma^2}{n})$')
ax.set_title(smp_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) # (目盛の共通化用)

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

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

ax.plot(
    z_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='standard 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.axvline(
    x=(mu_pop - x_bar) / sigma_smp, 
    color='red', linewidth=1.0, linestyle='--', 
    zorder=20
) # 母平均
ax.axvline(
    x=0.0, 
    color='black', linewidth=1.0, linestyle='--', 
    zorder=21
) # 標本平均
for z, lbl in zip([cr_bound_lower, cr_bound_upper], ['central bounds', None]):
    ax.axvline(
        x=z, 
        color='black', linewidth=1.0, linestyle='-.', 
        label=lbl, 
        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 = (\\bar{x} - \\mu) \\frac{\\sqrt{n}}{\\sigma}$')
ax2x.set_xticks(
    ticks =[cr_bound_lower, cr_bound_upper], 
    labels=['$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()

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

　以上で、分散が既知の正規分布の母平均の信頼区間の計算における3つの分布の対応関係を確認しました。

スポンサードリンク

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

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

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

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

　サンプルサイズ $N$ が大きくなるのに応じて、母平均 $\mu$ の信頼区間の範囲が小さくなっていくのを確認できます。

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

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

　平均 $\mu$ ・分散 $\sigma^2$ の正規分布 $x_1, \cdots, x_N \sim \mathcal{N}(\mu, \sigma^2)$ から $I$ 回(複数回)サンプリングを行い、それぞれの標本を用いて信頼区間を推定することを考えます。

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

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

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

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

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

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

参考文献

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

東京図書

Amazon

おわりに

　統計検定2級の勉強を始め(始めて半年ほどが過ぎてようやく真面目に取り組み始め)ました。3.4節の内容ですが、この本のシリーズ記事としては、初めて書いた記事です。
　3.4節については他の内容も書くつもりですが、他の節については今のところなのとも言えません。

　私はベイズ統計ネイティブな人でして、頻度統計をまともに勉強したのは初めてと言ってもいいほどなので、新鮮で筆が乗って勢いで書き始めました。
　概ねサクっと理解と共に執筆も進んだのですが、信頼区間と仮説検定とで用語が変わるものいい加減に書いていたのを信頼区間の範囲の表現に書き直したり、95%信頼区間の95%の解釈の説明に苦心したりで、最終段階で悩みました。
　(この記事に限らずですが)間違いがあればぜひ教えてください。よろしくお願いします。

　本でも記事でもこれまでにいくつもの解説がされている分野ですが、このブログらしく「数式・プログラム・図」を対応付ける解説ができたのではないでしょうか。たとえ他と差別化できなかったとしても、私は書かないと理解できない人なので、まぁ書くのですか。
　これまで通り、端的な解説をご希望の方には他を当たってもらうとして、クドクドと回り道をしながらじっくりと理解したい方がいれば読んでみてください。ただし読みやすいとは言っていませんが。

　2026年3月29日は、私立恵比寿中学の元メンバーの柏木ひなたさんの27歳のお誕生日です！

　誕生日当日に大阪でライブがあったというのに観に行けませんでした、、私はなんて弱い(懐具合)。。
　おめでとうございます。次の機会がある頃には私も強くなっていたい。

【次節の内容】

つづく

からっぽのしょこ

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

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