はじめに

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

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

【前節の内容】

【この節の内容】

はじめに
3.4.3 正規分布の母分散の信頼区間の実装
参考文献
おわりに

3.4.3 正規分布の母分散の信頼区間の実装

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

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

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

区間推定の実装

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

母集団の設定

　母集団のデータが従う分布(母分布・population distribution・正規分布) $x \sim \mathcal{N}(\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ の母数(パラメータ) $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2$ を設定します。母平均 $\mu_{\mathrm{pop}}$ も未知とします。
　正規分布のグラフ作成については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」を参照してください。

　母分布のパラメータ $\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^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_{\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}$ を計算します。

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

　母分布の確率変数 $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_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2)$ を計算します。
　正規分布の確率密度関数 sp.stats.norm.pdf() の確率変数の引数 x に $x$ 、平均の引数 loc に $\mu_{\mathrm{pop}}$ 、標準偏差の引数 scale に $\sigma_{\mathrm{pop}}$ を指定します。

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

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

# 母分布のラベルを作成
pop_param_lbl  = '$\\mu_{pop} = '+f'{mu_pop:.2f}, '
pop_param_lbl += '\\sigma_{pop} = '+f'{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
) # 母分布
for idx, x in enumerate([mu_pop, mu_pop+sigma_pop]):
    ax.axvline(
        x=x, 
        color='red', linewidth=1.0, linestyle='--', 
        label='population parameters' if idx == 0 else None, 
        zorder=11
    ) # 母数
sgm_x    = mu_pop + sigma_pop
sgm_dens = norm.pdf(x=sgm_x, loc=mu_pop, scale=sigma_pop)
ax.hlines(
    y=sgm_dens, xmin=mu_pop, xmax=sgm_x, 
    color='red', linewidth=2.0, 
    label='population sd', 
    zorder=12
) # 母標準偏差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[mu_pop, mu_pop+sigma_pop], 
    labels=['$\\mu_{pop}$', '$\\mu_{pop} + \\sigma_{pop}$']
) # パラメータのラベル
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()

　母分布(正規分布)を曲線、母平均(平均パラメータ)と母平均から母標準偏差(標準偏差パラメータ)1つ分離れた位置を垂直線(赤色・破線)、母標準偏差を線分(赤色)で示します。

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

標本の生成

　設定した母分布(正規分布) $x_n \sim \mathcal{N}(\mu_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^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_{\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_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_{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}$

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

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

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

# 度数軸の範囲を設定
freq_max = dens_max * class_size*N

# 標本のラベルを作成
obs_param_lbl  = f'$N = {N}, '
obs_param_lbl += '\\bar{x} = '+f'{x_bar:.2f}, '
obs_param_lbl += '\\hat{\\sigma} = '+f'{sigma_hat:.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軸の設定用

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.zeros(N), 
    color='black', alpha=0.33, s=50, clip_on=False, 
    label='sample values', 
    zorder=11
) # 標本
for idx, x in enumerate([x_bar, x_bar+sigma_hat]):
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='--', 
        label='sample parameters' if idx == 0 else None, 
        zorder=12
    ) # 標本統計量
sgm_x    = mu_pop + sigma_hat
sgm_dens = norm.pdf(x=sgm_x, loc=mu_pop, scale=sigma_hat) # (表示位置の設定用)
ax.hlines(
    y=sgm_dens, xmin=x_bar, xmax=sgm_x-mu_pop+x_bar, 
    color='black', linewidth=2.0, 
    label='sample sd', 
    zorder=11
) # 不偏標準偏差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[x_bar, x_bar+sigma_hat], 
    labels=['$\\bar{x}$', '$\\bar{x} + \\hat{\\sigma}$']
) # 統計量のラベル
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=0.0, ymax=dens_max)   # (目盛の共通化用)
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]) # 度数軸目盛
ax2y.set_ylim(ymin=0.0, ymax=freq_max) # (目盛の共通化用)

plt.show()

　標本(正規分布の乱数)の実現値を点、密度・度数を棒グラフ(青緑色)、標本平均と標本平均から不偏標準偏差1つ分離れた位置を垂直線(破線)、不偏標準偏差を線分で示します。

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

スケール変換統計量の分布の計算

　不偏分散 $\hat{\sigma}^2$ をスケール変換した $\chi^2$ を確率変数とする分布(スケール変換統計量の分布・distribution of a scale-transformed statistic・自由度 $N-1$ のカイ二乗分布) $\chi^2 \sim \mathrm{chi2}(N-1)$ における中央領域 $[\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}] = \chi^2_{1-\frac{\alpha}{2}} \leq \chi^2 \leq \chi^2_{\frac{\alpha}{2}}$ の下限・上限 $\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}$ を設定します。
　カイ二乗分布のグラフ作成については「いつか書く」を参照してください。

