機械学習と情報技術

ベイズモデル比較 — WAIC・LOO-CVの理論と実践

ベイズ統計

ベイズ統計でモデルを構築したとき、「このモデルはデータをよく説明しているのか」「2つの候補モデルのどちらが良いか」を判断する必要があります。しかし、この一見単純な問いには深い理論的背景があります。

頻度論的な統計ではAIC（赤池情報量基準）やBIC（ベイズ情報量基準）がモデル選択の定番ですが、ベイズ統計ではパラメータが事後分布として推定されるため、これらの基準をそのまま適用するには修正が必要です。そもそも「良いモデル」とは何を意味するのかという問い自体が、ベイズの枠組みでは異なる形で現れます。

ベイズモデル比較を理解すると、以下のような応用が開けます。

  • モデル選択: 複数の候補モデルから最も予測性能の高いモデルを選ぶ
  • モデル診断: 構築したモデルが新しいデータに対してどの程度の予測精度を持つか評価する
  • 正則化の度合いの調整: 事前分布のハイパーパラメータの選択基準として利用する
  • 機械学習: ベイズニューラルネットワークやガウス過程回帰モデルの比較

本記事の内容

  • ベイズモデル比較の2つのパラダイム — ベイズファクター vs 予測性能
  • 予測分布と期待対数予測密度（elpd）の理論
  • AICからWAICへ — 情報量基準のベイズ版
  • LOO-CV（Leave-One-Out交差検証）とPSIS近似の理論
  • WAICとLOO-CVの関係と使い分け
  • Pythonによる実装と比較

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

ベイズモデル比較の2つの考え方

ベイズファクターによるアプローチ

ベイズの枠組みでモデルを比較する最も古典的な方法は、ベイズファクター（Bayes factor）です。2つのモデル $M_1$ と $M_2$ に対して、ベイズファクター $B_{12}$ は

$$ \begin{equation} B_{12} = \frac{p(y | M_1)}{p(y | M_2)} = \frac{\int p(y | \theta_1, M_1) p(\theta_1 | M_1) d\theta_1}{\int p(y | \theta_2, M_2) p(\theta_2 | M_2) d\theta_2} \end{equation} $$

と定義されます。分子と分母はそれぞれのモデルの周辺尤度（marginal likelihood）であり、パラメータを事前分布で積分消去した量です。

ベイズファクターは事前確率比と掛け合わせることで事後確率比を与えます。

$$ \frac{P(M_1 | y)}{P(M_2 | y)} = B_{12} \cdot \frac{P(M_1)}{P(M_2)} $$

ベイズファクターには理論的な美しさがありますが、実用上いくつかの問題があります。

事前分布への敏感性: 周辺尤度はパラメータの事前分布に強く依存します。曖昧な（分散の大きい）事前分布を使うと、周辺尤度が小さくなる傾向があります。これは「バートレットのパラドクス」として知られています。無情報事前分布を使いたい場合に特に問題です。

計算の困難さ: 周辺尤度の計算には高次元の積分が必要であり、MCMCの出力から直接推定するのは困難です。調和平均推定量などのナイーブな方法は分散が無限大になりえます。

「M-open」問題: 真のデータ生成過程がどの候補モデルにも含まれない場合（M-open設定）、ベイズファクターによるモデル選択は必ずしも予測性能の良いモデルを選びません。

予測性能によるアプローチ

これらの問題を回避する代替アプローチが、予測性能の評価に基づくモデル比較です。「将来のデータに対してどれだけ良い予測ができるか」でモデルを比較します。

このアプローチでは、モデルの良さを期待対数予測密度（expected log pointwise predictive density, elpd）で測ります。elpd が大きいモデルほど、新しいデータに対する予測が良いことを意味します。

WAICとLOO-CVは、このelpd を効率的に推定するための2つの方法です。

ベイズファクターと予測性能アプローチの違いを理解した上で、次に予測性能の理論的な基盤を定式化しましょう。

期待対数予測密度（elpd）

予測分布の定義

