はじめに

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

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

【前節の内容】

【この節の内容】

はじめに
3.4.2 正規分布の母平均の信頼区間の実装：母分散が未知の場合
参考文献
おわりに

3.4.2 正規分布の母平均の信頼区間の実装：母分散が未知の場合

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

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

# ライブラリを読込
import numpy as np
from scipy.stats import norm, t
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 \gt 0$ を計算します。または、分散 $\sigma^2$ を指定して、標準偏差 $\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}$ を生成します。

# シードを設定:(アニメーションとの対応用)
np.random.seed(10)

# サンプルサイズを指定
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{\hat{\sigma}^2}{N})$ のパラメータ $\mu, \frac{\hat{\sigma}^2}{N}$ を求めます。

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

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

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

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

4.0
0.15021082364377658
0.38757041120779145

　標本 $\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} \\ &\approx \frac{\hat{\sigma}^2}{N} \\ \sigma_{\mathrm{smp}} &= \sqrt{\sigma_{\mathrm{smp}}^2} \\ &= \frac{\sigma_{\mathrm{pop}}}{\sqrt{N}} \\ &\approx \frac{\hat{\sigma}}{\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])

[9.34042108e-53 1.50774209e-52 2.43148349e-52 3.91741203e-52
 6.30537387e-52]

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

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

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

# 標本分布のラベルを作成
smp_param_lbl = f'$\\bar{{x}} = {x_bar:.2f}, '
smp_param_lbl += f'\\mu = {mu_smp:.2f}, \\frac{{\\hat{{\\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{\\hat{\\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})$ になります。ただし、母分散 $\sigma^2$ が未知なので、代わりに標本分散 $\hat{\sigma}^2$ を用いて、平均 $\mu$ ・分散 $\frac{\hat{\sigma}^2}{N}$ (標準偏差 $\frac{\hat{\sigma}}{\sqrt{N}}$ )の正規分布 $\mathcal{N}(\bar{x} \mid \mu, \frac{\hat{\sigma}^2}{N})$ として近似的に計算します。
　サンプルサイズ $N$ が大きいほど、標本分布の分散 $\frac{\hat{\sigma}^2}{N}$ (標準偏差 $\frac{\hat{\sigma}}{\sqrt{N}}$ )が小さくなります。標本分布の平均 $\mu$ は、 $N$ の影響を受けません。

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

標準化分布の計算

　標本平均 $\bar{x}$ を標準化した $t$ を確率変数とする分布(標準化分布・自由度 $N-1$ の標準t分布) $t \sim \mathrm{t}(N-1, 0, 1)$ における中央領域 $[t_{1-\frac{\alpha}{2}}, t_{\frac{\alpha}{2}}] = t_{1-\frac{\alpha}{2}} \leq t \leq t_{\frac{\alpha}{2}}$ の下限・上限 $t_{1-\frac{\alpha}{2}}, t_{\frac{\alpha}{2}}$ を設定します。
　t分布については「【Python】1次元スチューデントのt分布の作図 - からっぽのしょこ」を参照してください。

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

# 信頼係数を指定
gamma = 0.95

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

0.050000000000000044

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

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

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

-2.131449545559323 2.131449545559323

　連続型の確率分布に関して、確率密度関数を $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}$

の関係です。

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

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

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

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

　ただし、標準t分布の場合、 $t_{1-\frac{\alpha}{2}} = -t_{\frac{\alpha}{2}}$ です。

　t分布の分位点関数 sp.stats.t.ppd() のパーセンタイルの引数 q に $1-\frac{\alpha}{2}$ 、自由度の引数 df に $N-1$ 、位置パラメータの引数 loc に $0$ 、尺度パラメータの引数 scale に $1$ を指定します。

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

# t軸の範囲を設定
t_min = (x_min - x_bar) * np.sqrt(N) / sigma_pop
t_max = (x_max - x_bar) * np.sqrt(N) / sigma_pop
print(t_min, t_max)

# t軸の値を作成
t_vec = np.linspace(start=t_min, stop=t_max, num=1001)
print(t_vec[:5])

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

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

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

# 標準化分布の確率密度を計算
std_dens_vec = t.pdf(x=t_vec, df=N-1, loc=0.0, scale=1.0)
print(std_dens_vec[:5])

[1.02612482e-09 1.05484872e-09 1.08442363e-09 1.11487612e-09
 1.14623367e-09]

　 $t$ 軸の値ごとに、t分布に従う確率密度 $\mathrm{t}(t \mid N-1, 0, 1)$ を計算します。
　t分布の確率密度関数 sp.stats.t.pdf() の確率変数の引数 x に $t$ 、自由度の引数 df に $N-1$ 、位置パラメータの引数 loc に $0$ 、尺度パラメータの引数 scale に $1$ を指定します。

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

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

