統計的仮説検定を行うとき、「どのような検定統計量を構成すれば最も効率よく帰無仮説と対立仮説を区別できるか」は根本的な問いです。この問いに対する理論的に最も美しい回答を与えるのが 尤度比検定（Likelihood Ratio Test, LRT） です。
尤度比検定は、帰無仮説のもとでの最大尤度と、制約のない場合の最大尤度の比に基づいて仮説を検定します。ネイマン・ピアソンの基本補題により、単純仮説同士の検定では尤度比検定が 最強力検定 であることが示されます。さらに、ウィルクスの定理により、大標本では $-2\log\Lambda$ が $\chi^2$ 分布に従うことが保証されるため、複合仮説に対しても広く適用可能です。
本記事の内容
- 尤度比の定義と直感的な理解
- ネイマン・ピアソンの基本補題（最強力検定の証明スケッチ）
- 一般化尤度比検定統計量 $\Lambda$ の構成
- ウィルクスの定理（$-2\log\Lambda \to \chi^2$）の導出
- 正規分布の平均検定への適用
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
尤度比とは
直感的な理解
手元にデータ $\bm{x} = (x_1, x_2, \dots, x_n)$ があり、2つの仮説を比べたいとします。
- $H_0$: パラメータは $\theta = \theta_0$
- $H_1$: パラメータは $\theta = \theta_1$
どちらの仮説のもとでデータがより「もっともらしい」かを測るには、それぞれの仮説のもとでの 尤度（likelihood） を比べるのが自然です。尤度とは、パラメータを固定したときにデータが得られる確率（密度）のことでした。
$$ \text{尤度比} = \frac{L(\theta_0 \mid \bm{x})}{L(\theta_1 \mid \bm{x})} $$
この比が小さいほど、帰無仮説 $H_0$ のもとでのデータのもっともらしさが、対立仮説 $H_1$ に比べて低いことを意味します。
尤度関数の復習
データ $\bm{x} = (x_1, \dots, x_n)$ が独立同分布に従うとき、尤度関数は
$$ L(\theta \mid \bm{x}) = \prod_{i=1}^{n} f(x_i \mid \theta) $$
対数尤度は
$$ \ell(\theta \mid \bm{x}) = \sum_{i=1}^{n} \log f(x_i \mid \theta) $$
です。
ネイマン・ピアソンの基本補題
問題設定
単純仮説 同士の検定を考えます。
$$ H_0: \theta = \theta_0 \quad \text{vs} \quad H_1: \theta = \theta_1 $$
有意水準 $\alpha$ の検定のうち、検出力（$H_1$ が真のとき正しく $H_0$ を棄却する確率）を最大にする検定を 最強力検定（Most Powerful Test） と呼びます。
補題の主張
ネイマン・ピアソンの基本補題: 有意水準 $\alpha$ の最強力検定は、尤度比
$$ \Lambda_{\text{NP}}(\bm{x}) = \frac{L(\theta_0 \mid \bm{x})}{L(\theta_1 \mid \bm{x})} $$
がある閾値 $k$ 以下のとき $H_0$ を棄却する検定、すなわち
$$ \text{棄却域} = \{\bm{x} : \Lambda_{\text{NP}}(\bm{x}) \leq k\} $$
です。ここで $k$ は
$$ P(\Lambda_{\text{NP}}(\bm{X}) \leq k \mid H_0) = \alpha $$
を満たすように選びます。
証明のスケッチ
尤度比検定の棄却域を $R^* = \{\bm{x} : \Lambda_{\text{NP}}(\bm{x}) \leq k\}$ とし、同じ有意水準 $\alpha$ をもつ任意の他の検定の棄却域を $R$ とします。
最強力性を示すには、$P(\bm{X} \in R^* \mid H_1) \geq P(\bm{X} \in R \mid H_1)$ を証明すればよいです。
ステップ1: 棄却域の構成から、$\bm{x} \in R^*$ のとき $L(\theta_0 \mid \bm{x}) \leq k \cdot L(\theta_1 \mid \bm{x})$ です。$\bm{x} \notin R^*$ のとき $L(\theta_0 \mid \bm{x}) > k \cdot L(\theta_1 \mid \bm{x})$ です。
したがって、任意の $\bm{x}$ に対して
$$ \left[\mathbf{1}_{R^*}(\bm{x}) – \mathbf{1}_R(\bm{x})\right] \left[L(\theta_1 \mid \bm{x}) – \frac{1}{k} L(\theta_0 \mid \bm{x})\right] \geq 0 $$
が成り立ちます。ここで $\mathbf{1}_A(\bm{x})$ は集合 $A$ の指示関数です。
これを確認しましょう。
- $\bm{x} \in R^*$ のとき: $L(\theta_1 \mid \bm{x}) \geq \frac{1}{k} L(\theta_0 \mid \bm{x})$ なので、$\mathbf{1}_{R^*} – \mathbf{1}_R$ が1または0であれば積は $\geq 0$
- $\bm{x} \notin R^*$ のとき: $L(\theta_1 \mid \bm{x}) < \frac{1}{k} L(\theta_0 \mid \bm{x})$ なので、$\mathbf{1}_{R^*} - \mathbf{1}_R$ が0または $-1$ であれば積は $\geq 0$
ステップ2: 上の不等式を $\bm{x}$ の全空間で積分します。
$$ \int \left[\mathbf{1}_{R^*}(\bm{x}) – \mathbf{1}_R(\bm{x})\right] L(\theta_1 \mid \bm{x}) \, d\bm{x} \geq \frac{1}{k} \int \left[\mathbf{1}_{R^*}(\bm{x}) – \mathbf{1}_R(\bm{x})\right] L(\theta_0 \mid \bm{x}) \, d\bm{x} $$
左辺は
$$ P(\bm{X} \in R^* \mid H_1) – P(\bm{X} \in R \mid H_1) $$
右辺は
$$ \frac{1}{k}\left[P(\bm{X} \in R^* \mid H_0) – P(\bm{X} \in R \mid H_0)\right] = \frac{1}{k}(\alpha – \alpha) = 0 $$
ここで両方の検定が有意水準 $\alpha$ であることを使いました。
ステップ3: したがって
$$ P(\bm{X} \in R^* \mid H_1) – P(\bm{X} \in R \mid H_1) \geq 0 $$
すなわち、尤度比検定の検出力が他の任意の検定の検出力以上であることが示されました。
補題の意義
ネイマン・ピアソンの基本補題は、単純仮説 vs 単純仮説の枠組みでは、尤度比に基づく検定が最適である ことを保証します。この結果は、より一般的な検定を構成する際の理論的基盤となります。
一般化尤度比検定
複合仮説への拡張
実際の検定問題では、帰無仮説や対立仮説がパラメータの特定の値ではなく、パラメータ空間の部分集合として指定されることが多いです。
パラメータ空間を $\Theta$、帰無仮説に対応する部分空間を $\Theta_0 \subset \Theta$ とします。
$$ H_0: \theta \in \Theta_0 \quad \text{vs} \quad H_1: \theta \in \Theta \setminus \Theta_0 $$
一般化尤度比統計量
一般化尤度比統計量（Generalized Likelihood Ratio Statistic） は次のように定義されます。
$$ \Lambda(\bm{x}) = \frac{\sup_{\theta \in \Theta_0} L(\theta \mid \bm{x})}{\sup_{\theta \in \Theta} L(\theta \mid \bm{x})} $$
分子は帰無仮説の制約のもとでの最大尤度、分母は制約のない場合の最大尤度です。
最尤推定量（MLE）を用いて書くと、
$$ \Lambda(\bm{x}) = \frac{L(\hat{\theta}_0 \mid \bm{x})}{L(\hat{\theta} \mid \bm{x})} $$
ここで $\hat{\theta}_0 = \arg\max_{\theta \in \Theta_0} L(\theta \mid \bm{x})$ は制約付きMLE、$\hat{\theta} = \arg\max_{\theta \in \Theta} L(\theta \mid \bm{x})$ は制約なしMLEです。
$\Lambda$ の性質
- $0 \leq \Lambda(\bm{x}) \leq 1$（分子は分母以下）
- $\Lambda$ が小さいほど、帰無仮説のもとでの尤度が相対的に低い
- 棄却域は $\Lambda(\bm{x}) \leq c$（$c$ は有意水準から決定）
対数尤度比統計量
実用上は $\Lambda$ そのものではなく、$-2\log\Lambda$ を検定統計量として用います。
$$ -2\log\Lambda(\bm{x}) = -2\left[\ell(\hat{\theta}_0 \mid \bm{x}) – \ell(\hat{\theta} \mid \bm{x})\right] = 2\left[\ell(\hat{\theta} \mid \bm{x}) – \ell(\hat{\theta}_0 \mid \bm{x})\right] $$
$-2\log\Lambda \geq 0$ であり、大きいほど帰無仮説を棄却する方向です。
ウィルクスの定理
定理の主張
ウィルクスの定理（Wilks’ theorem）: 正則条件のもとで、帰無仮説 $H_0: \theta \in \Theta_0$ が真のとき、
$$ -2\log\Lambda(\bm{X}) \xrightarrow{d} \chi^2_r \quad (n \to \infty) $$
ここで $r = \dim(\Theta) – \dim(\Theta_0)$ は帰無仮説によって課される制約の数です。
正則条件
ウィルクスの定理が成り立つには、以下のような正則条件が必要です。
- パラメータ空間 $\Theta$ が開集合
- 真のパラメータ $\theta_0$ が $\Theta_0$ の内点
- 対数尤度関数が3回連続微分可能
- フィッシャー情報行列 $\bm{I}(\theta)$ が正定値
導出
パラメータ $\bm{\theta} = (\theta_1, \dots, \theta_p)$ の全体空間を $\Theta \subset \mathbb{R}^p$ とし、帰無仮説が $\theta_{q+1} = \theta_{q+1}^{(0)}, \dots, \theta_p = \theta_p^{(0)}$（$r = p – q$ 個の制約）の形をしているとします。
ステップ1: 対数尤度のテイラー展開
真のパラメータ値 $\bm{\theta}_0$ のまわりで、対数尤度 $\ell(\bm{\theta})$ を2次までテイラー展開します。
$$ \ell(\bm{\theta}) \approx \ell(\bm{\theta}_0) + (\bm{\theta} – \bm{\theta}_0)^\top \nabla\ell(\bm{\theta}_0) + \frac{1}{2}(\bm{\theta} – \bm{\theta}_0)^\top \nabla^2\ell(\bm{\theta}_0) (\bm{\theta} – \bm{\theta}_0) $$
ステップ2: MLEの漸近正規性
制約なしMLE $\hat{\bm{\theta}}$ の漸近正規性より、
$$ \sqrt{n}(\hat{\bm{\theta}} – \bm{\theta}_0) \xrightarrow{d} N(\bm{0}, \bm{I}(\bm{\theta}_0)^{-1}) $$
ここで $\bm{I}(\bm{\theta}_0)$ は1標本あたりのフィッシャー情報行列です。
また、MLEでは $\nabla\ell(\hat{\bm{\theta}}) = \bm{0}$ が成り立ちます。
ステップ3: 対数尤度比の2次近似
MLEのまわりでテイラー展開すると、
$$ \ell(\hat{\bm{\theta}}_0) \approx \ell(\hat{\bm{\theta}}) + \frac{1}{2}(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}})^\top \nabla^2\ell(\hat{\bm{\theta}})(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}}) $$
1次の項は $\nabla\ell(\hat{\bm{\theta}}) = \bm{0}$ なので消えます。
したがって、
$$ \begin{align} -2\log\Lambda &= 2[\ell(\hat{\bm{\theta}}) – \ell(\hat{\bm{\theta}}_0)] \\ &\approx -(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}})^\top \nabla^2\ell(\hat{\bm{\theta}})(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}}) \\ &= (\hat{\bm{\theta}}_0 – \hat{\bm{\theta}})^\top \left[-\nabla^2\ell(\hat{\bm{\theta}})\right](\hat{\bm{\theta}}_0 – \hat{\bm{\theta}}) \end{align} $$
ステップ4: 観測情報行列への収束
大数の法則により、観測情報行列は期待情報行列に収束します。
$$ -\frac{1}{n}\nabla^2\ell(\hat{\bm{\theta}}) \xrightarrow{p} \bm{I}(\bm{\theta}_0) $$
すなわち $-\nabla^2\ell(\hat{\bm{\theta}}) \approx n\bm{I}(\bm{\theta}_0)$ です。
ステップ5: 2次形式のχ²分布への収束
$\hat{\bm{\theta}}$ と $\hat{\bm{\theta}}_0$ の差は、$H_0$ のもとで $r$ 次元の制約を課したことの影響のみを反映します。
$\bm{Z} = \sqrt{n}\bm{I}(\bm{\theta}_0)^{1/2}(\hat{\bm{\theta}} – \bm{\theta}_0) \xrightarrow{d} N(\bm{0}, \bm{I}_p)$ とおくと、
$$ -2\log\Lambda \approx n(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}})^\top \bm{I}(\bm{\theta}_0)(\hat{\bm{\theta}}_0 – \hat{\bm{\theta}}) $$
帰無仮説のもとでは、$\hat{\bm{\theta}}_0$ と $\hat{\bm{\theta}}$ の差は $r$ 個の制約方向にのみ生じます。したがって、
$$ -2\log\Lambda \xrightarrow{d} \chi^2_r $$
ここで $r = p – q = \dim(\Theta) – \dim(\Theta_0)$ です。
定理の直感的な意味
ウィルクスの定理は、「制約の数だけ自由度をもつ $\chi^2$ 分布に収束する」ことを述べています。帰無仮説が $r$ 個のパラメータを固定すると、対数尤度比は $r$ 次元の正規ベクトルの2次形式として表され、これが $\chi^2_r$ に従います。
正規分布の平均検定への適用
問題設定
$X_1, X_2, \dots, X_n \sim N(\mu, \sigma^2)$（$\sigma^2$ は未知）として、
$$ H_0: \mu = \mu_0 \quad \text{vs} \quad H_1: \mu \neq \mu_0 $$
を検定します。
パラメータ空間
- 全体: $\Theta = \{(\mu, \sigma^2) : \mu \in \mathbb{R}, \sigma^2 > 0\}$、$\dim(\Theta) = 2$
- 帰無仮説: $\Theta_0 = \{(\mu_0, \sigma^2) : \sigma^2 > 0\}$、$\dim(\Theta_0) = 1$
- 制約の数: $r = 2 – 1 = 1$
尤度関数
$$ L(\mu, \sigma^2 \mid \bm{x}) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i – \mu)^2}{2\sigma^2}\right) $$
対数尤度は
$$ \ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi) – \frac{n}{2}\log\sigma^2 – \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i – \mu)^2 $$
制約なしMLE
$\partial\ell/\partial\mu = 0$ より $\hat{\mu} = \bar{x}$
$\partial\ell/\partial\sigma^2 = 0$ より $\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i – \bar{x})^2$
最大対数尤度は
$$ \ell(\hat{\mu}, \hat{\sigma}^2) = -\frac{n}{2}\log(2\pi) – \frac{n}{2}\log\hat{\sigma}^2 – \frac{n}{2} $$
制約付きMLE（$\mu = \mu_0$ の制約）
$\mu = \mu_0$ として $\sigma^2$ のみを最適化します。
$\partial\ell/\partial\sigma^2 = 0$ より $\hat{\sigma}_0^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i – \mu_0)^2$
最大対数尤度は
$$ \ell(\mu_0, \hat{\sigma}_0^2) = -\frac{n}{2}\log(2\pi) – \frac{n}{2}\log\hat{\sigma}_0^2 – \frac{n}{2} $$
尤度比統計量の導出
$$ \begin{align} -2\log\Lambda &= 2[\ell(\hat{\mu}, \hat{\sigma}^2) – \ell(\mu_0, \hat{\sigma}_0^2)] \\ &= 2\left[\left(-\frac{n}{2}\log\hat{\sigma}^2 – \frac{n}{2}\right) – \left(-\frac{n}{2}\log\hat{\sigma}_0^2 – \frac{n}{2}\right)\right] \\ &= 2 \cdot \left(-\frac{n}{2}\right)(\log\hat{\sigma}^2 – \log\hat{\sigma}_0^2) \\ &= n\log\frac{\hat{\sigma}_0^2}{\hat{\sigma}^2} \end{align} $$
ここで、
$$ \hat{\sigma}_0^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i – \mu_0)^2 $$
$$ \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i – \bar{x})^2 $$
$\sum(x_i – \mu_0)^2$ を分解します。
$$ \begin{align} \sum_{i=1}^{n}(x_i – \mu_0)^2 &= \sum_{i=1}^{n}[(x_i – \bar{x}) + (\bar{x} – \mu_0)]^2 \\ &= \sum_{i=1}^{n}(x_i – \bar{x})^2 + 2(\bar{x} – \mu_0)\sum_{i=1}^{n}(x_i – \bar{x}) + n(\bar{x} – \mu_0)^2 \\ &= \sum_{i=1}^{n}(x_i – \bar{x})^2 + n(\bar{x} – \mu_0)^2 \end{align} $$
3行目で、$\sum_{i=1}^{n}(x_i – \bar{x}) = 0$ を使いました。
したがって、
$$ \frac{\hat{\sigma}_0^2}{\hat{\sigma}^2} = \frac{\sum(x_i – \bar{x})^2 + n(\bar{x} – \mu_0)^2}{\sum(x_i – \bar{x})^2} = 1 + \frac{n(\bar{x} – \mu_0)^2}{\sum(x_i – \bar{x})^2} $$
$s^2 = \frac{1}{n-1}\sum(x_i – \bar{x})^2$（不偏分散）とおくと、$\sum(x_i – \bar{x})^2 = (n-1)s^2$ ですから、
$$ \frac{\hat{\sigma}_0^2}{\hat{\sigma}^2} = 1 + \frac{n(\bar{x} – \mu_0)^2}{(n-1)s^2} = 1 + \frac{t^2}{n-1} $$
ここで $t = \frac{\bar{x} – \mu_0}{s/\sqrt{n}}$ はt統計量です。
最終的に、
$$ -2\log\Lambda = n\log\left(1 + \frac{t^2}{n-1}\right) $$
$n$ が大きいとき、$\log(1 + x) \approx x$（$x$ が小さいとき）の近似を用いると、
$$ -2\log\Lambda \approx n \cdot \frac{t^2}{n-1} \approx t^2 $$
$t$ 分布は自由度が大きいとき標準正規分布に近づくので、$t^2$ は $\chi^2_1$ に近似的に従います。これはウィルクスの定理（$r = 1$）と整合しています。
具体例
数値例
工場で製造される部品の長さが、設計値 $\mu_0 = 100$ mm であるかを検定します。$n = 30$ 個の部品を測定し、$\bar{x} = 101.2$, $s = 3.0$ でした。
尤度比統計量の計算:
$$ t = \frac{\bar{x} – \mu_0}{s/\sqrt{n}} = \frac{101.2 – 100}{3.0/\sqrt{30}} = \frac{1.2}{0.5477} = 2.191 $$
$$ -2\log\Lambda = 30 \cdot \log\left(1 + \frac{2.191^2}{29}\right) = 30 \cdot \log\left(1 + 0.1655\right) = 30 \times 0.1533 = 4.600 $$
ウィルクスの定理より $-2\log\Lambda \sim \chi^2_1$ として p 値を計算すると、
$$ p = P(\chi^2_1 > 4.600) \approx 0.032 $$
有意水準5%で帰無仮説を棄却します。
参考: $t^2 = 2.191^2 = 4.800$ と比較すると、$-2\log\Lambda = 4.600$ は近い値です。$n$ が大きくなるほど両者は一致していきます。
Pythonでの実装
尤度比検定のスクラッチ実装
import numpy as np
from scipy import stats
def likelihood_ratio_test_normal_mean(data, mu_0):
"""
正規分布の平均に関する尤度比検定（分散未知）
Parameters
----------
data : array-like
標本データ
mu_0 : float
帰無仮説の母平均
Returns
-------
lrt_stat : float
-2 log Λ（尤度比統計量）
p_value_chi2 : float
χ²近似によるp値
t_stat : float
対応するt統計量
p_value_t : float
t検定によるp値
"""
data = np.array(data, dtype=float)
n = len(data)
x_bar = np.mean(data)
s2 = np.var(data, ddof=1) # 不偏分散
# t統計量
t_stat = (x_bar - mu_0) / np.sqrt(s2 / n)
# 尤度比統計量
lrt_stat = n * np.log(1 + t_stat**2 / (n - 1))
# χ²近似によるp値（ウィルクスの定理、自由度1）
p_value_chi2 = 1 - stats.chi2.cdf(lrt_stat, df=1)
# t検定によるp値（参考値）
p_value_t = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n - 1))
return lrt_stat, p_value_chi2, t_stat, p_value_t
# --- 具体例 ---
np.random.seed(42)
mu_0 = 100 # 帰無仮説の母平均
n = 30
# データ生成（真の平均は101.2）
data = np.random.normal(101.2, 3.0, n)
lrt_stat, p_chi2, t_stat, p_t = likelihood_ratio_test_normal_mean(data, mu_0)
print("=== 正規分布の平均に関する尤度比検定 ===")
print(f"標本サイズ: n = {n}")
print(f"標本平均: x̄ = {np.mean(data):.4f}")
print(f"標本標準偏差: s = {np.std(data, ddof=1):.4f}")
print(f"帰無仮説: μ₀ = {mu_0}")
print()
print(f"t統計量: t = {t_stat:.4f}")
print(f"t² = {t_stat**2:.4f}")
print(f"-2 log Λ = {lrt_stat:.4f}")
print()
print(f"χ²近似によるp値: {p_chi2:.4f}")
print(f"t検定によるp値: {p_t:.4f}")
print(f"判定 (α=0.05): {'棄却' if p_chi2 < 0.05 else '棄却しない'}")
一般化尤度比検定（2つの正規分布の平均の同時検定）
import numpy as np
from scipy import stats
def lrt_two_normal_means(data1, data2):
"""
2つの独立な正規分布の平均が等しいかの尤度比検定
H0: μ1 = μ2 (共通分散は未知)
Parameters
----------
data1, data2 : array-like
各群の標本データ
Returns
-------
lrt_stat : float
-2 log Λ
p_value : float
χ²近似によるp値
"""
x = np.array(data1, dtype=float)
y = np.array(data2, dtype=float)
n1, n2 = len(x), len(y)
n = n1 + n2
# 制約なしMLE
mu1_hat = np.mean(x)
mu2_hat = np.mean(y)
sigma2_hat = (np.sum((x - mu1_hat)**2) + np.sum((y - mu2_hat)**2)) / n
# 制約付きMLE (μ1 = μ2 = μ)
mu_hat = (n1 * mu1_hat + n2 * mu2_hat) / n
sigma2_hat_0 = (np.sum((x - mu_hat)**2) + np.sum((y - mu_hat)**2)) / n
# 対数尤度比統計量
lrt_stat = n * np.log(sigma2_hat_0 / sigma2_hat)
# p値（χ²近似、自由度1）
p_value = 1 - stats.chi2.cdf(lrt_stat, df=1)
return lrt_stat, p_value
# --- 具体例 ---
np.random.seed(42)
group_a = np.random.normal(50, 5, 25)
group_b = np.random.normal(53, 5, 30)
lrt_stat, p_value = lrt_two_normal_means(group_a, group_b)
# t検定との比較
t_stat, p_t = stats.ttest_ind(group_a, group_b)
print("=== 2群の平均の尤度比検定 ===")
print(f"群A: n₁ = {len(group_a)}, x̄₁ = {np.mean(group_a):.4f}")
print(f"群B: n₂ = {len(group_b)}, x̄₂ = {np.mean(group_b):.4f}")
print()
print(f"-2 log Λ = {lrt_stat:.4f}")
print(f"χ²近似 p値: {p_value:.4f}")
print()
print(f"t検定 t統計量: {t_stat:.4f}")
print(f"t検定 p値: {p_t:.4f}")
ウィルクスの定理の検証シミュレーション
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n_simulations = 10000
sample_sizes = [10, 30, 100, 500]
mu_0 = 0 # 帰無仮説のもとでの真の値
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()
for idx, n in enumerate(sample_sizes):
lrt_stats = []
for _ in range(n_simulations):
# H0のもとでデータ生成（μ=0, σ=1）
data = np.random.normal(mu_0, 1.0, n)
x_bar = np.mean(data)
s2 = np.var(data, ddof=1)
# t統計量
t = (x_bar - mu_0) / np.sqrt(s2 / n)
# 尤度比統計量
lrt = n * np.log(1 + t**2 / (n - 1))
lrt_stats.append(lrt)
lrt_stats = np.array(lrt_stats)
# ヒストグラム
x_range = np.linspace(0, 12, 200)
axes[idx].hist(lrt_stats, bins=80, density=True, alpha=0.6,
color='steelblue', label='$-2\\log\\Lambda$ (simulated)')
# 理論的なχ²(1)分布
axes[idx].plot(x_range, stats.chi2.pdf(x_range, df=1), 'r-',
linewidth=2.5, label='$\\chi^2_1$ (theoretical)')
axes[idx].set_title(f'n = {n}', fontsize=13)
axes[idx].set_xlabel('$-2\\log\\Lambda$', fontsize=12)
axes[idx].set_ylabel('Density', fontsize=12)
axes[idx].legend(fontsize=10)
axes[idx].set_xlim(0, 12)
axes[idx].set_ylim(0, 3)
axes[idx].grid(True, alpha=0.3)
plt.suptitle("Wilks' Theorem: Convergence of $-2\\log\\Lambda$ to $\\chi^2_1$ under $H_0$",
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
第1種の過誤率の確認
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n_simulations = 20000
sample_sizes = np.arange(5, 205, 5)
alpha = 0.05
mu_0 = 0
type1_lrt = []
type1_t = []
for n in sample_sizes:
reject_lrt = 0
reject_t = 0
for _ in range(n_simulations):
# H0のもとでデータ生成
data = np.random.normal(mu_0, 1.0, n)
x_bar = np.mean(data)
s2 = np.var(data, ddof=1)
t_stat = (x_bar - mu_0) / np.sqrt(s2 / n)
# 尤度比検定（χ²近似）
lrt = n * np.log(1 + t_stat**2 / (n - 1))
p_lrt = 1 - stats.chi2.cdf(lrt, df=1)
if p_lrt < alpha:
reject_lrt += 1
# t検定（正確な分布）
p_t = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=n - 1))
if p_t < alpha:
reject_t += 1
type1_lrt.append(reject_lrt / n_simulations)
type1_t.append(reject_t / n_simulations)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(sample_sizes, type1_lrt, 'b-o', markersize=4, linewidth=1.5,
label='LRT ($\\chi^2$ approximation)')
ax.plot(sample_sizes, type1_t, 'r-s', markersize=4, linewidth=1.5,
label='t-test (exact)')
ax.axhline(y=alpha, color='gray', linestyle='--', linewidth=1.5,
label=f'Nominal $\\alpha$ = {alpha}')
ax.set_xlabel('Sample Size $n$', fontsize=13)
ax.set_ylabel('Type I Error Rate', fontsize=13)
ax.set_title('Type I Error Rate: LRT vs t-test', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0.03, 0.08)
plt.tight_layout()
plt.show()
まとめ
本記事では、尤度比検定の理論と導出について解説しました。
- 尤度比 $\Lambda = L(\hat{\theta}_0)/L(\hat{\theta})$ は、帰無仮説と対立仮説のもとでのデータのもっともらしさの比
- ネイマン・ピアソンの基本補題: 単純仮説同士の検定では尤度比検定が最強力検定
- 一般化尤度比統計量: 複合仮説に対して $\Lambda = \sup_{\Theta_0} L / \sup_{\Theta} L$ と定義
- ウィルクスの定理: $-2\log\Lambda \xrightarrow{d} \chi^2_r$（$r$ は制約の数）
- 正規分布の平均検定では $-2\log\Lambda = n\log(1 + t^2/(n-1))$ となり、t検定と漸近的に等価
- 小標本では $\chi^2$ 近似が不正確になるため、正確なt検定を使うのが安全
次のステップとして、以下の記事も参考にしてください。