ベイズ統計における予測の中心概念は事後予測分布（posterior predictive distribution）です。観測データ $y = (y_1, \ldots, y_n)$ が与えられたとき、新しいデータ点 $\tilde{y}$ に対する予測分布は

$$ \begin{equation} p(\tilde{y} | y) = \int p(\tilde{y} | \theta) \, p(\theta | y) \, d\theta \end{equation} $$

パラメータ $\theta$ を事後分布で「平均化」した予測です。点推定（最尤推定値やMAP推定値）ではなく、パラメータの不確実性を織り込んだ予測になっています。

elpdの定義

期待対数予測密度（elpd）は、事後予測分布の「良さ」を真の分布 $f(\tilde{y})$ に対して測る量です。

$$ \begin{equation} \text{elpd} = \sum_{i=1}^{n} E_f\left[\ln p(\tilde{y}_i | y)\right] = \sum_{i=1}^{n} \int f(\tilde{y}_i) \ln p(\tilde{y}_i | y) \, d\tilde{y}_i \end{equation} $$

ここで期待値は、将来のデータ $\tilde{y}_i$ が真の分布 $f$ から生成されるとして取っています。

elpdとKLダイバージェンスの関係

elpd の理論的な意味は、KLダイバージェンスとの関係で明確になります。$i$ 番目のデータ点に対して

$$ \begin{align} E_f[\ln p(\tilde{y}_i | y)] &= E_f[\ln f(\tilde{y}_i)] – D_{\text{KL}}(f \| p(\cdot | y)) \end{align} $$

右辺の第1項 $E_f[\ln f(\tilde{y}_i)]$ は真の分布のエントロピーの負で、モデルに依存しない定数です。したがって、elpd を最大化するモデルは、事後予測分布と真の分布のKLダイバージェンスを最小化するモデルと一致します。

elpdの推定問題

elpd の定義には真の分布 $f$ が含まれるため、直接計算することはできません。実際に利用可能なのは手元のデータ $y_1, \ldots, y_n$ だけです。

素朴な推定量として、トレーニングデータでの対数予測密度の合計（対数事後予測密度、lppd）を考えます。

$$ \begin{equation} \text{lppd} = \sum_{i=1}^{n} \ln p(y_i | y) = \sum_{i=1}^{n} \ln \int p(y_i | \theta) \, p(\theta | y) \, d\theta \end{equation} $$

しかし、lppd はトレーニングデータ $y_i$ を予測にも評価にも使っているため、elpd の楽観的な過大推定になります。これは頻度論におけるトレーニング誤差と汎化誤差の関係と同じ構造です。

elpd とlppd の差を有効パラメータ数 $p_{\text{eff}}$ と呼びます。

$$ \begin{equation} \text{elpd} \approx \text{lppd} – p_{\text{eff}} \end{equation} $$

WAICとLOO-CVは、この $p_{\text{eff}}$ を推定する異なる方法を提供します。

lppd の楽観的バイアスをどう補正するか——このアイデアは、実は頻度論的なAICと全く同じ発想です。次のセクションで、AICからWAICへの自然な流れを見ていきましょう。

AICからWAICへ

AICの復習

赤池情報量基準（AIC）は、最尤推定に基づくモデル選択の基準です。

$$ \begin{equation} \text{AIC} = -2 \sum_{i=1}^{n} \ln p(y_i | \hat{\theta}_{\text{MLE}}) + 2k \end{equation} $$

ここで $\hat{\theta}_{\text{MLE}}$ は最尤推定値、$k$ はパラメータの数です。第1項は当てはまりの良さ（deviance）、第2項は複雑さのペナルティです。

AICの理論的根拠は、正則条件の下で対数尤度の楽観バイアスが漸近的に $k$（パラメータ数）に等しいという結果に基づいています。

AICの限界

AICにはいくつかの限界があります。

  1. 点推定への依存: 最尤推定値 $\hat{\theta}_{\text{MLE}}$ のみを使い、パラメータの不確実性を無視する
  2. 正則条件の仮定: 尤度関数が正則である（Fisher情報行列が正則、MLEが漸近正規）ことを仮定する
  3. 事前情報の非利用: ベイズ的な事前情報を考慮しない