　不偏分散(の標本) $\hat{\sigma}^2$ をスケール変換します。

# 標本統計量をスケール変換
chi2_obs = (N-1) * sigma2_hat / sigma2_pop
#chi2_obs = np.sum((x_n - x_bar)**2) / sigma2_pop
print(chi2_obs)

9.012649418626594

　不偏分散 $\hat{\sigma}^2$ を用いて、スケール変換した変数 $\chi^2$ を、次の式で計算します。

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

　母分散 $\sigma_{\mathrm{pop}}^2$ は未知なので、 $\chi^2_{\mathrm{obs}}$ も未知の値です。そのため推定には使いませんが、作図(真の値の可視化)用に作成しておきます。

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

# 信頼係数を指定
gamma = 0.95

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

0.050000000000000044

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

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

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

6.262137795043253 27.488392863442975

　連続型の確率分布に関して、確率密度関数を $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$ )となる分位点を $\chi^2_p$ で表すことにします(下側確率を基準にする場合は $\chi^2_{1-p}$ と表記する方が自然な表現だと思われます)。

$\displaystyle \chi^2_{p} = F^{-1}(1-p)$

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

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

　中心が0で左右対称な形状の標準正規分布などの場合とは異なり、カイ二乗分布の場合は、 $\chi^2_{1-\frac{\alpha}{2}} \neq -\chi^2_{\frac{\alpha}{2}}$ です。

　カイ二乗分布の分位点関数 sp.stats.chi2.ppd() のパーセンタイルの引数 q に $\frac{\alpha}{2}, 1-\frac{\alpha}{2}$ 、自由度の引数 df に $N-1$ を指定します。

　スケール変換統計量の分布の確率変数 $\chi^2$ の作図範囲を設定します。

# χ2軸の範囲を設定
chi2_min = 0.0
#chi2_max = 40.0
k = 4.0
u = 5.0
chi2_size   = np.sqrt(2.0 * (N-1)) # 基準値を指定
chi2_size  *= k # 定数倍
chi2_center = N-1 # 中心値を指定
chi2_min    = np.max([0.0, chi2_center-chi2_size]) # 正の値
chi2_max    = chi2_center + chi2_size
chi2_min    = np.floor(chi2_min /u)*u # u単位で切り下げ
chi2_max    = np.ceil(chi2_max /u)*u  # u単位で切り上げ
print(chi2_min, chi2_max)

# χ2軸の値を作成
chi2_vec = np.linspace(start=chi2_min, stop=chi2_max, num=1001)
print(chi2_vec[:5])

0.0 40.0
[0.   0.04 0.08 0.12 0.16]

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

　スケール変換統計量の分布の確率密度を計算します。

# 標準化分布の確率密度を計算
std_dens_vec = chi2.pdf(x=chi2_vec, df=N-1)
print(std_dens_vec[:5])

[0.00000000e+00 2.37053446e-15 2.10307788e-13 2.87582629e-12
 1.82885195e-11]

　 $\chi^2$ 軸の値ごとに、カイ二乗分布に従う確率密度 $\mathrm{chi2}(\chi^2 \mid N-1)$ を計算します。
　カイ二乗分布の確率密度関数 sp.stats.chi2.pdf() の確率変数の引数 x に $\chi^2$ 、自由度の引数 df に $N-1$ を指定します。

　スケール変換統計量の分布のグラフを作成します。

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

# 中央領域を計算
cr_chi2_vec = np.linspace(start=cr_bound_lower, stop=cr_bound_upper, num=501)
cr_dens_vec = chi2.pdf(x=cr_chi2_vec, df=N-1)

# 外側領域を計算
tail_chi2_vec    = np.hstack([
    np.linspace(start=chi2_min, stop=cr_bound_lower, num=251), 
    np.nan, # (塗りつぶしの分割用)
    np.linspace(start=cr_bound_upper, stop=chi2_max, num=251)
])
tail_dens_vec = chi2.pdf(x=tail_chi2_vec, df=N-1)

# 標準化分布のラベルを作成
std_param_lbl  = '$\\chi^2_{obs} = '+f'{chi2_obs:.2f}, '
std_param_lbl += '\\nu_{std} = N-1 = '+f'{N-1}$\n'
std_param_lbl += f'$\\alpha = {alpha:.2f}, '
std_param_lbl += '\\chi^2_{1-\\frac{\\alpha}{2}} = '+f'{cr_bound_lower:.2f}, '
std_param_lbl += '\\chi^2_{\\frac{\\alpha}{2}} = '+f'{cr_bound_upper:.2f}$'

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

