【統計検定対策】区間推定、信頼区間の求め方とpythonによる実装

確率・統計学

はじめに

前回までは、確率分布の母平均の点推定について焦点を当てていました。

母平均の最尤推定量である標本平均は一般には誤差を含み、母平均に一致しません。そこで母平均を高い確率で含む区間を求めることにします。

そのようなパラメータを含む区間を求める技術を区間推定と言います。

この記事では区間推定の仕方と、$t$ 分布による信頼区間の計算方法を述べます。

また、この記事で使用しているコードは以下のGoogle Colabで試すことができます。

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*}

▲標準正規分布の場合、区間 $[-1.96, 1.96]$ の面積は全体の $95\%$ となる

この区間

\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 \%$ の確率で正解であることを意味します。

このように、真の推定値が含まれる区間を確率的に推定することを区間推定と言います。

信頼区間について、 「母平均 $\mu$ がこの区間に収まる確率が $95\%$ である」とは言えません。母平均はあくまでも定まった値(定数)であり、確率変数ではありません。確率変数であるのは信頼区間に含まれている標本平均 $\bar{x}_n$ ですので、サンプリングの際、変動するのは信頼区間の方で母平均ではありません。 後ほど、シミュレーションを通して確認します。

ただし、上記の信頼区間 $\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$ 分布に従うかはこの記事の最後に述べます)。

$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$ 分布の理論を用いることにより、信頼区間を計算することができます。

$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()
▲サンプルサイズ $n = 5$ として $\mathcal{N}(50, 10)$ の正規乱数から、$\f{\bar{X} – \mu}{U / \sqrt{n}}$ を$20000$ 回計算し、ヒストグラムにした。この場合、$\f{\bar{X} – \mu}{U / \sqrt{n}}$ は自由度 $4$ の $t$ 分布に従う。それは、ヒストグラムが赤線(自由度 $4$ の $t$ 分布)にフィットしていることから確かめられる。

上図は赤線が自由度 $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)$ に従うとして、下記を実行します。

  1. 母集団から $100$ 個サンプリングし$95\%$ 信頼区間を計算する。
  2. 1の試行を $10000$ 回繰り返す。
  3. $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\%$ となります。これが信頼区間の統計的な意味です。

▲$95\%$ 信頼区間は母平均を含んだり、含まなかったりする様子。縦棒が $95\%$ 信頼区間を表しており、赤が母平均を含む、青が母平均を含まないを意味する。信頼区間は確率変数であるため、各試行で信頼区間は異なり、無限の試行では $95\%$ 信頼区間が母平均を含まない確率は $5\%$ となる。

$\frac{\bar{X} – \mu}{U / \sqrt{n}}$ が $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 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}$ が従う確率分布を求めるには、次の定理を用います。

定理1

標準正規分布に従う確率変数の組 $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}}}$ が従う確率分布を求めるには、次の定理を用います。

定理2

確率変数 $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*}


参考書籍

タイトルとURLをコピーしました