はじめに
前回までは、確率分布の母平均の点推定について焦点を当てていました。
母平均の最尤推定量である標本平均は一般には誤差を含み、母平均に一致しません。そこで母平均を高い確率で含む区間を求めることにします。
そのようなパラメータを含む区間を求める技術を区間推定と言います。
この記事では区間推定の仕方と、$t$ 分布による信頼区間の計算方法を述べます。
また、この記事で使用しているコードは以下のGoogle Colabで試すことができます。
\begin{align*}
\newcommand{\mat}[1]{\begin{pmatrix} #1 \end{pmatrix}}
\newcommand{\f}[2]{\frac{#1}{#2}}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\d}[2]{\frac{{\rm d}#1}{{\rm d}#2}}
\newcommand{\e}{{\rm e}}
\newcommand{\T}{\mathsf{T}}
\newcommand{\(}{\left(}
\newcommand{\)}{\right)}
\newcommand{\{}{\left\{}
\newcommand{\}}{\right\}}
\newcommand{\[}{\left[}
\newcommand{\]}{\right]}
\newcommand{\dis}{\displaystyle}
\newcommand{\eq}[1]{{\rm Eq}(\ref{#1})}
\newcommand{\n}{\notag\\}
\newcommand{\t}{\ \ \ \ }
\newcommand{\tt}{\t\t\t\t}
\newcommand{\argmax}{\mathop{\rm arg\, max}\limits}
\newcommand{\argmin}{\mathop{\rm arg\, min}\limits}
\def\l<#1>{\left\langle #1 \right\rangle}
\def\us#1_#2{\underset{#2}{#1}}
\def\os#1^#2{\overset{#2}{#1}}
\newcommand{\case}[1]{\{ \begin{array}{ll} #1 \end{array} \right.}
\newcommand{\s}[1]{{\scriptstyle #1}}
\definecolor{myblack}{rgb}{0.27,0.27,0.27}
\definecolor{myred}{rgb}{0.78,0.24,0.18}
\definecolor{myblue}{rgb}{0.0,0.443,0.737}
\definecolor{myyellow}{rgb}{1.0,0.82,0.165}
\definecolor{mygreen}{rgb}{0.24,0.47,0.44}
\newcommand{\c}[2]{\textcolor{#1}{#2}}
\newcommand{\ub}[2]{\underbrace{#1}_{#2}}
\end{align*}
大標本に対する区間推定
今回の目的は、母平均を高い確率で含む区間を求めることです。まずは大標本の場合で考えてみます。
サンプルサイズ $n$ が大きい標本は大標本と呼ばれます。大標本に対して、標本平均を $\bar{x}_n = \f{1}{n} \sum_{i=1}^n x^{(i)}$ とすると、
\begin{align*}
\f{\bar{x}_n – \mu}{\sigma / \sqrt{n}}
\end{align*}
は中心極限定理から、標準正規分布に従います（ここで、$\mu, \sigma^2$ は母平均、母分散です）。
すなわち、
\begin{align*}
P\(a < \f{\bar{x}_n – \mu}{\sigma / \sqrt{n}} < b \) \simeq \int_a^b \f{1}{\sqrt{2 \pi}} \e^{-\f{x^2}{2}} {\rm d}x
\end{align*}
が成り立ちます。求めたいのは母平均 $\mu$ を含む区間ですから、左辺の事象を $\mu$ を挟む形に変形します。
\begin{align*}
P\(\bar{x}_n – b \f{\sigma}{\sqrt{n}} < \mu < \bar{x}_n – a \f{\sigma}{\sqrt{n}} \) \simeq \int_a^b \f{1}{\sqrt{2 \pi}} \e^{-\f{x^2}{2}} {\rm d}x.
\end{align*}
さて、ここで $a = -1.96, b=1.96$ とすると、右辺の積分値は $0.95$ となることがわかっています。
\begin{align*}
P\(\bar{x}_n – 1.96 \f{\sigma}{\sqrt{n}} < \mu < \bar{x}_n + 1.96 \f{\sigma}{\sqrt{n}} \) \simeq \int_{-1.96}^{1.96} \f{1}{\sqrt{2 \pi}} \e^{-\f{x^2}{2}} {\rm d}x = 0.95.
\end{align*}
この区間
\begin{align*}
\bar{x}_n – 1.96 \f{\sigma}{\sqrt{n}} < \mu < \bar{x}_n + 1.96 \f{\sigma}{\sqrt{n}}
\end{align*}
は $95\%$ 信頼区間呼ばれ、「母平均 $\mu$ がこの区間にある」と言えば $95 \%$ の確率で正解であることを意味します。
このように、真の推定値が含まれる区間を確率的に推定することを区間推定と言います。
ただし、上記の信頼区間 $\bar{x}_n \pm 1.96 \f{\sigma}{\sqrt{n}} $ を計算するには、サンプルサイズ $n$ が大きい標本であり、母集団の標準偏差 $\sigma$ が既知である必要があります。そうではない場合に使用されるのが次の $t$ 分布の理論です。
$t$ 分布の理論
サンプルサイズが小さい場合は、中心極限定理による正規分布への近似では誤差が大きくなり実用上問題となります。また、一般には母集団の標準偏差 $\sigma$ は未知であるため、信頼区間を計算できません。
この問題を解決するのが、$t$ 分布の理論です。
先ほどまでは大標本である場合、標本平均を $\bar{X}$ として
\begin{align*}
\f{\bar{X} – \mu}{\sigma / \sqrt{n}}
\end{align*}
が標準正規分布に従うことを用いていました。
一般には $\sigma$ は未知であるので、不偏分散: $U^2 = \f{1}{n-1} \sum_{i=1}^n \(X_i – \bar{X} \)^2$ の平方根で置き換えます。
\begin{align*}
\f{\bar{X} – \mu}{U / \sqrt{n}}
\end{align*}
こうして定義される確率変数は、自由度 $n-1$ の $t$ 分布に従うことがわかっています（なぜ $t$ 分布に従うかはこの記事の最後に述べます）。
$\bar{X}$ を確率変数 $X_i \us{ \sim }_{\rm i.i.d.} \mathcal{N}(\mu, \sigma^2)$ の標本平均 $\bar{X} = \f{1}{n} \sum_{i=1}^n$ とするとき、
\begin{align*}
T = \f{\bar{X} – \mu}{U / \sqrt{n}}, \t U^2 = \f{1}{n-1} \sum_{i=1}^n \(X_i – \bar{X} \)^2
\end{align*}
で定義される確率変数 $T$ は自由度 $n-1$ の $t$ 分布: $t(n-1)$に従う。
自由度 $n$ の $t$ 分布の確率関数は次式で表されます。
\begin{align*}
f(x; n) = \f{1}{\sqrt{n}\, B\(\f{1}{2}, \f{n}{2}\)} \( 1 + \f{x^2}{n} \)^{-\f{n+1}{2}}, \t (-\infty < x < \infty).
\end{align*}
ここで、$B(x, y)$ はベータ関数で、ガンマ関数: $\Gamma(x)$ と次の関係があります。
\begin{align*}
B(x, y) = \f{\Gamma(x) \Gamma(y)}{\Gamma(x + y)}, \t \Gamma(x) = \int_0^{\infty} u^{x-1} \e^{-u} {\rm d}u.
\end{align*}
$t$ 分布の理論を用いることにより、信頼区間を計算することができます。
$t_{0.025}(n-1)$ を確率変数 $T$ がその値以上となる確率が $2.5\%$ となる点とすると（下図参照）、
\begin{align*}
P\(-t_{0.025}(n-1) < \f{\bar{x}_n – \mu}{U / \sqrt{n}} < t_{0.025}(n-1) \) = 0.95
\end{align*}
が成り立ちます。
よって、$t$ 分布の理論による $95\%$ 信頼区間は
\begin{align*}
P\(\bar{x}_n – t_{0.025}(n-1) \f{\sigma}{\sqrt{n}} < \mu < \bar{x}_n + t_{0.025}(n-1) \f{\sigma}{\sqrt{n}} \) = 0.95
\end{align*}
より、次のようになります。
\begin{align*}
\bar{x}_n – t_{0.025}(n-1) \f{\sigma}{\sqrt{n}} < \mu < \bar{x}_n + t_{0.025}(n-1) \f{\sigma}{\sqrt{n}} .
\end{align*}
この信頼区間の計算はライブラリを用いて簡単に行うことができます。
さて、$\f{\bar{X} – \mu}{U / \sqrt{n}}$ がなぜ $t$ 分布に従うかの証明は本記事の最後で述べるとして、本当に従っているのか、シミュレーションで確認してみます。
以下では、サンプルサイズ $n = 5$ として $\mathcal{N}(50, 10)$ の正規乱数を生成し、$\f{\bar{X} – \mu}{U / \sqrt{n}}$ を計算する過程を $20000$ 回おこなって、ヒストグラムを作成しています。
import numpy as np
import pandas as pd
import seaborn as sns
import scipy as sp
import scipy.stats as st
import matplotlib.pyplot as plt
population_mean = 50
population_var = 10
norm_dist = st.norm(loc=population_mean, scale=population_var)
sample_size = 5
t_samples = []
for _ in range(200000):
samples = norm_dist.rvs(size=sample_size)
sample_mean = st.tmean(samples)
unbiased_var = st.tvar(samples)
t = (sample_mean - population_mean) / np.sqrt(unbiased_var / sample_size)
t_samples.append(t)
# plot
sns.displot(t_samples, stat='density', alpha=0.6)
x = np.linspace(-6, 6, 300)
t_dist = st.t.pdf(x, sample_size-1)
norm_dist = st.norm.pdf(x, loc=0, scale=1)
plt.plot(x, t_dist, color=colors[0], label='t dist')
plt.plot(x, norm_dist, color=colors[2], linestyle='--', label='standard normal dist')
plt.xlim(-6, 6)
plt.legend(fontsize=8)
plt.show()
上図は赤線が自由度 $4$ の $t$ 分布、黄色の破線が標準正規分布を表しています。ヒストグラムをみれば、自由度 $4$ の $t$ 分布によくフィットしていることがわかります。
pythonによる信頼区間の計算
pythonで信頼区間を計算するには、scipy.stats ライブラリの t.interval 関数で行うことができます。
今回は実験として、考える母集団が平均 $50$、分散 $10$ の正規分布 $\mathcal{N}(50, 10)$ に従うとして、この母集団から $100$ 個サンプリングすることにします。
import numpy as np
import pandas as pd
import seaborn as sns
import scipy as sp
import scipy.stats as st
import matplotlib.pyplot as plt
population_mean = 50
population_var = 10
norm_dist = st.norm(loc=population_mean, scale=population_var)
sample_size = 10
alpha = 0.95
deg_of_freedom = sample_size - 1
samples = norm_dist.rvs(size=sample_size)
sample_mean = st.tmean(samples) # 標本平均
sample_var = st.tvar(samples) # 不偏分散
t_interval = st.t.interval(alpha, deg_of_freedom, loc=sample_mean, scale=np.sqrt(sample_var/sample_size))
print(f'95%信頼区間: {t_interval}')# 出力
95%信頼区間: (44.49084252434853, 58.87622334835179)信頼区間を計算すると、95%信頼区間: (44.49084252434853, 58.87622334835179) となりました。注意点として、信頼区間は確率変数ですので、実行の度に95%信頼区間の値は変動します。
次に、母平均と $95\%$ 信頼区間の関係性をシミュレーションで確認してみます。
シミュレーションの手順は、母集団が平均 $50$、分散 $10$ の正規分布 $\mathcal{N}(50, 10)$ に従うとして、下記を実行します。
- 母集団から $100$ 個サンプリングし$95\%$ 信頼区間を計算する。
- 1の試行を $10000$ 回繰り返す。
- $10000$ パターンの信頼区間の内、母平均がその区間に収まっているものをカウントする。
population_mean = 50
population_var = 10
norm_dist = st.norm(loc=population_mean, scale=population_var)
alpha = 0.95
sample_size = 10
deg_of_freedom = sample_size - 1
results = []
for _ in range(10000):
samples = norm_dist.rvs(size=sample_size)
sample_mean = st.tmean(samples)
sample_var = st.tvar(samples)
min_val_t, max_val_t = st.t.interval(alpha, deg_of_freedom, loc=sample_mean, scale=np.sqrt(sample_var/sample_size))
correct = min_val_t < population_mean < max_val_t
results.append({'min_val_t': min_val_t, 'max_val_t': max_val_t, 'correct': correct})
df_result = pd.DataFrame(results)
print(df_result)
result_series = df_result['correct'].value_counts()
print(result_series)
print(f'信頼区間が母平均を含む確率: {result_series[True]/result_series.sum()}')# 出力
min_val_t max_val_t correct
0 39.574942 57.104679 True
1 49.214656 57.570718 True
2 38.208271 59.257566 True
3 47.996852 60.674018 True
4 32.060261 48.765342 False
... ... ... ...
9995 44.120682 59.360833 True
9996 38.679976 55.922974 True
9997 41.938874 61.338453 True
9998 42.027188 52.692872 True
9999 40.332269 47.691603 False
True 9496
False 504
Name: correct, dtype: int64
信頼区間が母平均を含む確率: 0.9496先ほど述べた通り、信頼区間は確率変数ですので、実行の度に $95\%$ 信頼区間の値は変動します。そのため、$95\%$ 信頼区間は母平均を含んだり、含まなかったりします。
今回のシミュレーションでは、$95\%$ 信頼区間が母平均を含む確率は $94.96\%$ となりました。
この様に、無限回の試行を行えば、$95\%$ 信頼区間が母平均を含む確率は $95\%$ となります。これが信頼区間の統計的な意味です。
$\frac{\bar{X} – \mu}{U / \sqrt{n}}$ が $t$ 分布に従う理由
$\bar{X}$ を確率変数 $X_i \us{ \sim }_{\rm i.i.d.} \mathcal{N}(\mu, \sigma^2)$ の標本平均 $\bar{X} = \f{1}{n} \sum_{i=1}^n X_i$ とするとき、
\begin{align*}
T = \f{\bar{X} – \mu}{U / \sqrt{n}}, \t U^2 = \f{1}{n-1} \sum_{i=1}^n \(X_i – \bar{X} \)^2
\end{align*}
で定義される確率変数 $T$ は自由度 $n-1$ の $t$ 分布: $t(n-1)$に従う。
上記がなぜ成り立つのかを確認します。まず $T$ を変形します。
\begin{align*}
T &= \f{\bar{X} – \mu}{U / \sqrt{n}} \n
&=\f{\bar{X} – \mu}{\sigma / \sqrt{n}} \times \sqrt{\f{\sigma^2}{U^2}} \n
&=\f{\bar{X} – \mu}{\sigma / \sqrt{n}} \times \sqrt{\f{\sigma^2}{U^2} \cdot \f{n-1}{n-1}}
\end{align*}
ここで、$V = \f{\bar{X} – \mu}{\sigma / \sqrt{n}},\ W = (n-1)\f{U^2}{\sigma^2}$ とおきます。
\begin{align*}
T = \f{V}{\sqrt{\f{W}{n-1}}}.
\end{align*}
それでは、確率変数 $T$ が従う分布を求めるために、まずは、$V$ と $W$ が従う確率分布をそれぞれ求めます。
$V$ が従う確率分布
$V = \f{\bar{X} – \mu}{\sigma / \sqrt{n}}$ ですが、標本平均 $\bar{X} = \f{1}{n} \sum_{i=1}^n X_i$ の各々の確率変数は $X_i \us{ \sim }_{\rm i.i.d.} \mathcal{N}(\mu, \sigma^2)$ であるため、$V$ は標準正規分布 $\mathcal{N}(0, 1)$ に従います。
\begin{align*}
V = \f{\bar{X} – \mu}{\sigma / \sqrt{n}} \sim \mathcal{N}(0, 1).
\end{align*}
$W$ が従う確率分布
$W = (n-1)\f{U^2}{\sigma^2}$ が従う確率分布を求めるには、次の定理を用います。
標準正規分布に従う確率変数の組 $Y_i {\us \sim_{\rm i.i.d.}} \mathcal{N}(0, 1)$ を用いて作られる、確率変数 $X = \sum_{i=1}^n Y_i^2$ は自由度 $n$ のカイ二乗分布に従う。
\begin{align*}
X \sim \chi^2(n).
\end{align*}
上記の証明は▼こちらで行っています。
では、$W$ が従う確率分布を求めるためにまず、標準正規分布に従う確率変数 $Z_i \us{ \sim }_{\rm i.i.d.} \mathcal{N}(0, 1)$ を用いて、$\sum_{i=1}^n (Z_i – \bar{Z})^2$ が従う分布を調べます。
いきなりですが、次の等式を考えます。
\begin{align*}
\sum_{i=1}^n Z_i^2 = \sum_{i=1}^n \( Z_i – \bar{Z} \)^2 + (\sqrt{n} \bar{Z} )^2.
\end{align*}
そして、上式の積率母関数を比較します。定理1から、
\begin{align*}
&\sum_{i=1}^n Z_i^2 \sim \chi^2(n),\n
&(\sqrt{n} \bar{Z} )^2 \sim \chi^2(1), \t \because \sqrt{n} \bar{Z} \sim \mathcal{N}(0, 1).
\end{align*}
ですので、積率母関数は次式で書けます。
\begin{align*}
\( \f{1}{1 -2t} \)^{\f{n}{2}} = \( \sum_{i=1}^n \( Z_i – \bar{Z} \)^2 \s{の積率母関数} \) \cdot \( \f{1}{1 -2t} \)^{\f{1}{2}}.
\end{align*}
したがって、この等式より
\begin{align*}
\( \sum_{i=1}^n \( Z_i – \bar{Z} \)^2 \s{の積率母関数} \) = \( \f{1}{1 -2t} \)^{\f{n-1}{2}}
\end{align*}
\begin{align*}
\therefore \sum_{i=1}^n \( Z_i – \bar{Z} \)^2 \sim \chi^2 (n-1)
\end{align*}
が成り立ちます。
さて、確率変数 $X_i$ と $Z_i$ の関係は、$X_i = \mu + \sigma Z_i,\ \bar{X} = \mu + \sigma \bar{Z}$ なので、
\begin{align*}
\sum_{i=1}^n \( X_i – \bar{X} \)^2 &= \sum_{i=1}^n \( \mu + \sigma Z_i – (\mu + \sigma \bar{Z}) \)^2 \n
&= \sigma^2 \sum_{i=1}^n \( Z_i – \bar{Z} \)^2
\end{align*}
が成立します。
以上より、$W$ が従う確率分布は、
\begin{align*}
W &= \f{n-1}{\sigma^2} \cdot U^2 \n
&= \f{n-1}{\sigma^2} \cdot \f{1}{n-1} \ub{ \sum_{i=1}^n \(X_i – \bar{X} \)^2 }{=\sigma^2 \sum_{i=1}^n \( Z_i – \bar{Z} \)^2} \n
&= \sum_{i=1}^n \( Z_i – \bar{Z} \)^2
\end{align*}
したがって、$W$ は自由度 $n-1$ の カイ二乗分布に従います。
\begin{align*}
W = (n-1)\f{U^2}{\sigma^2} \sim \chi^2(n-1).
\end{align*}
$\frac{V}{\sqrt{\f{W}{n-1}}}$ が従う確率分布
$\frac{V}{\sqrt{\f{W}{n-1}}}$ が従う確率分布を求めるには、次の定理を用います。
確率変数 $Y$ が自由度 $n$ のカイ二乗分布、$Z$ が標準正規分布に従う $Y \sim \chi^2(n), Z \sim \mathcal{N}(0, 1)$ 時、確率変数 $X = \f{Z}{\sqrt{Y/n}}$ は自由度 $n$ の $t$ 分布に従う。
\begin{align*}
X = \f{Z}{\sqrt{\f{Y}{n}}} \sim t(n).
\end{align*}
上記の証明は▼こちらで行っています。
今、$V \sim \mathcal{N}(0, 1),\ W \sim \chi^2(n-1)$ ですので、定理2より $\frac{V}{\sqrt{\f{W}{n-1}}}$ は自由度 $n-1$ の$t$ 分布に従います。以上より、
\begin{align*}
\f{\bar{X} – \mu}{U / \sqrt{n}} \sim t(n-1).
\end{align*}
参考書籍