ax.plot(
    chi2_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='scale-transformed statistic distribution', 
    zorder=10
) # 標準化分布
ax.fill_between(
    x=cr_chi2_vec, y1=np.zeros_like(cr_chi2_vec), y2=cr_dens_vec, 
    facecolor='purple', alpha=0.5, 
    label='central region', 
    zorder=11
) # 中央領域
ax.fill_between(
    x=tail_chi2_vec, y1=np.zeros_like(tail_chi2_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=chi2_obs, 
    color='black', linewidth=1.0, linestyle='--', 
    zorder=14
) # 標準化変数
for idx, chi2_val in enumerate([cr_bound_lower, cr_bound_upper]):
    ax.axvline(
        x=chi2_val, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='central bounds' if idx == 0 else None, 
        zorder=15
    ) # 中央領域の境界値

ax.set_xlabel('$\\chi^2 = \\frac{(n-1) \\hat{\\sigma}^2}{\\sigma_{pop}^2}$')
ax2x.set_xticks(
    ticks =[chi2_obs, cr_bound_lower, cr_bound_upper], 
    labels=['$\\chi^2_{obs}$', '$\\chi^2_{1-\\frac{\\alpha}{2}}$', '$\\chi^2_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$chi2(\\chi^2 \\mid n-1)$')
ax.set_title(std_param_lbl, loc='left')
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=chi2_min, xmax=chi2_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=chi2_min, xmax=chi2_max) # (目盛の共通化用)

plt.show()

　スケール変換統計量の分布(カイ二乗分布)を曲線、スケール変換した不偏分散(の標本)を垂直線(破線)、境界値を垂直線(鎖線)、中央領域(境界値の内側)を線分と塗りつぶし(紫色)、外側の領域(境界値の外側)を塗りつぶし(青色)で示します。

　スケール変換変数 $\chi^2$ の理論分布は、スケール変換により、自由度 $N-1$ のカイ二乗分布 $\mathrm{chi2}(\chi^2 \mid N-1)$ になります。

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

$\displaystyle \begin{aligned} P(\chi^2_{1-\frac{\alpha}{2}} \leq \chi^2 \leq \chi^2_{\frac{\alpha}{2}}) &= \int_{\chi^2_{1-\frac{\alpha}{2}}}^{\chi^2_{\frac{\alpha}{2}}} \mathrm{chi2}(\chi^2 \mid N-1) \mathrm{d} \chi^2 \\ &= 1 - \alpha \end{aligned}$

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

$\displaystyle \begin{aligned} P(\chi^2 \lt \chi^2_{1-\frac{\alpha}{2}})+ P(\chi^2 \gt \chi^2_{\frac{\alpha}{2}}) &= 1 - P(\chi^2_{1-\frac{\alpha}{2}} \leq \chi^2 \leq \chi^2_{\frac{\alpha}{2}}) \\ &= \alpha \end{aligned}$

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

$\displaystyle \begin{aligned} P(\chi^2 \lt \chi^2_{1-\frac{\alpha}{2}}) &= \int_{-\infty}^{\chi^2_{1-\frac{\alpha}{2}}} \mathrm{chi2}(\chi^2 \mid N-1) \mathrm{d} \chi^2 \\ &= \frac{\alpha}{2} \\ P(\chi^2 \gt \chi^2_{\frac{\alpha}{2}}) &= \int_{\chi^2_{\frac{\alpha}{2}}}^{\infty} \mathrm{chi2}(\chi^2 \mid N-1) \mathrm{d} \chi^2 \\ &= \frac{\alpha}{2} \end{aligned}$

です。

　標本 $\mathbf{X}$ から求めた不偏分散 $\hat{\sigma}_{\mathrm{obs}}^2$ をスケール変換した確率変数(標本) $\chi^2_{\mathrm{obs}}$ が中央領域(紫色の範囲) $[\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}]$ に含まれる確率 $P(\chi^2_{1-\frac{\alpha}{2}} \leq \chi^2_{\mathrm{obs}} \leq \chi^2_{\frac{\alpha}{2}})$ が信頼係数 $1-\alpha$ に対応します。言い換えると、外側の領域(青色の範囲) $[-\infty, \chi^2_{1-\frac{\alpha}{2}}], [\chi^2_{\frac{\alpha}{2}}, \infty]$ に含まれる(中央領域に含まれない)確率 $P(\chi^2_{\mathrm{obs}} \lt \chi^2_{1-\frac{\alpha}{2}}) + P(\chi^2_{\mathrm{obs}} \gt \chi^2_{\frac{\alpha}{2}})$ が有意水準 $\alpha$ に対応します。

$\displaystyle \begin{aligned} P(\chi^2_{1-\frac{\alpha}{2}} \leq \chi^2_{\mathrm{obs}} \leq \chi^2_{\frac{\alpha}{2}}) &= 1 - \alpha \\ P(\chi^2_{\mathrm{obs}} \lt \chi^2_{1-\frac{\alpha}{2}}) + P(\chi^2_{\mathrm{obs}} \gt \chi^2_{\frac{\alpha}{2}}) &= \alpha \end{aligned}$

