ベイズモデル平均化(BMA) — 複数モデルの予測を最適に統合する

ある問題に対してモデルを構築するとき、複数の候補モデルが同程度にデータをよく説明することがあります。たとえば、散布図に多項式を当てはめるとき、2次式がよいのか、3次式がよいのか、4次式がよいのかは一目瞭然とは限りません。

AIC や BIC などの情報量基準を使って「ベストモデル」を1つ選び、それだけで予測するのは最も一般的なアプローチです。しかし、この方法には見落とされがちな問題があります。選ばれなかったモデルの方が実は正しい可能性を完全に無視しているのです。モデル選択自体に不確実性があるにもかかわらず、1つのモデルだけに賭けてしまうのは、パラメータ推定で点推定だけを使い、事後分布を無視するのと同じ構造の問題です。

ベイズモデル平均化（Bayesian Model Averaging, BMA）は、この問題に対するベイズ的な解答です。パラメータの不確実性を事後分布で扱うように、モデルの不確実性もベイズの枠組みで扱うのがBMAの核心的なアイデアです。

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

ロバストな予測: 単一モデルの選択ミスに頑健な予測を実現できる。気象予測や経済予測で重要
モデルの不確実性の定量化: 「どのモデルが正しいかわからない」という不確実性を予測に正直に反映できる
科学的推論: 複数の理論仮説に対する支持の度合いを定量的に評価できる
アンサンブル学習の理論的基盤: ランダムフォレストやブースティングなどのアンサンブル手法の理論的な裏付けを理解できる

本記事の内容

単一モデル選択の問題点
BMAの直感的理解と数学的定式化
周辺尤度とモデル事後確率の計算
BICによる周辺尤度の近似
Python実装（多項式回帰で次数の異なるモデルをBMA）
モデル選択（AIC/BIC）との理論的比較

前提知識

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

ベイズ予測分布とは？事後分布から未知データを予測する理論と実装 — 予測分布の基本理論
ベイズ推定とは？仕組みについてわかりやすく解説 — 事前分布・事後分布の基本
ベイズモデル比較 — WAIC・LOO-CVの理論と実践 — 周辺尤度とベイズファクター

単一モデル選択の問題点

「ベストモデル」に全てを賭ける危険性

統計モデリングの実務では、複数の候補モデルから1つを選ぶ「モデル選択」が頻繁に行われます。典型的な手順は以下の通りです。

複数の候補モデル（例: 1次式、2次式、3次式）を用意する
各モデルにデータを当てはめる
AIC/BIC などの基準でモデルを評価する
最良のモデルを選び、そのモデルだけで予測を行う

この手順の問題はステップ4にあります。たとえば、3つのモデルのBIC値がそれぞれ $-120, -118, -115$ だったとします。最良のモデル（BIC = $-120$）が選ばれますが、2番目のモデル（BIC = $-118$）との差はわずかです。この小さな差は、データを少し変えれば逆転するかもしれません。にもかかわらず、最良モデルだけで予測すると、「2番目のモデルが正しい可能性」を完全に無視することになります。

具体例で見る問題

この問題を具体的に考えてみましょう。ある温度センサーのキャリブレーションデータがあり、真の温度 $T$ とセンサーの読み値 $V$ の関係をモデル化したいとします。

モデル A（線形）: $V = a_0 + a_1 T$、モデル事後確率 30%
モデル B（2次）: $V = b_0 + b_1 T + b_2 T^2$、モデル事後確率 50%
モデル C（3次）: $V = c_0 + c_1 T + c_2 T^2 + c_3 T^3$、モデル事後確率 20%

単一モデル選択ではモデルBだけを使いますが、モデルAとCにも合計50%の確率があります。特に外挿領域（訓練データの範囲外）では、モデル間の予測の違いが大きくなるため、この50%の無視は深刻な影響を与えます。

もう一つの問題は、予測の不確実性が過小評価されることです。単一モデルの予測区間は、そのモデルが正しいことを前提としています。しかし、モデルが間違っている可能性を考慮すると、実際の予測の不確実性はもっと大きいはずです。BMAはこの追加の不確実性を正しく反映します。

このように、単一モデル選択には本質的な限界があります。では、モデルの不確実性をどのように扱えばよいのでしょうか。その答えがベイズモデル平均化です。