ベイズの枠組みでは事後分布全体を使って予測するため、AICをそのまま適用するのは不適切です。

WAIC（Widely Applicable Information Criterion）

渡辺（2010）が提案したWAIC（広く適用可能な情報量基準、Watanabe-Akaike Information Criterion）は、AICのベイズ版として位置づけられます。特に、AICが要求する正則条件を必要とせず、特異モデル（パラメータが識別不能なモデル）にも適用できます。

WAICは次のように定義されます。

$$ \begin{equation} \text{WAIC} = -2(\text{lppd} – p_{\text{WAIC}}) \end{equation} $$

ここで lppd は前述の対数事後予測密度、$p_{\text{WAIC}}$ は有効パラメータ数の推定量です。

有効パラメータ数の2つの定義

WAICの有効パラメータ数には2つの定義があります。

$p_{\text{WAIC1}}$（推定量1）:

$$ \begin{equation} p_{\text{WAIC1}} = 2\sum_{i=1}^{n}\left[\ln E_{\text{post}}[p(y_i | \theta)] – E_{\text{post}}[\ln p(y_i | \theta)]\right] \end{equation} $$

これは「対数尤度の事後期待値」と「尤度の事後期待値の対数」の差の2倍です。イェンセンの不等式から、各項は常に非負です。

$p_{\text{WAIC2}}$（推定量2）:

$$ \begin{equation} p_{\text{WAIC2}} = \sum_{i=1}^{n} \text{Var}_{\text{post}}[\ln p(y_i | \theta)] \end{equation} $$

こちらは各データ点の対数尤度の事後分散の和です。

$p_{\text{WAIC2}}$ の方が安定性が高いことが多く、実用的には $p_{\text{WAIC2}}$ が推奨されます。Gelman et al. (2014) も $p_{\text{WAIC2}}$ を推奨しています。

WAICの直感的な理解

$p_{\text{WAIC2}}$ の直感的な意味を考えましょう。

あるデータ点 $y_i$ に対して、対数尤度 $\ln p(y_i | \theta)$ が事後分布の下で大きく変動する（分散が大きい）場合、そのデータ点に対する予測はパラメータの値に強く依存しており、「有効パラメータ」を多く消費しています。

逆に、対数尤度の分散が小さい場合は、どのパラメータ値でも似たような予測を行うため、事実上パラメータを使っていないのと同じです。

この解釈は、AICのペナルティ $k$ よりも精緻です。AICではパラメータの数を数えるだけですが、WAICでは各パラメータが実際にどの程度使われているかを考慮します。強い事前分布でパラメータが拘束されている場合、有効パラメータ数はパラメータの実数よりも小さくなります。

MCMCからWAICを計算する

WAICの実用的な利点は、MCMCの事後サンプルから直接計算できることです。$S$ 個の事後サンプル $\theta^{(1)}, \ldots, \theta^{(S)}$ があるとき

$$ \text{lppd} \approx \sum_{i=1}^{n} \ln\left(\frac{1}{S}\sum_{s=1}^{S} p(y_i | \theta^{(s)})\right) $$

$$ p_{\text{WAIC2}} \approx \sum_{i=1}^{n} \text{Var}_s\left[\ln p(y_i | \theta^{(s)})\right] $$

ここで $\text{Var}_s$ はサンプルに基づく分散です。

周辺尤度の計算を必要としないため、ベイズファクターよりもはるかに簡単に計算できます。

WAICが「トレーニングデータでの予測密度からバイアスを引く」という情報量基準的なアプローチだとすれば、LOO-CVは「1つのデータ点を除いて予測し、その精度を直接評価する」という交差検証的なアプローチです。

LOO-CV（Leave-One-Out 交差検証）

理論的な定義

