尤度比検定

これまで扱ってきたZ検定やt検定、母比率の検定では、それぞれの問題設定に特化した検定統計量を用いてきた。しかし、統計モデルが複雑になると、個別の検定統計量を導くのは容易ではない。もっと汎用的な検定の枠組みはないだろうか？

その答えが尤度比検定（likelihood ratio test, LRT）である。最尤推定の考え方を検定に応用したこの手法は、「制約なしのモデルと制約ありのモデルで尤度がどれだけ異なるか」を比較する。1938年にアメリカの数学者サミュエル・ウィルクス（Samuel S. Wilks, 1906–1964）が漸近的な分布を明らかにして以来、パラメトリックモデルにおける仮説検定の標準的な構成法として広く用いられている。

尤度比検定とは

基本的な考え方

尤度比検定の考え方はシンプルである。2つのモデルの「データへの当てはまり」を比較する。

① 帰無仮説のもとで最もデータに合うモデル（パラメータに制約あり）
② 制約なしで最もデータに合うモデル（パラメータを自由に動かせる）

①と②の当てはまりに大差がなければ「帰無仮説の制約があっても十分データを説明できる」ので帰無仮説を採択する。逆に、制約を外すと大幅に当てはまりが改善されるなら、帰無仮説の制約は不適切と判断して棄却する。この「当てはまりの差」を尤度関数の比で測定するのが尤度比検定である。

コイン投げで理解する

具体例で見てみよう。コインを100回投げて60回表が出たとする。表が出る確率を \theta とする。

帰無仮説 H_0: \theta = 0.5 のもとでの尤度（このデータが生じる確率の大きさ）を L(0.5)、制約なしで尤度が最大となる \hat{\theta} = 0.6（最尤推定量）での尤度を L(0.6) とする。尤度比は

この比が1に近ければ「\theta = 0.5 でも十分データを説明できる」、大きければ「\theta = 0.5 ではデータをうまく説明できない」と判断する。

一般的な定義

一般に、パラメータ \boldsymbol{\theta} をもつ統計モデルの尤度関数を L(\boldsymbol{\theta}) と書くとき、帰無仮説 H_0（パラメータに何らかの制約を課す仮説）に対する尤度比は次のように定義される。

分子は全パラメータ空間で尤度を最大化したもの、分母は帰無仮説の制約のもとで最大化したものである。定義より \lambda_n \geq 1 であり、帰無仮説が正しければ \lambda_n は1に近くなる。

ウィルクスの定理

「どれくらい大きければ棄却するか」の基準

尤度比 \lambda_n が大きいほど帰無仮説に反する証拠が強いことはわかった。しかし、実際に検定を行うには「どれくらい大きければ棄却するのか」という基準が必要である。そのためには、帰無仮説が正しいときに \lambda_n がどのような分布に従うかを知る必要がある。

これに答えたのが、ウィルクスが1938年に証明した定理である。結論は非常にシンプルで、「尤度比の対数の2倍」を計算すれば、それが近似的にカイ二乗分布に従うというものである。

自由度の求め方

自由度 p は「帰無仮説がいくつのパラメータを固定しているか」で決まる。例えば、先ほどのコイン投げの例では、帰無仮説 H_0: \theta = 0.5 はパラメータ \theta を 1つの値に固定しているので p = 1 であり、\Lambda は自由度1のカイ二乗分布に従う。もし2つのパラメータ（例えば正規分布の平均と分散）を同時に固定する仮説なら p = 2 となる。

実際の検定手順

対数尤度 \ell(\cdot) = \log L(\cdot) を使うと、検定統計量は次のように書ける。

ここで \ell(\hat{\theta}) は制約なしでの最大対数尤度、\ell(\hat{\theta}_0) は帰無仮説の制約下での最大対数尤度である。この \Lambda の値がカイ二乗分布の上側 \alpha 点を超えれば、帰無仮説を棄却する。

計算例：母比率の尤度比検定

解答

X \sim \text{Bin}(100, \theta) とする。対数尤度は（定数項を除いて）

制約なしの最尤推定量は \hat{\theta} = \dfrac{60}{100} = 0.6、帰無仮説のもとでは \theta_0 = 0.5。検定統計量は

自由度1のカイ二乗分布の上側5%点は \chi^2_{0.05}(1) = 3.841 である。