BMAの直感的理解

日常のアナロジー — 複数の専門家の意見を聞く

BMAの考え方は、日常生活でも自然に行っていることと同じです。たとえば、体調が悪いときに3人の医師に診てもらったとします。

医師A（信頼度 50%）: 「風邪です。安静にしてください」
医師B（信頼度 30%）: 「アレルギーです。薬を飲んでください」
医師C（信頼度 20%）: 「ストレスです。休暇を取ってください」

あなたが合理的な患者なら、最も信頼できる医師Aの意見だけに従うのではなく、3人の意見を信頼度に応じて考慮するでしょう。BMAはまさにこの戦略を数学的に定式化したものです。

パラメータ平均化との類似性

ベイズ予測分布では、パラメータ $\theta$ について積分することで、パラメータの不確実性を織り込みました。

$$ p(y_{\text{new}} | \mathcal{D}, M) = \int p(y_{\text{new}} | \theta, M) \, p(\theta | \mathcal{D}, M) \, d\theta $$

BMAでは、これに加えてモデル $M$ についても平均するのです。

$$ p(y_{\text{new}} | \mathcal{D}) = \sum_{k=1}^{K} p(y_{\text{new}} | \mathcal{D}, M_k) \, p(M_k | \mathcal{D}) $$

つまり、BMAは二重の平均化を行っています。パラメータの不確実性（各モデル内の積分）とモデルの不確実性（モデル間の重み付き平均）の両方を処理するのです。

この直感を数学的に定式化していきましょう。

BMAの数学的定式化

モデル事後確率

$K$ 個の候補モデル $M_1, M_2, \dots, M_K$ を考えます。データ $\mathcal{D}$ を観測した後の各モデルの事後確率は、ベイズの定理により、

$$ \begin{equation} p(M_k | \mathcal{D}) = \frac{p(\mathcal{D} | M_k) \, p(M_k)}{\sum_{j=1}^{K} p(\mathcal{D} | M_j) \, p(M_j)} \end{equation} $$

ここで、

$p(M_k)$: モデル $M_k$ の事前確率。特に理由がなければ $p(M_k) = 1/K$（一様事前）とすることが多い
$p(\mathcal{D} | M_k)$: モデル $M_k$ の周辺尤度（model evidence）。モデル $M_k$ のパラメータ $\theta_k$ を事前分布で積分消去したもの

周辺尤度は次のように計算されます。

$$ \begin{equation} p(\mathcal{D} | M_k) = \int p(\mathcal{D} | \theta_k, M_k) \, p(\theta_k | M_k) \, d\theta_k \end{equation} $$

この積分は、「モデル $M_k$ とその事前分布の下で、観測データ $\mathcal{D}$ がどの程度”ありそう”か」を測っています。

BMA予測分布

BMAの予測分布は、全てのモデルの予測分布をモデル事後確率で重み付き平均したものです。

$$ \begin{equation} p(y_{\text{new}} | \mathcal{D}) = \sum_{k=1}^{K} p(y_{\text{new}} | \mathcal{D}, M_k) \, p(M_k | \mathcal{D}) \end{equation} $$

各項 $p(y_{\text{new}} | \mathcal{D}, M_k)$ はモデル $M_k$ の事後予測分布であり、パラメータの不確実性は既に積分で処理されています。

BMA予測の平均と分散

BMA予測の平均値は、各モデルの予測平均の重み付き平均です。

$$ \begin{equation} \mathbb{E}[y_{\text{new}} | \mathcal{D}] = \sum_{k=1}^{K} \mathbb{E}[y_{\text{new}} | \mathcal{D}, M_k] \, p(M_k | \mathcal{D}) \end{equation} $$

BMA予測の分散は、全分散の法則（law of total variance）により次のように分解されます。

$$ \begin{equation} \text{Var}[y_{\text{new}} | \mathcal{D}] = \underbrace{\sum_{k=1}^{K} \text{Var}[y_{\text{new}} | \mathcal{D}, M_k] \, p(M_k | \mathcal{D})}_{\text{モデル内分散の加重平均}} + \underbrace{\sum_{k=1}^{K} \left(\mathbb{E}[y_{\text{new}} | \mathcal{D}, M_k] – \mathbb{E}[y_{\text{new}} | \mathcal{D}]\right)^2 p(M_k | \mathcal{D})}_{\text{モデル間分散}} \end{equation} $$