Leave-One-Out交差検証（LOO-CV）は、最も直接的な予測性能の推定方法です。各データ点 $y_i$ を1つずつ取り除き、残りのデータ $y_{-i} = (y_1, \ldots, y_{i-1}, y_{i+1}, \ldots, y_n)$ でモデルを推定して $y_i$ を予測します。

LOO-CVの推定量は次のように定義されます。

$$ \begin{equation} \text{elpd}_{\text{loo}} = \sum_{i=1}^{n} \ln p(y_i | y_{-i}) = \sum_{i=1}^{n} \ln \int p(y_i | \theta) \, p(\theta | y_{-i}) \, d\theta \end{equation} $$

ここで $p(\theta | y_{-i})$ は $y_i$ を除いたデータに基づく事後分布です。

LOO-CVの魅力は、各データ点 $y_i$ が評価のためにのみ使われ、予測モデルの構築には使われないことです。これにより、lppd のような楽観バイアスを回避できます。

計算の困難さと解決策

素朴にLOO-CVを計算するには、$n$ 回のMCMCを実行する必要があります。$n = 1000$ のデータに対して1回のMCMCに1分かかるとすれば、約16時間かかります。これは明らかに実用的ではありません。

この計算コストの問題を解決するのがPSIS-LOO（Pareto Smoothed Importance Sampling LOO）です。

PSIS-LOOの理論

PSIS-LOOは、重点サンプリング（importance sampling）を用いて、完全データの事後分布 $p(\theta | y)$ からのサンプルを使って、各 $p(\theta | y_{-i})$ からのサンプルを近似します。

完全データの事後分布からのサンプル $\theta^{(s)} \sim p(\theta | y)$ に対して、重点重み（importance weight）は

$$ \begin{equation} w_i^{(s)} = \frac{p(\theta^{(s)} | y_{-i})}{p(\theta^{(s)} | y)} \propto \frac{1}{p(y_i | \theta^{(s)})} \end{equation} $$

最後の比例関係は、$p(\theta | y) \propto p(y_i | \theta) p(\theta | y_{-i})$ から得られます。$y_i$ を除いた事後分布は、完全データの事後分布を尤度 $p(y_i | \theta)$ で「割り戻す」ことで得られるのです。

これにより、LOO-CVの各項は次のように推定できます。

$$ \begin{equation} p(y_i | y_{-i}) \approx \frac{\sum_{s=1}^{S} w_i^{(s)} p(y_i | \theta^{(s)})}{\sum_{s=1}^{S} w_i^{(s)}} \end{equation} $$

あるいは、より安定な推定として

$$ \begin{equation} p(y_i | y_{-i}) \approx \frac{1}{\frac{1}{S}\sum_{s=1}^{S} \frac{1}{p(y_i | \theta^{(s)})}} \end{equation} $$

重点重みの問題とパレート平滑化

素朴な重点サンプリング推定量には問題があります。いくつかの重み $w_i^{(s)} = 1/p(y_i | \theta^{(s)})$ が非常に大きくなることがあり、推定量の分散が爆発します。これは特に、データ点 $y_i$ がモデルにとって「外れ値的」な場合に起こります。

Vehtari, Gelman, Gabry (2017) が提案したパレート平滑化は、この問題を次のように解決します。

  1. 各データ点 $i$ について、$S$ 個の重み $w_i^{(1)}, \ldots, w_i^{(S)}$ を計算する
  2. 重みの上位部分（上位20%程度）を一般化パレート分布にフィットさせる
  3. フィットしたパレート分布の量子関数を用いて、極端に大きな重みを「平滑化」（truncation + smoothing）する

フィットしたパレート分布の形状パラメータ $\hat{k}$ は、推定の信頼性を示す診断指標として利用できます。

$\hat{k}$ の範囲 解釈
$\hat{k} < 0.5$ 良好。推定量は安定
$0.5 \leq \hat{k} < 0.7$ 注意が必要。推定量はやや不安定
$0.7 \leq \hat{k} < 1$ 問題あり。推定量の分散が大きい可能性
$\hat{k} \geq 1$ 信頼できない。推定量の分散が無限大の可能性