　母分散 $\sigma_{\mathrm{pop}}^2$ が未知なので、スケール変換後の値 $\chi^2_{\mathrm{obs}}$ は、本来は分からない値です(真の値の可視化にのみ使っています)。しかし、母数に関わらず、 $\chi^2_{\mathrm{obs}}$ が従う分布はパラメータが固定(自由度 $N-1$ のカイ二乗分布)になるので、 $\chi^2_{\mathrm{obs}}$ が含まれる確率は分かります。

　不偏分散 $\hat{\sigma}^2$ やスケール変換した変数 $\chi^2$ は確率変数(標本として得られる実現値)なので、それぞれ対応する分布(理論分布)に従って確率的に値が決まるとみなせます。一方、次に求める信頼区間は、標本として得られた実現値(与えられた値) $\hat{\sigma}^2_{\mathrm{obs}}$ や $\chi^2_{\mathrm{obs}}$ に基づいて確定的に値が決まり(一意に定まり)ます。

信頼区間の計算

　母分散 $\sigma_{\mathrm{pop}}^2$ の信頼区間 $[L, U] = L \leq \sigma^2 \leq U$ の下限・上限 $L, U$ を求めます。

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

# 信頼区間の範囲を計算
ci_bound_lower = (N-1) * sigma2_hat / cr_bound_upper
ci_bound_upper = (N-1) * sigma2_hat / cr_bound_lower
print(ci_bound_lower, ci_bound_upper)

1.31148437282597 5.756915426396709

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

$\displaystyle \begin{aligned} L &= \frac{ (N-1) \hat{\sigma}^2 }{ \chi^2_{\frac{\alpha}{2}} } \\ U &= \frac{ (N-1) \hat{\sigma}^2 }{ \chi^2_{1-\frac{\alpha}{2}} } \end{aligned}$

　母分布のパラメータ $\sigma^2$ の作図範囲を設定します。

# σ2軸の範囲を設定
sigma2_min = 0.0
#sigma2_max = 9.0
k = 3.0
u = 1.0
sigma2_size   = sigma2_pop/(N-1) * np.sqrt(2.0*(N-1)) # 基準値を指定
sigma2_size  *= k # 定数倍
sigma2_center = sigma2_pop # 中心値を指定
sigma2_max    = sigma2_center + sigma2_size
sigma2_max    = np.ceil(sigma2_max /u)*u  # u単位で切り上げ
print(sigma2_min, sigma2_max)

0.0 9.0

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

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

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

# 信頼区間のラベルを作成
ci_param_lbl  = f'$\\alpha = {alpha:.2f}, '
ci_param_lbl += 'L = \\frac{(N-1) \\hat{\\sigma}^2}{\\chi^2_{\\frac{\\alpha}{2}}} = '+f'{ci_bound_lower:.2f}, '
ci_param_lbl += 'U = \\frac{(N-1) \\hat{\\sigma}^2}{\\chi^2_{1-\\frac{\\alpha}{2}}} = '+f'{ci_bound_upper:.2f}$'

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