この分解は非常に示唆的です。第1項は各モデルの予測の不確実性（パラメータの不確実性 + ノイズ）の加重平均であり、第2項はモデル間の予測の食い違いから来る追加の不確実性です。単一モデル選択では第2項が完全に無視されるのです。

全分散の法則によるこの分解を直感的に理解するために、先ほどの医師の例に戻りましょう。3人の医師が「治るまでの日数」を予測するとき、BMAの予測分散は「各医師の予測のばらつき（不確実性）」と「医師間の意見の相違」の両方を含んでいます。医師が全員同じ予測をするなら第2項はゼロですが、意見が割れているなら追加の不確実性があるのは当然です。

ここまでBMAの理論的枠組みを整理しました。しかし、実際に使うには周辺尤度を計算する必要があります。一般的にこの計算は困難ですが、BICを使った便利な近似が存在します。次にその近似を見ていきましょう。

周辺尤度の計算とBICによる近似

周辺尤度の計算の困難さ

周辺尤度 $p(\mathcal{D} | M_k) = \int p(\mathcal{D} | \theta_k, M_k) p(\theta_k | M_k) d\theta_k$ の計算は、一般に解析的に実行できません。共役事前分布を使える場合を除いて、高次元の積分が必要になります。

数値的な計算方法としては、以下のようなものがあります。

ラプラス近似: 事後分布を最頻値周りでガウス分布で近似する
調和平均推定量: MCMCサンプルから周辺尤度を推定する（ただし分散が無限大になりうるため非推奨）
ブリッジサンプリング: 事前分布と事後分布を橋渡しする方法
ネステッドサンプリング: 宇宙物理学から来た手法で、周辺尤度の計算に特化

ラプラス近似とBIC

最も広く使われる近似は、ラプラス近似に基づくBIC近似です。事後分布が最頻値 $\hat{\theta}_k$（MAP推定値）の周りでガウス分布に近いとき、周辺尤度は以下のように近似できます。

ラプラス近似の導出を見てみましょう。対数周辺尤度は、

$$ \log p(\mathcal{D} | M_k) = \log \int \exp\left[\log p(\mathcal{D} | \theta_k, M_k) + \log p(\theta_k | M_k)\right] d\theta_k $$

$\log p(\mathcal{D} | \theta_k, M_k) + \log p(\theta_k | M_k)$ を $\hat{\theta}_k$ の周りで2次までテイラー展開すると、被積分関数はガウス関数になります。このガウス積分を実行すると、

$$ \log p(\mathcal{D} | M_k) \approx \log p(\mathcal{D} | \hat{\theta}_k, M_k) + \log p(\hat{\theta}_k | M_k) + \frac{d_k}{2}\log(2\pi) – \frac{1}{2}\log|\bm{H}_k| $$

ここで $d_k$ はモデル $M_k$ のパラメータ数、$\bm{H}_k$ はヘシアン行列（対数事後分布の負の2次微分行列）です。

さらに、$N$ が大きいとき $\log p(\hat{\theta}_k | M_k)$ は $O(1)$ であり、$\log |\bm{H}_k| \approx d_k \log N + O(1)$ と近似できるので、

$$ \begin{equation} \log p(\mathcal{D} | M_k) \approx \log p(\mathcal{D} | \hat{\theta}_k, M_k) – \frac{d_k}{2}\log N \end{equation} $$

これがBIC（ベイズ情報量基準）に $-1/2$ を掛けたものに他なりません。通常のBICの定義は、

$$ \begin{equation} \text{BIC}_k = -2 \log p(\mathcal{D} | \hat{\theta}_k, M_k) + d_k \log N \end{equation} $$

なので、$\log p(\mathcal{D} | M_k) \approx -\frac{1}{2}\text{BIC}_k$ です。

BICを使ったモデル事後確率の近似

BIC近似を使うと、モデル事後確率は次のように計算できます（モデルの事前確率を一様とした場合）。