$\hat{k}$ が大きいデータ点は、そのデータ点がモデルにとって「影響力の大きい外れ値」であることを示唆しています。これはモデル診断としても有用な情報です。

WAICとLOO-CVの関係

漸近的な等価性

WAICとLOO-CVは、漸近的には同じ量を推定しています。データ数 $n$ が大きいとき

$$ \text{WAIC} \approx -2 \cdot \text{elpd}_{\text{loo}} $$

つまり、両者は漸近的に等価です。

しかし有限サンプルでは差が生じます。特に以下の場合にLOO-CVの方が信頼性が高い傾向があります。

影響力の大きいデータ点がある場合: LOO-CVは各データ点を直接除外して評価するため、外れ値の影響をより適切に捉えます。WAICの $p_{\text{WAIC}}$ は局所的な近似に基づくため、外れ値に対して不安定になることがあります。

モデルが誤特定されている場合: 真のデータ生成過程がモデルに含まれない場合、WAICの漸近近似の精度が低下します。

事後分布が正規近似で良く近似されない場合: WAICは暗黙的に事後分布の局所的な性質を使っているため、多峰性の事後分布などでは不安定になります。

実用的な使い分け

特性 WAIC PSIS-LOO
計算コスト 低い（MCMCサンプルの後処理のみ） やや高い（パレート平滑化の追加処理）
信頼性 一般的に良好、ただし外れ値に弱い 一般的に WAIC より安定
診断情報 なし $\hat{k}$ 診断が利用可能
ソフトウェア 多くのベイズパッケージで利用可能 ArviZ, loo パッケージ

実用的には、PSIS-LOOを第一選択とし、$\hat{k}$ 診断で推定の信頼性を確認するのが推奨されます。$\hat{k}$ が大きいデータ点があれば、そのデータ点がモデルに合っていないことを示唆しているため、モデルの改善のヒントにもなります。

WAICとLOO-CVの理論的な関係を理解したところで、次にPythonでの具体的な実装を通じて、これらの量がどのように計算されるかを確認しましょう。

Pythonによる実装

線形回帰モデルでの比較

ここでは、真の関係が2次関数であるデータに対して、1次・2次・3次の多項式回帰モデルを当てはめ、WAIC と LOO-CV によるモデル比較を行います。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# --- 合成データの生成（真の関係: y = 2 + 0.5x - 0.3x^2 + noise）---
n = 50
x = np.random.uniform(-3, 3, n)
y_true = 2 + 0.5 * x - 0.3 * x**2
sigma_true = 1.0
y = y_true + np.random.normal(0, sigma_true, n)