ax.axvline(
    x=sigma2_pop, 
    color='red', linewidth=1.0, linestyle='--', 
    label='population variance', 
    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, sgm2 in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.axvline(
        x=sgm2, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='confidence bounds' if idx == 0 else None, 
        zorder=12
    ) # 信頼区間の境界値
    
ax.set_xlabel('$\\sigma^2$')
ax2x.set_xticks(
    ticks =[sigma2_pop, ci_bound_lower, ci_bound_upper], 
    labels=['$\\sigma_{pop}^2$', '$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=sigma2_min, xmax=sigma2_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=sigma2_min, xmax=sigma2_max) # (目盛の共通化用)

plt.show()

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

　信頼区間は、 $\frac{(N-1) \hat{\sigma}^2}{\chi^2_{\frac{\alpha}{2}}}$ から $\frac{(N-1) \hat{\sigma}^2}{\chi^2_{1-\frac{\alpha}{2}}}$ の範囲です。 $(N-1) \hat{\sigma}^2 = \sum_{n=1}^N (x_n - \bar{x})^2$ は、偏差 $x_n - \bar{x}$ の2乗和を表します。
　中央領域の区間 $[\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}]$ と信頼区間 $[\frac{(N-1) \hat{\sigma}^2}{\chi^2_{\frac{\alpha}{2}}}, \frac{(N-1) \hat{\sigma}^2}{\chi^2_{1-\frac{\alpha}{2}}}]$ とで $\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}$ の位置関係が入れ替わっていることに注意してください。
　不偏分散(標本の分散) $\hat{\sigma}^2$ が大きいほど、 $(N-1) \hat{\sigma}^2$ が大きくなり、信頼区間の範囲(幅) $U-L$ が大きくなります。

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

スポンサードリンク

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

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

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

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

# σ軸の範囲を設定
sigma_min = x_min - mu_pop
sigma_max = x_max - mu_pop

# σ軸の値を作成
sigma_vec = np.linspace(start=sigma_min, stop=sigma_max, num=1001)
sigma_vec = sigma_vec[sigma_vec > 0] # 正の値を抽出

# σ2軸の値を作成
sigma2_vec = np.linspace(start=sigma2_min, stop=sigma2_max, num=1001)
sigma2_vec = sigma2_vec[sigma2_vec > 0] # 正の値を抽出

# 信頼区間の表示位置を指定
margin_ratio = 0.9

# ラベルを作成
pop_param_lbl  = f'$N = {N}$\n'
pop_param_lbl += '$\\mu_{pop} = '+f'{mu_pop:.2f}, '
pop_param_lbl += '\\sigma_{pop} = '+f'{sigma_pop:.2f}$\n'
pop_param_lbl += '$\\bar{x} = '+f'{x_bar:.2f}, '
pop_param_lbl += '\\hat{\\sigma} = '+f'{sigma_hat:.2f}$\n'
ci_param_lbl   = '$\\sigma_{pop}^2 = '+f'{sigma2_pop:.2f}, '
ci_param_lbl  += '\\hat{\\sigma}^2 = '+f'{sigma2_hat:.2f}$\n'
ci_param_lbl  += '$L = \\frac{(N-1) \\hat{\\sigma}^2}{\\chi^2_{\\frac{\\alpha}{2}}} = '+f'{ci_bound_lower:.2f}, '
ci_param_lbl  += 'U = \\frac{(N-1) \\hat{\\sigma}^2}{\\chi^2_{1-\\frac{\\alpha}{2}}} = '+f'{ci_bound_upper:.2f}$'
std_param_lbl  = '$\\chi^2 = '+f'{chi2_obs:.2f}$\n'
std_param_lbl += f'$\\alpha = {alpha:.2f}, '
std_param_lbl += '\\chi^2_{1-\\frac{\\alpha}{2}} = '+f'{cr_bound_lower:.2f}, '
std_param_lbl += '\\chi^2_{\\frac{\\alpha}{2}} = '+f'{cr_bound_upper:.2f}$'

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

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

# 母分布を描画
ax   = axes[0, 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.zeros(N), 
    color='black', alpha=0.33, s=50, clip_on=False, 
    zorder=9
) # 標本
ax.plot(
    x_vec, pop_dens_vec, 
    color='black', linewidth=1.0, 
    label='population distribution', 
    zorder=10
) # 母分布
for idx, x in enumerate([mu_pop, mu_pop+sigma_pop, x_bar, x_bar+sigma_hat]):
    ax.axvline(
        x=x, 
        color=['red', 'red', 'black', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 20, 21, 21][idx]
    ) # 母数・標本統計量
for sgm2 in [ci_bound_lower, ci_bound_upper]:
    x = mu_pop + np.sqrt(sgm2) # 並行移動
    ax.axvline(
        x=x, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値
sgm_x    = mu_pop + sigma_pop
sgm_dens = norm.pdf(x=sgm_x, loc=mu_pop, scale=sigma_pop)
ax.hlines(
    y=sgm_dens, xmin=mu_pop, xmax=sgm_x, 
    color='red', linewidth=1.0, 
    label='population sd', 
    zorder=31
) # 母標準偏差

ax.set_xlabel('$x$')
ax2x.set_xticks(
    ticks =[mu_pop, mu_pop+sigma_pop, x_bar, x_bar+sigma_hat], 
    labels=['$\\mu_{pop}$', '$\\mu_{pop} + \\sigma_{pop}$', '$\\bar{x}$', '$\\bar{x} + \\hat{\\sigma}_{obs}$']
) # パラメータのラベル
ax.set_ylabel('$N(x \\mid \\mu_{pop}, \\sigma_{pop}^2)$')
ax.set_title(pop_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=x_min, xmax=x_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=x_min, xmax=x_max) # (目盛の共通化用)

# ラベルの装飾を調整(重なる対策用)
halignments = [
    'left' if mu_pop >= x_bar else 'right', 
    'left' if mu_pop+sigma_pop >= x_bar+sigma_hat else 'center', 
    'center', 
    'center' if mu_pop+sigma_pop >= x_bar+sigma_hat else 'left'
] # 標準位置を指定
rotations   = [0, 30, 0, 30] # 表示角度を指定
labels      = ax2x.get_xticklabels() # 軸情報を取得
for label, ha, r in zip(labels, halignments, rotations):
    label.set_ha(ha)
    label.set_rotation(r)

##### σ to σ2の作図 -----

# 標準偏差と分散の関係を描画
ax   = axes[1, 0]
ax2x = ax.twiny()
ax2y = ax.twinx()

ax.plot(
    sigma_vec, sigma_vec**2, 
    color='black', linewidth=1.0, 
    zorder=10
) # 変換曲線
for idx, sgm in enumerate([sigma_pop, sigma_hat]):
    ax.vlines(
        x=sgm, ymin=sgm**2, ymax=sigma2_max, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        label=['population parameters', 'sample parameters'][idx], 
        zorder=[20, 21][idx]
    ) # 母・標本標準偏差
    ax.hlines(
        y=sgm**2, xmin=sgm, xmax=sigma_max, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母・不偏分散
for idx, sgm2 in enumerate([ci_bound_lower, ci_bound_upper]):
    ax.vlines(
        x=np.sqrt(sgm2), ymin=sgm2, ymax=sigma2_max, 
        color='black', linewidth=1.0, linestyle='-.', 
        label='confidence bounds' if idx == 0 else None, 
        zorder=22
    ) # 信頼区間の境界値
    ax.hlines(
        y=sgm2, xmin=np.sqrt(sgm2), xmax=sigma_max, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 信頼区間の境界値
tmp_val = margin_ratio * sigma_max # (表示位置の設定用)
ax.vlines(
    x=tmp_val, ymin=ci_bound_lower, ymax=ci_bound_upper, 
    color='purple', linewidth=2.0, 
    label='confidence interval', 
    zorder=30
) # 信頼区間

ax.set_xlabel('$\\sigma$')
ax2x.set_xticks(
    ticks =[sigma_pop, sigma_hat, np.sqrt(ci_bound_lower), np.sqrt(ci_bound_upper)], 
    labels=['$\\sigma_{pop}$', '$\\hat{\\sigma}_{obs}$', '$\\sqrt{L}$', '$\\sqrt{U}$'], 
    rotation=60
) # パラメータのラベル
ax.set_ylabel('$\\sigma^2$')
ax2y.set_yticks(
    ticks =[sigma2_pop, sigma2_hat, ci_bound_lower, ci_bound_upper], 
    labels=['$\\sigma_{pop}^2$', '$\\hat{\\sigma}_{obs}^2$', '$L$', '$U$']
) # パラメータのラベル
ax.set_title(ci_param_lbl, loc='left')
ax.legend(loc='upper left', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=sigma_min, xmax=sigma_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=sigma_min, xmax=sigma_max) # (目盛の共通化用)
ax.set_ylim(ymin=sigma2_min, ymax=sigma2_max)   # (目盛の共通化用)
ax2y.set_ylim(ymin=sigma2_min, ymax=sigma2_max) # (目盛の共通化用)

##### σ2 to χ2の作図 -----

# 分散と標準化変数の関係を描画
ax = axes[1, 1]

ax.plot(
    (N-1)*sigma2_vec/sigma2_pop, sigma2_vec, 
    color='maroon', linewidth=1.0, 
    label='$f(\\hat{\\sigma}^2) = c \\hat{\\sigma}^2$', 
    zorder=10
) # 変換直線
ax.plot(
    (N-1)*sigma2_hat/sigma2_vec, sigma2_vec, 
    color='indigo', linewidth=1.0, 
    label='$f(\\sigma_{pop}^2) = c \\frac{1}{\\sigma_{pop}^2}$', 
    zorder=11
) # 変換曲線
for idx, sgm2 in enumerate([sigma2_pop, sigma2_hat]):
    chi2_val = (N-1) * sgm2 / sigma2_pop # スケール変換
    ax.hlines(
        y=sgm2, xmin=chi2_min, xmax=chi2_val, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母・不偏分散
    ax.vlines(
        x=chi2_val, ymin=sigma2_min, ymax=sgm2, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母数・標準化変数
for idx, chi2_val in enumerate([cr_bound_lower, cr_bound_upper]):
    sgm2 = (N-1) * sigma2_hat / chi2_val # 逆変換
    ax.hlines(
        y=sgm2, xmin=chi2_min, xmax=chi2_val, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 中央領域の境界値
    ax.vlines(
        x=chi2_val, ymin=sigma2_min, ymax=sgm2, 
        color='black', linewidth=1.0, linestyle='-.', 
        zorder=22
    ) # 中央領域の境界値

std_def_lbl = '$\\chi^2 = (n-1) \\hat{\\sigma}^2 \\frac{1}{\\sigma_{pop}^2}, '
std_def_lbl += '\\chi^2 = \\frac{(n-1)}{\\sigma_{pop}^2} \\hat{\\sigma}^2$'
ax.set_xlabel(std_def_lbl)
inv_def_lbl  = '$\\sigma_{pop}^2 = (n-1) \\hat{\\sigma}^2 \\frac{1}{\\chi^2}, '
inv_def_lbl += '\\hat{\\sigma}^2 = \\frac{\\sigma_{pop}^2}{(n-1)} \\chi^2$'
ax.set_ylabel(inv_def_lbl)
ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=chi2_min, xmax=chi2_max)
ax.set_ylim(ymin=sigma2_min, ymax=sigma2_max)

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

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

ax.fill_between(
    x=cr_chi2_vec, y1=np.zeros_like(cr_chi2_vec), y2=cr_dens_vec, 
    facecolor='purple', alpha=0.5, 
    label='central region', 
    zorder=9
) # 中央領域
ax.plot(
    chi2_vec, std_dens_vec, 
    color='black', linewidth=1.0, 
    label='scale-transformed\nstatistic distribution', 
    zorder=10
) # 標準化分布
for idx, sgm2 in enumerate([sigma2_pop, sigma2_hat]):
    chi2_val = (N-1) * sgm2 / sigma2_pop # スケール変換
    ax.axvline(
        x=chi2_val, 
        color=['red', 'black'][idx], linewidth=1.0, linestyle='--', 
        zorder=[20, 21][idx]
    ) # 母・不偏分散
for idx, chi2_val in enumerate([cr_bound_lower, cr_bound_upper]):
    ax.axvline(
        x=chi2_val, 
        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('$\\chi^2 = \\frac{(n-1) \\hat{\\sigma}^2}{\\sigma_{pop}^2}$')
ax2x.set_xticks(
    ticks =[chi2_obs, cr_bound_lower, cr_bound_upper], 
    labels=['$\\chi^2_{obs}$', '$\\chi^2_{1-\\frac{\\alpha}{2}}$', '$\\chi^2_{\\frac{\\alpha}{2}}$']
) # 中央領域のラベル
ax.set_ylabel('$chi2(\\chi^2 \\mid n-1)$')
ax.set_title(std_param_lbl, loc='left')
#ax.legend(loc='upper right', prop={'size': 8})
ax.grid()
ax.set_xlim(xmin=chi2_min, xmax=chi2_max)   # (目盛の共通化用)
ax2x.set_xlim(xmin=chi2_min, xmax=chi2_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}
)

##### その他の設定 -----

# 不要なサブプロットを非表示
axes[0, 1].axis('off')
axes[2, 0].axis('off')

plt.show()

　中段左図は、標準偏差 $\sigma_{\mathrm{pop}}, \hat{\sigma}$ と分散 $\sigma_{\mathrm{pop}}^2, \hat{\sigma}^2$ の関係を表す曲線です。母分散 $\sigma_{\mathrm{pop}}^2$ の信頼区間 $[L, U]$ を(縦方向の)線分(紫色)で示します。
　中段右図は、分散 $\sigma_{\mathrm{pop}}, \hat{\sigma}$ と分散をスケール変換した変数 $\chi^2$ の関係を表す直線と曲線です。(この2色は、凡例との対応のみで、他の図との対応関係はないです。)
　不偏分散 $\hat{\sigma}^2$ から $\chi^2$ への変換では、母分散(真の値) $\sigma_{\mathrm{pop}}^2$ を固定して(定数として)計算します。 $\hat{\sigma}^2$ の関数 $\chi^2 = \frac{n-1}{\sigma_{\mathrm{pop}}^2} \hat{\sigma}^2$ なので、直線(線形な変換)になります。
　分位点 $\chi^2_{\alpha}$ から $L, U$ への変換では、不偏分散(観測値) $\hat{\sigma}^2$ を固定して(定数として)計算します。 $\chi^2$ の関数 $\sigma_{\mathrm{pop}}^2 = (n-1) \hat{\sigma}^2 \frac{1}{\chi^2}$ なので、曲線(非線形な変換)になります。

　カイ二乗分布における $100 (1-\alpha) \%$ 中央領域(中央 $1-\alpha$ 領域)の境界値 $\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}$ と、母分散 $\sigma_{\mathrm{pop}}^2$ の $100 (1-\alpha) \%$ 信頼区間の境界値 $L, U$ が、対応するのを確認できます。中央領域の区間 $[\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}]$ と信頼区間 $[\frac{(N-1) \hat{\sigma}^2}{\chi^2_{\frac{\alpha}{2}}}, \frac{(N-1) \hat{\sigma}^2}{\chi^2_{1-\frac{\alpha}{2}}}]$ とで $\chi^2_{1-\frac{\alpha}{2}}, \chi^2_{\frac{\alpha}{2}}$ の位置関係が入れ替わるのも図から分かります。
　母数 $\sigma_{\mathrm{pop}}, \sigma_{\mathrm{pop}}^2$ と標本統計量 $\hat{\sigma}, \hat{\sigma}^2$ の位置関係や、 $\bar{x}$ の不偏分散 $\hat{\sigma}^2$ の大きさ、上側確率が $1-\frac{\alpha}{2}$ (下側確率が $\frac{\alpha}{2}$ )となる分位点 $\chi^2_{1-\frac{\alpha}{2}}$ と上側確率が $\frac{\alpha}{2}$ (下側確率が $1-\frac{\alpha}{2}$ )となる分位点 $\chi^2_{\frac{\alpha}{2}}$ の大きさなどが、母分散 $\sigma_{\mathrm{pop}}^2$ と信頼区間 $[L, U]$ の位置関係に影響するのが分かります。

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

スポンサードリンク

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

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

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

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

　サンプルサイズ $N$ が大きくなり不偏分散 $\hat{\sigma}^2$ (不偏標準偏差 $\hat{\sigma}$ )が小さくなるのに応じて、母分散 $\sigma_{\mathrm{pop}}^2$ の信頼区間の範囲が小さくなっていくのを確認できます。 11

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

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

　平均 $\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$ 個の実現値をまとめて標本 $\mathbf{X}_i = \{x_{i,1}, \cdots, x_{i,N}\}$ とします。 $i$ 回目の標本 $\mathbf{X}_i$ から求めた不偏分散を $\hat{\sigma}_i^2$ 、確率変数 $\hat{\sigma}_i^2$ を標準化(スケール)した確率変数を $\chi^2_i$ 、信頼区間を $[L_i, U_i]$ で表します。
　 $\chi^2_i$ が従う分布(理論分布)は、自由度 $N-1$ のカイ二乗分布 $\chi^2_i \sim \mathrm{chi2}(N-1)$ になるため、サンプルサイズ $N$ が一定であれば試行ごとに変わりません。標準化によって標本の影響も受けなくなります。サンプルサイズの影響は受けます。 $\chi^2_i$ も確率変数なので、試行ごとに(標本によって)実現値は変わります。

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

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

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

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

参考文献

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

東京図書

Amazon

おわりに

　前回の記事からひと月ほど空いての更新になってしまいました。
　前記事にも書きましたが、本から逸脱する表記を使って再構成しようかと悩みながら、(既に書き終わっていたこの記事を投稿せずに、)3.5節(2つの母集団に関する信頼区間)の解説をいくつか書いていたところ、独自の表記(母分布や標本分布ごとにパラメータをおくなど)を用いないと扱いきれなくなり、本とのシームレスな対応は諦めて独自路線を突き進むことにしました。そんなこんなで、3.4節(1つの母集団に関する信頼区間)の解説から(つまりこれまでに書いていた全ての内容を)全面的に書き直していました。そして、書き直し終わりました。
　分布ごとにパラメータをおいたことで個々の計算(操作の意図)を分かりやすく整理できたと思うのですが、全体を通して眺めると、似たような記号(パラメータ)が増えたことで判別しにくくなったり、母集団との各分布や統計量などとの関係性が分かりにくくなったりしたような気がします。その上、両方の表記を図に載せたことで情報過多になってしまったと思われます。
　解説の読み方や図の見方が伝われば理解しやすくなったとは思っているのですが、これが私の表現力の限界です。ちなみに、3.5節に進むと図に載せる情報がさらに増えていきます。よければ頑張って読んでください。

　といった苦労とは別に、分散に関するアレコレを図で表現するのが悩ましかったです。標準偏差であれば分布上に描き示せるのですが。ちなみに、3.5.3項の母分散の比の描き表し方をイメージできそうになかったので保留にして、そこまでの記事を投稿することにした次第です。先に、F分布の解説記事を書けば何か思い付きそうな気もしますが、どうでしょうか。その前に、カイ二乗分布の解説記事を書く必要がありそうです。

　図を考える中で、母分散と不偏分散のどちらが固定して扱う(計算する)のかに悩んだのですが、見方(操作)によってどちらを固定するかが変わるんですね。この違いが、標本統計量を含む区間の式やその確率と、真の値を含む区間の式やそれが確率ではない理由、に対応するということですね。実感してよく分かりました。

　2026年5月5日は、えびちゅうこと私立恵比寿中学のメジャーデビュー14周年とココユノノカの加入5周年の日です！！！

　エビ中がここまで活動を続けていなければ、ゆのぴが加入していなければ、私がスタプラにハマることはなかったです。そしてこのブログを続けるモチベーションも維持できなかったことでしょう。
　おめでとうございます。ありがとうございます。これからもよろしくお願いします。

【次節の内容】

からっぽのしょこ

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

【Python】3.4.3：正規分布の母分散の信頼区間の実装【統計検定2級のノート】

はじめに

3.4.3 正規分布の母分散の信頼区間の実装

区間推定の実装

母集団の設定

標本の生成

スケール変換統計量の分布の計算

信頼区間の計算

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

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

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

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

参考文献

おわりに