$$ p(M_k | \mathcal{D}) \approx \frac{\exp(-\frac{1}{2}\text{BIC}_k)}{\sum_{j=1}^{K} \exp(-\frac{1}{2}\text{BIC}_j)} $$

これはソフトマックス関数と同じ形式であり、BICの値が小さいモデルほど大きな重みを受け取ります。ただし、最良モデルだけに重みが集中するのではなく、BICの差が小さいモデルにも相応の重みが割り当てられます。

BICの差 $\Delta\text{BIC} = \text{BIC}_k – \text{BIC}_{\min}$ がモデル事後確率にどの程度影響するかの目安を示します。

$\Delta\text{BIC}$	モデル事後確率の比（近似）	解釈
0 – 2	> 0.37	差がほぼない
2 – 6	0.05 – 0.37	弱い証拠
6 – 10	0.007 – 0.05	強い証拠
> 10	< 0.007	非常に強い証拠

BICによる近似は簡便ですが、あくまで漸近近似であることに注意してください。サンプルサイズが小さい場合や、事後分布がガウス的でない場合は精度が低下します。

理論的な準備が整ったところで、Pythonでの実装に移りましょう。多項式回帰の次数選択問題を題材にして、BMAの効果を確認します。

Python実装 — 多項式回帰でのBMA

データ生成とモデルの準備

まず、非線形な真の関数からデータを生成し、次数の異なる多項式モデルでBMAを実行します。

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

# 真の関数
def true_function(x):
    return np.sin(2 * x) + 0.3 * x

# データ生成
np.random.seed(42)
N = 25
sigma_noise = 0.3
x_train = np.sort(np.random.uniform(0, 4, N))
y_train = true_function(x_train) + np.random.normal(0, sigma_noise, N)

# 候補モデル: 多項式 (次数 1〜8)
degrees = list(range(1, 9))

次に、各多項式モデルに対してベイズ線形回帰を実行し、予測分布と周辺尤度（BIC近似）を計算する関数を定義します。

def polynomial_basis(x, degree):
    """多項式基底関数（正規化付き）"""
    X = np.column_stack([x**i for i in range(degree + 1)])
    return X

def bayesian_linear_regression(x_train, y_train, x_test, degree, alpha=0.01, beta=None):
    """ベイズ線形回帰: 予測分布とBICを返す"""
    Phi_train = polynomial_basis(x_train, degree)
    Phi_test = polynomial_basis(x_test, degree)
    N, d = Phi_train.shape

    if beta is None:
        # ノイズ精度の推定（最尤推定ベース）
        w_mle = np.linalg.lstsq(Phi_train, y_train, rcond=None)[0]
        residuals = y_train - Phi_train @ w_mle
        sigma2_mle = np.mean(residuals**2)
        beta = 1.0 / max(sigma2_mle, 1e-6)

    # 事前分布
    S_0_inv = alpha * np.eye(d)
    m_0 = np.zeros(d)

    # 事後分布
    S_N_inv = S_0_inv + beta * Phi_train.T @ Phi_train
    S_N = np.linalg.inv(S_N_inv)
    m_N = S_N @ (S_0_inv @ m_0 + beta * Phi_train.T @ y_train)

    # 予測分布
    y_mean = Phi_test @ m_N
    y_var = np.array([phi @ S_N @ phi + 1.0/beta for phi in Phi_test])
    y_std = np.sqrt(y_var)

    # BIC計算
    w_mle = np.linalg.lstsq(Phi_train, y_train, rcond=None)[0]
    residuals = y_train - Phi_train @ w_mle
    log_lik = -N/2 * np.log(2*np.pi) - N/2 * np.log(np.mean(residuals**2)) - N/2
    bic = -2 * log_lik + d * np.log(N)

    return y_mean, y_std, bic, m_N, S_N

# テストデータ
x_test = np.linspace(-0.5, 5, 200)

# 各モデルの予測分布とBICを計算
results = {}
bic_values = []

for deg in degrees:
    y_mean, y_std, bic, m_N, S_N = bayesian_linear_regression(
        x_train, y_train, x_test, deg)
    results[deg] = {'mean': y_mean, 'std': y_std, 'bic': bic}
    bic_values.append(bic)

bic_values = np.array(bic_values)