# 中央領域を計算
cr_t_vec    = np.linspace(start=cr_bound_lower, stop=cr_bound_upper, num=501)
cr_dens_vec = t.pdf(x=cr_t_vec, df=N-1, loc=0.0, scale=1.0)

# 外側領域を計算
tail_t_vec    = np.hstack([
    np.linspace(start=t_min, stop=cr_bound_lower, num=251), 
    np.nan, # (塗りつぶしの分割用)
    np.linspace(start=cr_bound_upper, stop=t_max, num=251)
])
tail_dens_vec = t.pdf(x=tail_t_vec, df=N-1, loc=0.0, scale=1.0)

# 標準化分布のラベルを作成
std_param_lbl  = f'$\\alpha = {alpha:.2f}, '
std_param_lbl += f't_{{1-\\frac{{\\alpha}}{{2}}}} = {cr_bound_lower:.2f}, '
std_param_lbl += f't_{{\\frac{{\\alpha}}{{2}}}} = {cr_bound_upper:.2f}$'

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

ax.plot(
    t_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='standard distribution', 
    zorder=10
) # 標準化分布
ax.fill_between(
    x=cr_t_vec, y1=np.zeros_like(cr_t_vec), y2=cr_dens_vec, 
    facecolor='purple', alpha=0.5, 
    label='central region', 
    zorder=11
) # 中央領域
ax.fill_between(
    x=tail_t_vec, y1=np.zeros_like(tail_t_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 t_val, lbl in zip([cr_bound_lower, cr_bound_upper], ['central bounds', None]):
    ax.axvline(
        x=t_val, 
        color='black', linewidth=1.0, linestyle='-.', 
        label=lbl, 
        zorder=13
    ) # 中央領域の境界値

ax.set_xlabel('$t = (\\bar{x} - \\mu) \\frac{\\sqrt{n}}{\\hat{\\sigma}}$')
ax2x.set_xticks(
    ticks =[cr_bound_lower, cr_bound_upper], 
    labels=['$t_{1-\\frac{\\alpha}{2}}$', '$t_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$t(t \\mid n-1, 0, 1)$')
ax.set_title(std_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=t_min, xmax=t_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=t_min, xmax=t_max) # (目盛の共通化用)

plt.show()

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

　標準化により、自由度 $N-1$ ・位置パラメータ $0$ ・尺度パラメータ $1$ のt分布(自由度 $N-1$ の標準t分布) $\mathrm{t}(t \mid N-1, 0, 1)$ になります。
　標準t分布は、中心が0で左右対称な形状なので、 $t_{1-\frac{\alpha}{2}} = -t_{\frac{\alpha}{2}}$ であり、 $|t_{1-\frac{\alpha}{2}}| = t_{\frac{\alpha}{2}}$ です。

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

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

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

$\displaystyle \begin{aligned} P(t \lt t_{1-\frac{\alpha}{2}})+ P(t \gt t_{\frac{\alpha}{2}}) &= \frac{\alpha}{2} + \frac{\alpha}{2} \\ &= \alpha \end{aligned}$

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

$\displaystyle \begin{aligned} P(t \lt t_{1-\frac{\alpha}{2}}) &= \int_{-\infty}^{t_{1-\frac{\alpha}{2}}} \mathrm{t}(t \mid N-1, 0, 1) \mathrm{d} t \\ &= \frac{\alpha}{2} \\ P(t \gt t_{\frac{\alpha}{2}}) &= \int_{t_{\frac{\alpha}{2}}}^{\infty} \mathrm{t}(t \mid N-1, 0, 1) \mathrm{d} t \\ &= \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_hat / np.sqrt(N)
ci_bound_upper = x_bar + cr_bound_upper * sigma_hat / np.sqrt(N)
print(ci_bound_lower, ci_bound_upper)

3.5453035726082454 5.197477126290419

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

$\displaystyle \begin{aligned} L &= \bar{x} + t_{1-\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{N}} \\ &= \bar{x} - t_{\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{N}} \\ U &= \bar{x} + t_{\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{N}} \\ &= \bar{x} - t_{1-\frac{\alpha}{2}} \frac{\hat{\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} + t_{1-\\frac{\\alpha}{2}} \\frac{\\hat{\\sigma}}{\\sqrt{N}}$\n'
ci_def_lbl   += '$U = \\bar{x} + t_{\\frac{\\alpha}{2}} \\frac{\\hat{\\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 t_{\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{N}}$ の範囲です。 $t_{\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{N}}$ は、確率変数 $\bar{x}$ の標準偏差 $\frac{\hat{\sigma}}{\sqrt{N}}$ の $t_{\frac{\alpha}{2}}$ 個分のサイズを表します。
　サンプルサイズ $N$ が大きいほど、 $\frac{\hat{\sigma}}{\sqrt{N}}$ が小さくなり、信頼区間の範囲(幅) $U-L = 2 t_{\frac{\alpha}{2}} \frac{\hat{\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{{\\hat{{\\sigma}}}}{{\\sqrt{{N}}}} = {sigma_smp:.2f}$'
std_param_lbl  = f'$\\alpha = {alpha:.2f}, '
std_param_lbl += f't_{{1-\\frac{{\\alpha}}{{2}}}} = {cr_bound_lower:.2f}, '
std_param_lbl += f't_{{\\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 (unknown 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} - t_{\\frac{\\alpha}{2}} \\frac{\\hat{\\sigma}}{\\sqrt{N}}$', 
        '$\\bar{x}$', 
        '$\\bar{x} + t_{\\frac{\\alpha}{2}} \\frac{\\hat{\\sigma}}{\\sqrt{N}}$'
    ]
) # 信頼区間のラベル
ax.set_ylabel('$N(\\bar{x} \\mid \\mu, \\frac{\\hat{\\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(
    t_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='standard distribution', 
    zorder=10
) # 標準化分布
ax.fill_between(
    x=cr_t_vec, y1=np.zeros_like(cr_t_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 t_val, lbl in zip([cr_bound_lower, cr_bound_upper], ['central bounds', None]):
    ax.axvline(
        x=t_val, 
        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('$t = (\\bar{x} - \\mu) \\frac{\\sqrt{n}}{\\hat{\\sigma}}$')
ax2x.set_xticks(
    ticks =[cr_bound_lower, cr_bound_upper], 
    labels=['$t_{1-\\frac{\\alpha}{2}}$', '$t_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$t(t \\mid n-1, 0, 1)$')
ax.set_title(std_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=t_min, xmax=t_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=t_min, xmax=t_max) # (目盛の共通化用)

plt.show()

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

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

スポンサードリンク

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

　最後は、正規分布と信頼区間の対応関係について、設定を変更したときの信頼区間の変化をアニメーションで確認します。
　作図コードについては「Toukei-Kentei/code/confidence_interval/estimate_mean_unknown_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$ を標準化した確率変数を $t_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$ が一定であれば試行ごとに変わりません。ただし実際には、母分散 $\sigma^2$ が未知のため標本分散 $\hat{\sigma}^2$ で代用して近似分布 $\bar{x}_i \sim \mathcal{N}(\mu, \frac{\hat{\sigma}^2}{N})$ として扱っているので、試行ごとに(標本によって)変わります。 $\bar{x}_i$ は確率変数なので、試行ごとに(標本によって)実現値は変わります。
　さらに、 $t_i$ が従う分布(理論分布)も、標準化により、自由度 $N-1$ ・平均 $0$ ・分散 $1$ のt分布 $t_i \sim \mathrm{t}(N-1, 0, 1)$ になるため、サンプルサイズ $N$ が一定であれば試行ごとに変わりません。標準化によって標本の影響も受けなくなります。サンプルサイズの影響は受けます。 $\bar{x}_i$ の値によってグラフの表示範囲や表示位置が変わっていますが、 $t$ 軸に対する確率密度(分布の形状)や中央領域は固定です。 $t_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

おわりに

　思いの外サクサク進んでおり、この記事の投稿時点で母比率の信頼区間(4記事目)が概ね書き終わり(理解でき)ました。ただ、現時点で既にこの記事(2記事目)も含めてここまでの記事の構成を書き直したい部分が生じております。ここで修正作業に入ると進捗が悪くなるのは経験から分かっているので、とりあえず初稿として投稿しました。修正したとてその先を進めると修正の修正が発生するのは必至で、内容が関連する仮説検定までは進めた方がいいのも分かっているのですが、はてさて。
　他の視点の悩みとしては再構成したい内容というのが、本の数式表記では混乱しそうなのでオリジナルの数式表記に置き換えて解説の手順を整理したいというものです。しかし、本と表現が変わること自体が混乱を招きかねず悩ましいです。いやでもよく考えたら、このブログでオリジナル表記を使うのは日常茶飯事だったかもしれません。
　手元のノートで書き直してみて様子見し(お茶を濁しておき)ます。

　原作を踏襲しつつオリジナルの要素を付け加えて話を発展させる、どこかで聞いた文化だ。

　2026年4月9日は、ukkaの結城りなさんの23歳のお誕生日です。

　サムネのセンターの方です。おめでとうございます！
　私がukkaを追い始めるきっかけになった方です。過去記事を確認したところその頃からもう2年が経ったようです。
　近頃言及したグループやメンバーはこんな話ばかりで誠に悲しい限りなんですが、来月末でグループの解散が決まっておりまして、、、運動・スポーツが得意なアイドルとして活躍の場が広がっていたところだったのに、本当にもう何と言ったらいいのか。

【次節の内容】

つづく

からっぽのしょこ

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

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

はじめに

3.4.2 正規分布の母平均の信頼区間の実装：母分散が未知の場合

区間推定の実装

母集団の設定

標本の生成

標本分布の計算

標準化分布の計算

信頼区間の計算

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

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

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

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

参考文献

おわりに