帰無仮説を棄却する。このコインは公正ではないと判断される（P値 = 0.045）。

なお、同じデータに対して母比率の検定のスコア統計量を計算すると Z^2 = 4.000 であり、尤度比統計量 \Lambda = 4.027 と近い値になる。漸近的にはこれらは等価である。

計算例：正規分布の平均の尤度比検定

解答

正規分布の対数尤度は（定数項を除いて）

\hat{\mu} = \bar{x} = 52、\mu_0 = 50 を代入して

\Lambda = 4.000 > 3.841 なので帰無仮説を棄却する（P値 = 0.046）。

この結果はZ検定の検定統計量 Z = \dfrac{\bar{x} - \mu_0}{\sigma / \sqrt{n}} = \dfrac{52 - 50}{5/5} = 2.000 の2乗に一致する。すなわち、正規分布の平均に関して、尤度比検定はZ検定と完全に等価である。

3つの古典的検定

パラメトリックモデルにおける仮説検定には、次の3つの古典的なアプローチがある。いずれも漸近的にカイ二乗分布に従い、大標本では等価な結論を与える。

先ほどのコイン投げの例（n=100, x=60, \theta_0=0.5）では、スコア検定が Z^2 = 4.000、尤度比検定が \Lambda = 4.027、ワルド検定が Z^2 = 4.167 と、いずれも近い値となっている。

練習問題

まとめ

Python実装

尤度比検定は対数尤度の差から簡潔に計算できる。以下では母比率とポアソン分布のパラメータに対する尤度比検定を実装する。

import numpy as np
from scipy import stats

def binomial_lrt(x, n, theta_0):
    """二項分布の尤度比検定"""
    theta_hat = x / n
    if theta_hat == 0 or theta_hat == 1:
        return np.inf, 0.0
    log_lr = x * np.log(theta_hat / theta_0) \
           + (n - x) * np.log((1 - theta_hat) / (1 - theta_0))
    lrt_stat = 2 * log_lr
    p_value = 1 - stats.chi2.cdf(lrt_stat, df=1)
    return lrt_stat, p_value

def poisson_lrt(sum_x, n, lambda_0):
    """ポアソン分布の尤度比検定"""
    lambda_hat = sum_x / n
    log_lr = sum_x * np.log(lambda_hat / lambda_0) \
           + n * (lambda_0 - lambda_hat)
    lrt_stat = 2 * log_lr
    p_value = 1 - stats.chi2.cdf(lrt_stat, df=1)
    return lrt_stat, p_value

# ========== 例題1: コイン投げ ==========
print("【母比率の尤度比検定】")
print("n = 100, x = 60, H₀: θ = 0.5")
lrt, p = binomial_lrt(60, 100, 0.5)
print(f"Λ = {lrt:.4f}, P値 = {p:.4f}")
print(f"結論: 帰無仮説を棄却" if lrt > 3.841 else "結論: 帰無仮説を棄却しない")

# スコア検定・ワルド検定との比較
theta_hat, theta_0, n = 0.6, 0.5, 100
z_score = (theta_hat - theta_0) / np.sqrt(theta_0*(1-theta_0)/n)
z_wald = (theta_hat - theta_0) / np.sqrt(theta_hat*(1-theta_hat)/n)
print(f"\n【3つの検定の比較】")
print(f"スコア検定: Z² = {z_score**2:.4f}")
print(f"尤度比検定: Λ  = {lrt:.4f}")
print(f"ワルド検定: Z² = {z_wald**2:.4f}")

print()
print("=" * 50)
print()

# ========== 例題2: 正規分布の平均 ==========
print("【正規分布の平均の尤度比検定】")
n2, xbar, mu0, sigma = 25, 52, 50, 5
lrt2 = n2 * (xbar - mu0)**2 / sigma**2
p2 = 1 - stats.chi2.cdf(lrt2, df=1)
print(f"n = {n2}, x̄ = {xbar}, μ₀ = {mu0}, σ = {sigma}")
print(f"Λ = n(x̄ - μ₀)²/σ² = {lrt2:.4f}")
print(f"P値 = {p2:.4f}")
z = (xbar - mu0) / (sigma / np.sqrt(n2))
print(f"Z検定の Z² = {z**2:.4f}（尤度比と一致）")