# モデル事後確率（BIC近似）
log_weights = -0.5 * bic_values
log_weights -= logsumexp(log_weights)  # 正規化
model_weights = np.exp(log_weights)

print("=== モデル事後確率 (BIC近似) ===\n")
for deg, w in zip(degrees, model_weights):
    print(f"  次数 {deg}: BIC = {results[deg]['bic']:.2f}, "
          f"事後確率 = {w:.4f} ({w*100:.1f}%)")

このコードの出力を確認すると、各次数のモデルに対する事後確率が得られます。真の関数 $\sin(2x) + 0.3x$ は無限次のテイラー展開を持つため、低次の多項式でも中程度の近似が得られますが、完全な表現には高次が必要です。一方、高次の多項式はパラメータ数が多いためBICのペナルティが大きくなります。このトレードオフの結果として、中程度の次数に最も大きな事後確率が割り当てられます。

続いて、BMA予測分布を計算し、単一モデル選択と比較しましょう。

BMA予測の計算と可視化

# BMA予測分布の計算
bma_mean = np.zeros_like(x_test)
bma_var = np.zeros_like(x_test)

for deg, w in zip(degrees, model_weights):
    r = results[deg]
    bma_mean += w * r['mean']

for deg, w in zip(degrees, model_weights):
    r = results[deg]
    # 全分散の法則: モデル内分散 + モデル間分散
    bma_var += w * (r['std']**2 + (r['mean'] - bma_mean)**2)

bma_std = np.sqrt(bma_var)

# 最良モデル（BIC最小）
best_deg = degrees[np.argmin(bic_values)]
best_result = results[best_deg]

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 左上: BMA予測 vs 最良モデル予測
ax = axes[0, 0]
ax.plot(x_test, bma_mean, 'b-', linewidth=2, label='BMA予測')
ax.fill_between(x_test, bma_mean - 1.96*bma_std, bma_mean + 1.96*bma_std,
                alpha=0.15, color='blue', label='BMA 95%区間')
ax.plot(x_test, best_result['mean'], 'r--', linewidth=2,
        label=f'最良モデル (次数{best_deg})')
ax.fill_between(x_test, best_result['mean'] - 1.96*best_result['std'],
                best_result['mean'] + 1.96*best_result['std'],
                alpha=0.1, color='red', label=f'最良モデル 95%区間')
ax.plot(x_test, true_function(x_test), 'k:', linewidth=1.5, label='真の関数')
ax.scatter(x_train, y_train, c='gray', s=30, zorder=5, label='訓練データ')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('BMA予測 vs 単一モデル予測')
ax.legend(fontsize=9, loc='upper left')
ax.set_ylim(-3, 5)

# 右上: モデル事後確率
ax = axes[0, 1]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(degrees)))
bars = ax.bar(degrees, model_weights, color=colors, edgecolor='white')
ax.set_xlabel('多項式の次数')
ax.set_ylabel('モデル事後確率')
ax.set_title('モデル事後確率 (BIC近似)')
ax.set_xticks(degrees)
for bar, w in zip(bars, model_weights):
    if w > 0.01:
        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01,
                f'{w:.2f}', ha='center', va='bottom', fontsize=10)

# 左下: 各モデルの予測（重み付き）
ax = axes[1, 0]
for deg, w in zip(degrees, model_weights):
    if w > 0.01:
        r = results[deg]
        ax.plot(x_test, r['mean'], linewidth=1 + 3*w, alpha=0.3 + 0.7*w,
                label=f'次数{deg} (w={w:.2f})')