# 可視化
fig, ax = plt.subplots(figsize=(8, 5))
x_plot = np.linspace(-3.5, 3.5, 200)
ax.scatter(x, y, color="steelblue", alpha=0.7, s=40, label="Data")
ax.plot(x_plot, 2 + 0.5*x_plot - 0.3*x_plot**2, "r--", linewidth=2,
        label="True function")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_title("Synthetic Data (True: Quadratic)", fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("model_comparison_data.png", dpi=150, bbox_inches="tight")
plt.show()

データの散布図と真の関数を比較すると、2次関数（放物線）の周りにノイズが載った構造が確認できます。正しいモデル選択基準を使えば、2次モデルが選ばれるはずです。1次モデルは過小適合（underfitting）、高次モデルは過剰適合（overfitting）のリスクがあります。

# --- ベイズ線形回帰（共役事前分布）---
def bayesian_linear_regression(X, y, n_samples=4000):
    """正規-逆ガンマ事前分布による共役ベイズ線形回帰"""
    n, p = X.shape
    # 曖昧な事前分布
    mu_0 = np.zeros(p)
    Lambda_0 = 0.01 * np.eye(p)  # 精度行列（弱い事前分布）
    a_0, b_0 = 1.0, 1.0  # 逆ガンマのパラメータ

    # 事後パラメータ
    Lambda_n = Lambda_0 + X.T @ X
    Sigma_n = np.linalg.inv(Lambda_n)
    mu_n = Sigma_n @ (Lambda_0 @ mu_0 + X.T @ y)
    a_n = a_0 + n / 2
    b_n = b_0 + 0.5 * (y @ y + mu_0 @ Lambda_0 @ mu_0 - mu_n @ Lambda_n @ mu_n)

    # 事後サンプルの生成
    sigma2_samples = 1.0 / np.random.gamma(a_n, 1.0 / b_n, n_samples)
    beta_samples = np.array([
        np.random.multivariate_normal(mu_n, s2 * Sigma_n)
        for s2 in sigma2_samples
    ])

    return beta_samples, sigma2_samples

# --- 3つのモデルの設計行列 ---
models = {
    "Linear (degree 1)": np.column_stack([np.ones(n), x]),
    "Quadratic (degree 2)": np.column_stack([np.ones(n), x, x**2]),
    "Cubic (degree 3)": np.column_stack([np.ones(n), x, x**2, x**3]),
}

results = {}
for name, X in models.items():
    beta_samples, sigma2_samples = bayesian_linear_regression(X, y)
    results[name] = {
        "X": X,
        "beta_samples": beta_samples,
        "sigma2_samples": sigma2_samples,
    }

# --- WAIC の計算 ---
def compute_waic(X, y, beta_samples, sigma2_samples):
    """WAICを計算する"""
    S = len(sigma2_samples)
    n = len(y)

    # 各サンプル・各データ点の対数尤度
    log_lik = np.zeros((S, n))
    for s in range(S):
        mu_pred = X @ beta_samples[s]
        for i in range(n):
            log_lik[s, i] = stats.norm.logpdf(y[i], mu_pred[i],
                                               np.sqrt(sigma2_samples[s]))

    # lppd: 対数事後予測密度
    # log(mean(exp(log_lik))) を数値安定に計算
    max_log_lik = np.max(log_lik, axis=0)
    lppd_i = max_log_lik + np.log(np.mean(np.exp(log_lik - max_log_lik), axis=0))
    lppd = np.sum(lppd_i)

    # p_WAIC2: 有効パラメータ数（分散ベース）
    p_waic = np.sum(np.var(log_lik, axis=0, ddof=1))

    # WAIC
    waic = -2 * (lppd - p_waic)

    return waic, lppd, p_waic, log_lik

# --- PSIS-LOO の計算 ---
def compute_loo(log_lik, k_threshold=0.7):
    """簡易版 PSIS-LOO を計算する"""
    S, n = log_lik.shape

    loo_i = np.zeros(n)
    k_hat = np.zeros(n)

    for i in range(n):
        # 重点重み（対数スケール）
        log_w = -log_lik[:, i]  # w_i^(s) = 1/p(y_i|theta^(s))
        log_w -= np.max(log_w)
        w = np.exp(log_w)

        # パレート平滑化の簡易版:
        # 上位20%の重みをソートしてパレート分布をフィット
        M = max(int(0.2 * S), 5)
        sorted_w = np.sort(w)
        tail = sorted_w[-M:]

        if np.min(tail) > 0 and len(tail) > 2:
            # 簡易的なパレート形状パラメータの推定
            # (Method of Moments)
            log_tail = np.log(tail / tail[0])
            k_hat[i] = np.mean(log_tail)
        else:
            k_hat[i] = 0.0

        # 正規化された重み
        w_norm = w / np.sum(w)

        # LOO予測密度
        loo_i[i] = np.log(np.sum(w_norm * np.exp(log_lik[:, i])))

    elpd_loo = np.sum(loo_i)
    loo = -2 * elpd_loo

    return loo, elpd_loo, k_hat

# 全モデルについて計算
print(f"{'Model':<25} {'WAIC':>8} {'p_WAIC':>8} {'LOO':>8} {'p_eff':>8}")
print("-" * 60)

for name, res in results.items():
    waic, lppd, p_waic, log_lik = compute_waic(
        res["X"], y, res["beta_samples"], res["sigma2_samples"])
    loo, elpd_loo, k_hat = compute_loo(log_lik)

    results[name]["waic"] = waic
    results[name]["p_waic"] = p_waic
    results[name]["loo"] = loo
    results[name]["lppd"] = lppd
    results[name]["elpd_loo"] = elpd_loo
    results[name]["k_hat"] = k_hat

    print(f"{name:<25} {waic:>8.2f} {p_waic:>8.2f} {loo:>8.2f} "
          f"{lppd - elpd_loo:>8.2f}")

# --- 結果の可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) WAIC と LOO の比較
ax = axes[0]
model_names = list(results.keys())
waics = [results[m]["waic"] for m in model_names]
loos = [results[m]["loo"] for m in model_names]

x_pos = np.arange(len(model_names))
width = 0.35
bars1 = ax.bar(x_pos - width/2, waics, width, color="steelblue", alpha=0.8,
               label="WAIC")
bars2 = ax.bar(x_pos + width/2, loos, width, color="coral", alpha=0.8,
               label="LOO")
ax.set_ylabel("Information Criterion (lower is better)", fontsize=11)
ax.set_title("WAIC vs LOO-CV", fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels(["Linear", "Quadratic", "Cubic"], fontsize=10)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")

# (b) 有効パラメータ数
ax = axes[1]
p_waics = [results[m]["p_waic"] for m in model_names]
n_params = [2, 3, 4]  # 実際のパラメータ数

ax.bar(x_pos - width/2, n_params, width, color="lightgray", alpha=0.8,
       edgecolor="gray", label="Actual params")
ax.bar(x_pos + width/2, p_waics, width, color="steelblue", alpha=0.8,
       label=r"$p_{\mathrm{WAIC}}$")
ax.set_ylabel("Number of parameters", fontsize=11)
ax.set_title("Actual vs Effective Parameters", fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels(["Linear", "Quadratic", "Cubic"], fontsize=10)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")

# (c) 予測曲線の比較
ax = axes[2]
ax.scatter(x, y, color="gray", alpha=0.5, s=20)
x_plot = np.linspace(-3.5, 3.5, 200)
colors = ["#e74c3c", "#2ecc71", "#3498db"]

for j, (name, res) in enumerate(results.items()):
    X_plot = np.column_stack([x_plot**d for d in range(res["X"].shape[1])])
    beta_mean = np.mean(res["beta_samples"], axis=0)
    y_pred = X_plot @ beta_mean
    ax.plot(x_plot, y_pred, color=colors[j], linewidth=2, label=name.split(" ")[0])

ax.plot(x_plot, 2 + 0.5*x_plot - 0.3*x_plot**2, "k--", linewidth=1.5,
        label="True", alpha=0.5)
ax.set_xlabel("x", fontsize=11)
ax.set_ylabel("y", fontsize=11)
ax.set_title("Model Predictions", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-3.5, 3.5)

plt.tight_layout()
plt.savefig("model_comparison_results.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから、WAICとLOO-CVによるモデル比較の結果が読み取れます。

  1. 左図（WAIC vs LOO-CV）: 両基準とも2次モデル（Quadratic）が最小の値を示しており、真のモデルを正しく選択しています。WAICとLOO-CVの値はほぼ一致しており、漸近的な等価性が確認できます。1次モデルは過小適合によるペナルティ、3次モデルは若干の過剰適合によるペナルティを受けています

  2. 中央図（有効パラメータ数）: $p_{\text{WAIC}}$ は実際のパラメータ数にほぼ等しい値を示しています。弱い事前分布（曖昧な事前分布）を使っているため、事前分布によるパラメータの拘束がほとんどなく、有効パラメータ数がフルのパラメータ数に近くなっています。強い事前分布を使えば $p_{\text{WAIC}}$ はパラメータの実数より小さくなります

  3. 右図（予測曲線）: 1次モデルは直線的な関係しか捉えられず、データの湾曲を見逃しています。2次モデルは真の関数に非常に近い予測曲線を描きます。3次モデルはデータ範囲内では2次モデルとほぼ同じ予測をしていますが、余分なパラメータのペナルティを受けています

BICとの比較

BICの定義と性質

AIC/WAICとは異なるモデル選択の哲学を持つのがBIC（ベイズ情報量基準）です。

$$ \begin{equation} \text{BIC} = -2 \ln p(y | \hat{\theta}_{\text{MLE}}) + k \ln n \end{equation} $$

BICのペナルティ $k \ln n$ はAICの $2k$ より厳しく（$n > 7$ のとき）、より単純なモデルを選ぶ傾向があります。

BICは真のモデルの選択（model selection consistency）を目標としています。つまり、データ数 $n \to \infty$ のとき、候補モデルの中に真のモデルが含まれていれば、BICは確率1で真のモデルを選びます。

一方、AIC/WAICは予測性能の最大化を目標としています。真のモデルが候補に含まれていなくても（M-open設定）、最も予測が良いモデルを選びます。

目的に応じた使い分け

目的 推奨基準
真のモデルを特定したい（因果推論、科学的解釈） BIC、ベイズファクター
予測精度を最大化したい（予測、意思決定） WAIC、LOO-CV
モデルの誤特定が疑われる（M-open設定） LOO-CV（診断付き）

実用的には、両方の基準が同じモデルを選べば安心ですが、異なるモデルを選ぶ場合は目的に応じて判断します。

実用上の注意点

WAICの数値安定性

WAICの lppd の計算では、対数尤度のlog-sum-exp計算が必要です。素朴に実装すると数値オーバーフローやアンダーフローが起きるため、log-sum-exp のトリックを使う必要があります。

$$ \ln\left(\frac{1}{S}\sum_{s=1}^{S} p(y_i | \theta^{(s)})\right) = \max_s \ln p(y_i | \theta^{(s)}) + \ln\left(\frac{1}{S}\sum_{s=1}^{S} \exp\left(\ln p(y_i | \theta^{(s)}) – \max_s \ln p(y_i | \theta^{(s)})\right)\right) $$

$\hat{k}$ 診断の活用

PSIS-LOOの $\hat{k}$ 診断は、モデル改善のヒントを与えてくれます。$\hat{k} > 0.7$ のデータ点がある場合

  1. そのデータ点が外れ値かどうかを確認する
  2. モデルの裾の重さが不十分かどうか確認する（正規分布 → t分布への変更など）
  3. そのデータ点の影響を調べるためにケース削除分析を行う

MCMCサンプル数の影響

WAICもLOO-CVも、MCMCの事後サンプルに基づく推定量です。サンプル数 $S$ が少ないと推定の精度が低下します。一般的には $S \geq 1000$（理想的には $S \geq 4000$）が推奨されます。

特にWAICの $p_{\text{WAIC2}}$ は分散の推定を含むため、サンプル数の影響を受けやすくなります。

まとめ

本記事では、ベイズモデル比較の理論と実践について解説しました。

  • elpd（期待対数予測密度）はベイズモデルの予測性能を測る理論的な基盤であり、KLダイバージェンスの最小化と等価である
  • WAICは lppd（対数事後予測密度）から有効パラメータ数 $p_{\text{WAIC}}$ を引くことで elpd を推定する情報量基準である
  • LOO-CVは各データ点を除外して予測性能を直接評価する交差検証法であり、PSIS による近似で効率的に計算できる
  • PSIS-LOO の $\hat{k}$ 診断は推定の信頼性とモデルの適切さの両方を評価できる
  • WAICとLOO-CVは漸近的に等価だが、有限サンプルではPSIS-LOO の方が安定している
  • BICは真のモデル選択を目的とし、WAIC/LOO-CVは予測性能の最大化を目的とする。目的に応じた使い分けが重要

次のステップとして、以下の記事も参考にしてください。

関連記事

前の記事

経験ベイズ法の理論と応用 — データから事前分布を推定する

関連タグ: LOO-CV WAIC ベイズ統計 モデル比較 予測性能 情報量基準