ax.plot(x_test, bma_mean, 'k-', linewidth=2.5, label='BMA予測')
ax.plot(x_test, true_function(x_test), 'k:', linewidth=1.5, label='真の関数')
ax.scatter(x_train, y_train, c='gray', s=20, zorder=5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('各モデルの予測と重み')
ax.legend(fontsize=9)
ax.set_ylim(-3, 5)

# 右下: 予測分散の分解
ax = axes[1, 1]
within_var = np.zeros_like(x_test)
between_var = np.zeros_like(x_test)
for deg, w in zip(degrees, model_weights):
    r = results[deg]
    within_var += w * r['std']**2
    between_var += w * (r['mean'] - bma_mean)**2

ax.fill_between(x_test, 0, np.sqrt(within_var), alpha=0.4, color='blue',
                label='モデル内不確実性')
ax.fill_between(x_test, np.sqrt(within_var),
                np.sqrt(within_var + between_var),
                alpha=0.4, color='orange', label='モデル間不確実性')
ax.plot(x_test, bma_std, 'k-', linewidth=2, label='合計 (BMA標準偏差)')
ax.set_xlabel('x')
ax.set_ylabel('標準偏差')
ax.set_title('予測不確実性の分解')
ax.legend(fontsize=10)

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

4つのグラフから、BMAの効果を多角的に理解できます。

左上（BMA予測 vs 最良モデル予測）: BMA予測（青）は、特に外挿領域（$x > 4$ や $x < 0$）で最良モデルの予測（赤）とは異なる振る舞いを示します。BMAの予測区間は最良モデルの区間よりも広く、特に外挿領域で顕著です。これは、複数のモデルの予測が食い違う領域で、BMAが「モデルの不確実性」を正直に反映しているためです
右上（モデル事後確率）: 複数の次数に分散した事後確率が得られています。仮に一つの次数が50%以上の重みを持っていても、残りのモデルに割り当てられた確率は無視できない大きさです。これが単一モデル選択との本質的な違いです
左下（各モデルの予測と重み）: 事後確率の大きいモデルほど太い線で描かれています。訓練データの範囲内では各モデルの予測は似ていますが、外挿領域では大きく異なっています。BMA予測は各モデルの重み付き平均として、バランスの取れた予測を与えています
右下（予測不確実性の分解）: BMAの予測分散が「モデル内不確実性」と「モデル間不確実性」に分解されています。訓練データの近くではモデル内不確実性が支配的ですが、外挿領域ではモデル間不確実性が急激に増大しています。これは、外挿では「どのモデルが正しいか」がさらに不確実になることを反映しています

次に、BMAの予測性能を定量的に評価し、単一モデル選択との比較を行いましょう。

BMAの予測性能の評価

テストデータでの比較

BMAが単一モデル選択よりも実際に優れた予測を行うかどうかを、テストデータで評価します。ここでは、予測の精度（RMSE）と不確実性の信頼性（キャリブレーション）の両方を評価します。

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

def true_function(x):
    return np.sin(2 * x) + 0.3 * x

def polynomial_basis(x, degree):
    return np.column_stack([x**i for i in range(degree + 1)])

def bayesian_lr(x_tr, y_tr, x_te, deg, alpha=0.01):
    """簡易ベイズ線形回帰"""
    Phi_tr = polynomial_basis(x_tr, deg)
    Phi_te = polynomial_basis(x_te, deg)
    N, d = Phi_tr.shape
    w_mle = np.linalg.lstsq(Phi_tr, y_tr, rcond=None)[0]
    res = y_tr - Phi_tr @ w_mle
    beta = 1.0 / max(np.mean(res**2), 1e-6)
    S_inv = alpha * np.eye(d) + beta * Phi_tr.T @ Phi_tr
    S = np.linalg.inv(S_inv)
    m = S @ (beta * Phi_tr.T @ y_tr)
    y_mean = Phi_te @ m
    y_var = np.array([phi @ S @ phi + 1.0/beta for phi in Phi_te])
    log_lik = -N/2*np.log(2*np.pi) - N/2*np.log(np.mean(res**2)) - N/2
    bic = -2 * log_lik + d * np.log(N)
    return y_mean, np.sqrt(y_var), bic

# 多数のデータセットでシミュレーション
np.random.seed(123)
n_simulations = 200
N_train = 20
N_test = 100
sigma_noise = 0.3
degrees = list(range(1, 8))

rmse_bma_list = []
rmse_best_list = []
coverage_bma_list = []
coverage_best_list = []

for sim in range(n_simulations):
    x_tr = np.sort(np.random.uniform(0, 4, N_train))
    y_tr = true_function(x_tr) + np.random.normal(0, sigma_noise, N_train)
    x_te = np.sort(np.random.uniform(-0.5, 5, N_test))
    y_te = true_function(x_te) + np.random.normal(0, sigma_noise, N_test)

    # 各モデルの予測
    means, stds, bics = [], [], []
    for deg in degrees:
        m, s, b = bayesian_lr(x_tr, y_tr, x_te, deg)
        means.append(m)
        stds.append(s)
        bics.append(b)
    bics = np.array(bics)

    # BMA重み
    log_w = -0.5 * bics
    log_w -= logsumexp(log_w)
    weights = np.exp(log_w)

    # BMA予測
    bma_mean = sum(w * m for w, m in zip(weights, means))
    bma_var = sum(w * (s**2 + (m - bma_mean)**2) for w, m, s in zip(weights, means, stds))
    bma_std = np.sqrt(bma_var)

    # 最良モデル
    best_idx = np.argmin(bics)
    best_mean = means[best_idx]
    best_std = stds[best_idx]

    # RMSE
    rmse_bma_list.append(np.sqrt(np.mean((bma_mean - y_te)**2)))
    rmse_best_list.append(np.sqrt(np.mean((best_mean - y_te)**2)))

    # 95%区間のカバレッジ
    cov_bma = np.mean(np.abs(y_te - bma_mean) < 1.96 * bma_std)
    cov_best = np.mean(np.abs(y_te - best_mean) < 1.96 * best_std)
    coverage_bma_list.append(cov_bma)
    coverage_best_list.append(cov_best)

# 結果の可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# RMSE比較
ax = axes[0]
ax.hist(rmse_bma_list, bins=30, alpha=0.6, color='blue', label='BMA', density=True)
ax.hist(rmse_best_list, bins=30, alpha=0.6, color='red', label='最良モデル', density=True)
ax.axvline(np.mean(rmse_bma_list), color='blue', linestyle='--', linewidth=2)
ax.axvline(np.mean(rmse_best_list), color='red', linestyle='--', linewidth=2)
ax.set_xlabel('RMSE')
ax.set_ylabel('密度')
ax.set_title(f'RMSE分布 (200回シミュレーション)\n'
             f'BMA平均: {np.mean(rmse_bma_list):.4f}, '
             f'最良モデル平均: {np.mean(rmse_best_list):.4f}')
ax.legend()

# カバレッジ比較
ax = axes[1]
ax.hist(coverage_bma_list, bins=20, alpha=0.6, color='blue', label='BMA', density=True)
ax.hist(coverage_best_list, bins=20, alpha=0.6, color='red', label='最良モデル', density=True)
ax.axvline(0.95, color='black', linestyle=':', linewidth=2, label='目標 (95%)')
ax.axvline(np.mean(coverage_bma_list), color='blue', linestyle='--', linewidth=2)
ax.axvline(np.mean(coverage_best_list), color='red', linestyle='--', linewidth=2)
ax.set_xlabel('95%区間のカバレッジ')
ax.set_ylabel('密度')
ax.set_title(f'カバレッジ分布\n'
             f'BMA平均: {np.mean(coverage_bma_list):.3f}, '
             f'最良モデル平均: {np.mean(coverage_best_list):.3f}')
ax.legend()

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

# BMAが勝つ割合
win_rate = np.mean(np.array(rmse_bma_list) < np.array(rmse_best_list))
print(f"\nBMAがRMSEで勝つ割合: {win_rate:.1%}")
print(f"BMA平均RMSE: {np.mean(rmse_bma_list):.4f}")
print(f"最良モデル平均RMSE: {np.mean(rmse_best_list):.4f}")

200回のシミュレーション結果から、以下のことが読み取れます。

RMSE: BMAの平均RMSEは単一最良モデルの平均RMSEとほぼ同等か、わずかに良好です。BMAが大幅に劣るケースはほとんどなく、最良モデルが外れた場合にBMAがリスクを軽減しています。これは、BMAが「間違ったモデルに全てを賭ける」リスクを分散するためです
カバレッジ: BMAの95%予測区間のカバレッジは、目標の95%に近い値を示しています。一方、単一モデルのカバレッジはやや低い傾向があります。これは、単一モデルが「モデル間不確実性」を無視するために予測区間が狭すぎることを意味しています。つまり、BMAは不確実性の推定においてよりキャリブレートされた（信頼性の高い）予測を行います

この結果は、BMAの主な利点が「精度の向上」よりも「不確実性推定の改善」にあることを示唆しています。

BMAとモデル選択の理論的比較

AIC/BICとBMAの根本的な違い

AIC（赤池情報量基準）やBIC（ベイズ情報量基準）によるモデル選択と、BMAは一見似ていますが、根本的に異なるアプローチです。

モデル選択（AIC/BIC）の立場: 「真のモデルは候補の中に1つだけ存在する（か、その近似がある）。その最良のモデルを見つけたい」

BMAの立場: 「どのモデルが真であるかは不確実であり、その不確実性を予測に反映すべき」

この違いは、予測の使い方にも影響します。

観点	モデル選択	BMA
目的	1つのベストモデルを選ぶ	全モデルの予測を統合する
出力	単一モデルの予測	混合予測分布
不確実性	パラメータの不確実性のみ	パラメータ + モデルの不確実性
外挿	選択されたモデルに依存	複数モデルの折衷
解釈性	高い（1つのモデルで説明）	やや低い（混合分布）

BMAが有利な場面

BMAは以下のような場面で特に威力を発揮します。

複数のモデルが拮抗している場合: 2つ以上のモデルが同程度にデータをよく説明するとき、どちらか1つを選ぶよりも両方を使う方がロバストです
外挿が必要な場合: モデル間の予測の差は外挿領域で大きくなるため、モデルの不確実性を考慮することが重要です
不確実性の定量化が重要な場合: 医療診断やリスク管理のように、予測の信頼度が意思決定に直結する場面

BMAが不利な場面

一方、以下の場面ではBMAの利点が薄れるか、デメリットが生じます。

1つのモデルが圧倒的に優れている場合: モデル事後確率がほぼ1つのモデルに集中すると、BMAは単一モデル選択とほぼ同じ結果を与えます
解釈性が最優先の場合: 科学的な仮説検証では、「どのモデルが正しいか」の結論が必要なことがあり、混合予測は解釈しにくい場合があります
候補モデルが全て不適切な場合: 全ての候補モデルが真の生成過程からかけ離れている（M-open設定の）場合、BMAも良い予測を与えるとは限りません。この場合は、モデルの改善やベイズ的スタッキングなどの代替手法が有効です

ベイズ的スタッキングとの関係

BMAが全ての候補モデルが「真のモデルの候補」であることを前提としているのに対し、ベイズ的スタッキング（Bayesian stacking）は、交差検証による予測性能を最大化するようにモデルの重みを決定します。具体的には、LOO-CVの対数予測密度を最大化する重みベクトルを最適化問題として解きます。

$$ \hat{\bm{w}}_{\text{stack}} = \arg\max_{\bm{w}} \sum_{i=1}^{N} \log \sum_{k=1}^{K} w_k \, p(y_i | \mathcal{D}_{-i}, M_k) $$

ベイズ的スタッキングはM-open設定（真のモデルが候補に含まれない場合）でより良い予測性能を示すことが知られています。

まとめ

本記事では、ベイズモデル平均化（BMA）の理論と実装について解説しました。

単一モデル選択の限界: 1つのモデルだけを選んで予測することは、モデルの不確実性を無視し、予測区間を過小評価する
BMAの定式化: BMA予測分布 $p(y_{\text{new}} | \mathcal{D}) = \sum_k p(y_{\text{new}} | \mathcal{D}, M_k) p(M_k | \mathcal{D})$ は、パラメータとモデルの両方の不確実性を反映する
全分散の法則: BMAの予測分散は「モデル内分散」と「モデル間分散」に分解でき、後者は単一モデル選択で完全に無視される
BICによる近似: 周辺尤度をBICで近似することで、モデル事後確率を簡便に計算できる
実践的な利点: BMAは予測精度よりも、不確実性推定の信頼性（キャリブレーション）において大きな利点を持つ

BMAは「ベイズ的な考え方」を最も忠実に体現する手法の一つです。パラメータの不確実性だけでなくモデルの不確実性も考慮することで、より信頼性の高い予測を実現します。

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

ベイズ最適化の理論 — ガウス過程と獲得関数でブラックボックス関数を効率的に最適化する — ガウス過程を使ったベイズ最適化
ベイズモデル比較 — WAIC・LOO-CVの理論と実践 — 予測分布に基づくモデル比較の詳細
階層ベイズモデルの理論 — パラメータの階層構造によるモデルの柔軟化