あなたが天気予報を見るとき、「明日の気温は18度です」という点予測と、「明日の気温は15度から21度の間に95%の確率で収まります」という区間予測では、どちらが意思決定に役立つでしょうか。傘を持っていくか、コートを着るかを判断するには、予測の「確信度」がわかることが重要です。
統計モデルを構築してパラメータを推定した後、そのモデルで未知のデータを予測したい場面は頻繁に訪れます。最尤推定（MLE）やMAP推定では、パラメータを一つの値に固定して予測を行います。しかし、データが少ないとき、パラメータの推定には大きな不確実性が伴います。この不確実性を無視して一つのパラメータ値だけで予測すると、過信した予測になってしまうのです。
ベイズ予測分布はこの問題を解決します。パラメータの事後分布を使って、あらゆるパラメータの可能性を「重み付き平均」することで、パラメータの不確実性を織り込んだ予測を実現します。
ベイズ予測分布を理解すると、以下のような応用が開けます。
- 信頼性の高い予測区間: 単なる点予測ではなく、予測の不確実性を定量化できる。金融リスク管理や医療診断で重要
- モデル比較: 予測分布を通じて周辺尤度を計算し、モデルの良さを比較できる（ベイズファクター）
- オンライン学習: データが逐次的に到来する場合、事後分布を更新するたびに予測分布も自動的に更新される
- 能動学習・実験計画: 予測の不確実性が大きい領域を優先的に観測するための指標として利用できる
本記事の内容
- 点予測の限界と予測分布のモチベーション
- ベイズ予測分布の数学的定義と導出
- 正規分布の予測分布（正規-逆ガンマ共役からStudent-t予測分布へ）
- ベイズ線形回帰の予測分布
- Pythonでの実装と可視化（予測区間の描画）
- モデル比較への応用（周辺尤度との関係）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ベイズ推定とは？仕組みについてわかりやすく解説 — 事前分布・事後分布の基本
- ベイズ線形回帰の理論と導出をわかりやすく解説しPythonで実装する — パラメータの事後分布の導出
- 正規分布とは？定義やグラフの見方を徹底解説 — 正規分布の基本性質
点予測の限界 — なぜ予測分布が必要なのか
パラメータの不確実性が予測に影響する
統計モデルにおいて、データからパラメータを推定した後に新しいデータを予測する場面を考えましょう。たとえば、コインを5回投げて3回表が出たとき、次の投げの結果を予測したいとします。
最尤推定では、表の確率 $\theta$ を $\hat{\theta}_{\text{MLE}} = 3/5 = 0.6$ と推定し、次に表が出る確率を $0.6$ と予測します。しかし、たった5回の試行で $\theta = 0.6$ と断定してよいのでしょうか。実際には $\theta$ は $0.3$ かもしれないし、$0.7$ かもしれません。データが少ないため、$\theta$ の推定には大きな不確実性があります。
この問題を直感的に理解するために、2つの状況を比較しましょう。
- 状況A: コインを5回投げて3回表 → $\hat{\theta} = 0.6$
- 状況B: コインを5000回投げて3000回表 → $\hat{\theta} = 0.6$
どちらも点推定は $0.6$ ですが、状況Bの方がはるかに確信度が高いはずです。点予測だけでは、この確信度の違いを表現できません。
点予測が過信を生む具体例
点予測の危険性をもう少し具体的に見てみましょう。パラメータ $\theta$ が不確実なとき、予測分布の形状は点予測から得られるものとは大きく異なります。
たとえば、正規分布 $\mathcal{N}(\mu, \sigma^2)$ からデータをサンプリングする問題で、$\sigma^2$ は既知として $\mu$ だけが未知だとします。$\mu$ の点推定値 $\hat{\mu}$ を使った予測は $\mathcal{N}(\hat{\mu}, \sigma^2)$ となりますが、$\mu$ の不確実性を考慮した予測分布は、$\sigma^2$ よりも大きな分散を持ちます。つまり、パラメータの不確実性を無視すると、予測の裾が実際よりも軽くなり、極端な値の出現確率を過小評価するのです。
これは金融のリスク管理では致命的です。VaR（Value at Risk）を計算するとき、裾の確率を過小評価すると、実際のリスクを低く見積もってしまいます。
このように、点予測には本質的な限界があります。では、パラメータの不確実性を予測に正しく織り込むにはどうすればよいのでしょうか。その答えがベイズ予測分布です。
ベイズ予測分布の数学的定義
直感的な理解 — パラメータを「平均する」
ベイズ予測分布の核心的なアイデアは、実はとても自然です。パラメータ $\theta$ がどの値をとるか不確実なら、あらゆる $\theta$ の可能性を、その確からしさで重み付けて平均するというものです。
日常的なアナロジーで考えてみましょう。明日の通勤時間を予測したいとします。通勤経路が3つあり、それぞれの経路を選ぶ確率が $0.5, 0.3, 0.2$ だとします。各経路での通勤時間の分布がわかっていれば、明日の通勤時間の予測分布は、3つの経路の分布を確率で重み付けた「混合分布」になります。ベイズ予測分布も同じ発想です — ただし、離散的な「経路」の代わりに、連続的な「パラメータ値」について積分するのです。
数学的定義
データ $\mathcal{D} = \{x_1, x_2, \dots, x_N\}$ が観測された下で、新しいデータ点 $x_{\text{new}}$ の予測分布は次のように定義されます。
$$ \begin{equation} p(x_{\text{new}} | \mathcal{D}) = \int p(x_{\text{new}} | \theta) \, p(\theta | \mathcal{D}) \, d\theta \end{equation} $$
この式の各項を解釈しましょう。
- $p(x_{\text{new}} | \theta)$: パラメータ $\theta$ が与えられたときの、新しいデータ $x_{\text{new}}$ の尤度。これはモデルの仮定そのものです
- $p(\theta | \mathcal{D})$: データ $\mathcal{D}$ を観測した後のパラメータの事後分布。ベイズの定理で計算されます
- 積分 $\int \cdot \, d\theta$: パラメータ空間全体にわたって、上の2つの量の積を足し合わせる操作
つまり、予測分布は「パラメータ $\theta$ が $\theta$ である確率」$\times$「$\theta$ の下での予測」を、すべての $\theta$ について積算したものです。これは周辺化（marginalization）と呼ばれる操作で、ベイズ統計の最も重要な道具の一つです。
点予測との比較
点推定（MLE や MAP）では、$\hat{\theta}$ を一つ選んで予測を $p(x_{\text{new}} | \hat{\theta})$ とします。これは、事後分布 $p(\theta | \mathcal{D})$ をデルタ関数 $\delta(\theta – \hat{\theta})$ で近似することに相当します。
$$ p(x_{\text{new}} | \mathcal{D}) \approx \int p(x_{\text{new}} | \theta) \, \delta(\theta – \hat{\theta}) \, d\theta = p(x_{\text{new}} | \hat{\theta}) $$
事後分布が鋭いピーク（データが十分に多い場合）を持つなら、この近似は良好です。しかし、事後分布が広がっている（データが少ない、あるいはパラメータが弱く識別される）場合、この近似は不適切です。予測分布は事後分布の広がり具体的に反映するため、データ量に応じて自動的に不確実性を調整するという優れた性質を持っています。
ベイズ予測分布の定義を理解したところで、次に具体的な分布モデルでこの積分がどのような結果を与えるかを見ていきましょう。まずは最も基本的な正規分布のケースを扱います。
正規分布の予測分布 — 正規-逆ガンマ共役モデル
問題設定
正規分布 $\mathcal{N}(\mu, \sigma^2)$ から $N$ 個のデータ $x_1, \dots, x_N$ が得られたとき、次の観測値 $x_{\text{new}}$ の予測分布を求めましょう。ここでは $\mu$ と $\sigma^2$ の両方が未知とします。
事前分布の設定
正規分布の平均と分散に対する共役事前分布は正規-逆ガンマ分布です。
$$ \begin{equation} \mu | \sigma^2 \sim \mathcal{N}\left(\mu_0, \frac{\sigma^2}{\kappa_0}\right) \end{equation} $$
$$ \begin{equation} \sigma^2 \sim \text{Inv-Gamma}\left(\frac{\nu_0}{2}, \frac{\nu_0 \sigma_0^2}{2}\right) \end{equation} $$
ここで $\mu_0, \kappa_0, \nu_0, \sigma_0^2$ はハイパーパラメータです。$\mu_0$ は $\mu$ の事前的な推測値、$\kappa_0$ は「事前の仮想的なサンプルサイズ」と解釈でき、$\nu_0$ は事前の自由度、$\sigma_0^2$ は分散の事前スケールです。
なぜこの形の事前分布を使うのかというと、$\mu$ の不確実性は $\sigma^2$ に依存するのが自然だからです。$\sigma^2$ が大きい（データのばらつきが大きい）なら、$\mu$ の推定も不確実になります。この依存関係を $\mu | \sigma^2$ の条件付き分布として表現しています。
事後分布の導出
共役性により、事後分布も正規-逆ガンマ分布になります。$N$ 個のデータの標本平均を $\bar{x} = \frac{1}{N}\sum_{i=1}^{N} x_i$、標本分散を $s^2 = \frac{1}{N}\sum_{i=1}^{N} (x_i – \bar{x})^2$ とすると、事後ハイパーパラメータは以下のように更新されます。
$$ \kappa_N = \kappa_0 + N $$
まず、事前の仮想サンプルサイズ $\kappa_0$ に実際のサンプルサイズ $N$ を加えます。
$$ \mu_N = \frac{\kappa_0 \mu_0 + N \bar{x}}{\kappa_N} $$
事後平均 $\mu_N$ は、事前平均 $\mu_0$ と標本平均 $\bar{x}$ の重み付き平均になっています。重みはそれぞれのサンプルサイズ $\kappa_0$ と $N$ です。
$$ \nu_N = \nu_0 + N $$
$$ \nu_N \sigma_N^2 = \nu_0 \sigma_0^2 + N s^2 + \frac{\kappa_0 N}{\kappa_N}(\bar{x} – \mu_0)^2 $$
最後の項 $\frac{\kappa_0 N}{\kappa_N}(\bar{x} – \mu_0)^2$ は、事前平均と標本平均のずれから生じる追加の変動を表しています。
予測分布の導出 — Student-t分布が現れる
ここからが予測分布の核心部分です。新しいデータ $x_{\text{new}}$ の予測分布を求めるために、$\mu$ と $\sigma^2$ を積分消去します。
$$ p(x_{\text{new}} | \mathcal{D}) = \int \int p(x_{\text{new}} | \mu, \sigma^2) \, p(\mu, \sigma^2 | \mathcal{D}) \, d\mu \, d\sigma^2 $$
この二重積分を段階的に実行しましょう。まず、$\mu$ について積分します。$x_{\text{new}} | \mu, \sigma^2 \sim \mathcal{N}(\mu, \sigma^2)$ であり、$\mu | \sigma^2, \mathcal{D} \sim \mathcal{N}(\mu_N, \sigma^2/\kappa_N)$ なので、$\mu$ を積分消去すると、正規分布の畳み込みの公式により、
$$ p(x_{\text{new}} | \sigma^2, \mathcal{D}) = \mathcal{N}\left(x_{\text{new}} \,\middle|\, \mu_N, \, \sigma^2\left(1 + \frac{1}{\kappa_N}\right)\right) $$
となります。分散が $\sigma^2$ から $\sigma^2(1 + 1/\kappa_N)$ に増えていることに注目してください。$1/\kappa_N$ の項は、$\mu$ の不確実性から来る追加の分散です。データが増えると $\kappa_N$ が大きくなり、この追加分散は小さくなります。
次に、$\sigma^2$ について積分します。$\sigma^2 | \mathcal{D} \sim \text{Inv-Gamma}(\nu_N/2, \nu_N\sigma_N^2/2)$ なので、正規分布と逆ガンマ分布の混合を計算する必要があります。この積分を実行すると、結果はStudent-t分布になります。
$$ \begin{equation} p(x_{\text{new}} | \mathcal{D}) = t_{\nu_N}\left(x_{\text{new}} \,\middle|\, \mu_N, \, \sigma_N^2\left(1 + \frac{1}{\kappa_N}\right)\right) \end{equation} $$
ここで $t_{\nu_N}$ は自由度 $\nu_N$ のStudent-t分布です。
Student-t分布が現れる直感的な理由
なぜ正規分布ではなくStudent-t分布が現れるのでしょうか。直感的には、$\sigma^2$ が未知であることが「裾の重さ」をもたらしています。
$\sigma^2$ が既知であれば、$\mu$ だけの不確実性を積分消去しても結果は正規分布です。しかし、$\sigma^2$ が未知のとき、$\sigma^2$ が大きい値をとる可能性もあります。$\sigma^2$ が大きいということは、データのばらつきが大きいことを意味し、極端な値が出やすくなります。この「$\sigma^2$ が大きい可能性」が、予測分布の裾を正規分布よりも重くするのです。
データが増えると自由度 $\nu_N$ が大きくなり、Student-t分布は正規分布に近づきます。これは、データが多いほど $\sigma^2$ の推定が正確になり、$\sigma^2$ の不確実性が減るためです。
予測分布の性質
この予測分布には以下の重要な性質があります。
- 裾が重い: Student-t分布は正規分布よりも裾が重く、極端な値の確率を高く見積もります。これはパラメータの不確実性を正直に反映しています
- データが増えると正規分布に近づく: $N \to \infty$ で $\nu_N \to \infty$ となり、Student-t分布は正規分布に収束します
- 分散が大きい: 点予測の分散 $\hat{\sigma}^2$ よりも、予測分布の分散 $\sigma_N^2(1 + 1/\kappa_N)$ の方が大きくなります。追加の分散は $\mu$ の不確実性を反映しています
正規分布モデルでの予測分布が理解できたところで、次はより実用的なベイズ線形回帰モデルでの予測分布を導出しましょう。
ベイズ線形回帰の予測分布
問題設定
ベイズ線形回帰では、入力 $\bm{x}$ と出力 $y$ の関係を次のようにモデル化します。
$$ \begin{equation} y = \bm{w}^T \bm{\phi}(\bm{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \beta^{-1}) \end{equation} $$
ここで $\bm{\phi}(\bm{x})$ は基底関数ベクトル、$\bm{w}$ は重みパラメータ、$\beta$ は観測ノイズの精度（分散の逆数）です。ノイズの精度 $\beta$ は既知とします。
重みベクトル $\bm{w}$ に正規分布の事前分布 $p(\bm{w}) = \mathcal{N}(\bm{w} | \bm{m}_0, \bm{S}_0)$ を設定し、$N$ 個のデータ $\mathcal{D} = \{(\bm{x}_n, y_n)\}_{n=1}^{N}$ を観測した後の事後分布は、ベイズ線形回帰の記事で導出した通り、
$$ p(\bm{w} | \mathcal{D}) = \mathcal{N}(\bm{w} | \bm{m}_N, \bm{S}_N) $$
です。ここで $\bm{S}_N^{-1} = \bm{S}_0^{-1} + \beta \bm{\Phi}^T \bm{\Phi}$、$\bm{m}_N = \bm{S}_N(\bm{S}_0^{-1}\bm{m}_0 + \beta \bm{\Phi}^T \bm{y})$ です。
予測分布の導出
新しい入力 $\bm{x}_*$ に対する出力 $y_*$ の予測分布を求めます。
$$ p(y_* | \bm{x}_*, \mathcal{D}) = \int p(y_* | \bm{x}_*, \bm{w}) \, p(\bm{w} | \mathcal{D}) \, d\bm{w} $$
被積分関数の各項はどちらも正規分布です。$p(y_* | \bm{x}_*, \bm{w}) = \mathcal{N}(y_* | \bm{w}^T \bm{\phi}(\bm{x}_*), \beta^{-1})$ であり、$p(\bm{w} | \mathcal{D}) = \mathcal{N}(\bm{w} | \bm{m}_N, \bm{S}_N)$ です。
正規分布の線形変換の性質を使います。$y_* = \bm{w}^T \bm{\phi}(\bm{x}_*) + \epsilon$ において、$\bm{w}$ は正規分布に従い、$\epsilon$ も正規分布に従い、両者は独立です。したがって、$y_*$ も正規分布に従います。
$y_*$ の平均は、$\bm{w}$ の平均を代入して、
$$ \mathbb{E}[y_*] = \mathbb{E}[\bm{w}^T \bm{\phi}(\bm{x}_*) + \epsilon] = \bm{m}_N^T \bm{\phi}(\bm{x}_*) $$
$y_*$ の分散は、$\bm{w}$ の不確実性からくる分散とノイズの分散の和です。
$$ \text{Var}[y_*] = \bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*) + \beta^{-1} $$
第1項 $\bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*)$ はパラメータ $\bm{w}$ の不確実性から来る分散であり、認識的不確実性（epistemic uncertainty）に対応します。第2項 $\beta^{-1}$ はデータ自体のノイズから来る分散で、偶然的不確実性（aleatoric uncertainty）に対応します。
したがって、ベイズ線形回帰の予測分布は次のようになります。
$$ \begin{equation} p(y_* | \bm{x}_*, \mathcal{D}) = \mathcal{N}\left(y_* \,\middle|\, \bm{m}_N^T \bm{\phi}(\bm{x}_*), \, \bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*) + \beta^{-1}\right) \end{equation} $$
予測分布の幾何学的解釈
この予測分布には美しい幾何学的構造があります。予測分散 $\sigma_*^2 = \bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*)$ は、入力 $\bm{x}_*$ の位置によって変化します。
- 訓練データの近くでは $\sigma_*^2$ が小さい: パラメータの不確実性が、訓練データによってよく制約されているため
- 訓練データから離れると $\sigma_*^2$ が大きくなる: パラメータの不確実性が予測に大きく影響するため
これは直感に合っています。データがある場所では自信を持って予測でき、データがない場所では「よくわからない」と正直に認めるのです。この性質は、ガウス過程回帰にも共通する特徴です。
予測分布の理論を理解したところで、いよいよPythonでの実装に移りましょう。まずは正規分布モデルでのStudent-t予測分布を実装し、次にベイズ線形回帰の予測分布を可視化します。
Pythonでの実装 — 正規分布モデルの予測分布
Student-t予測分布の実装
まず、正規-逆ガンマ共役モデルにおけるStudent-t予測分布を実装します。少量データの場合と大量データの場合で、予測分布がどのように変化するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 真のパラメータ
true_mu = 5.0
true_sigma = 2.0
# ハイパーパラメータ（弱い事前分布）
mu_0 = 0.0
kappa_0 = 1.0
nu_0 = 1.0
sigma2_0 = 1.0
np.random.seed(42)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
sample_sizes = [5, 30, 200]
for ax, N in zip(axes, sample_sizes):
# データ生成
data = np.random.normal(true_mu, true_sigma, N)
x_bar = np.mean(data)
s2 = np.var(data)
# 事後ハイパーパラメータの更新
kappa_N = kappa_0 + N
mu_N = (kappa_0 * mu_0 + N * x_bar) / kappa_N
nu_N = nu_0 + N
nu_N_sigma2_N = nu_0 * sigma2_0 + N * s2 + (kappa_0 * N / kappa_N) * (x_bar - mu_0)**2
sigma2_N = nu_N_sigma2_N / nu_N
# 予測分布: Student-t
pred_scale = np.sqrt(sigma2_N * (1 + 1 / kappa_N))
pred_dist = stats.t(df=nu_N, loc=mu_N, scale=pred_scale)
# MLE点予測
mle_dist = stats.norm(loc=x_bar, scale=np.sqrt(s2))
# 真の分布
true_dist = stats.norm(loc=true_mu, scale=true_sigma)
# プロット
x = np.linspace(-2, 12, 300)
ax.plot(x, pred_dist.pdf(x), 'b-', linewidth=2, label='予測分布 (Student-t)')
ax.plot(x, mle_dist.pdf(x), 'r--', linewidth=2, label='MLE予測')
ax.plot(x, true_dist.pdf(x), 'k:', linewidth=2, label='真の分布')
ax.fill_between(x, pred_dist.pdf(x), alpha=0.2, color='blue')
ax.set_title(f'N = {N}', fontsize=14)
ax.set_xlabel('x')
ax.set_ylabel('密度')
ax.legend(fontsize=9)
ax.set_ylim(0, 0.35)
plt.tight_layout()
plt.savefig("predictive_distribution_student_t.png", dpi=150, bbox_inches="tight")
plt.show()
上のグラフから、以下の重要な特徴が読み取れます。
- $N = 5$ のとき、予測分布（Student-t）は MLE 予測よりも裾が重い: データが少ないため、パラメータの不確実性が大きく、予測分布はそれを正直に反映して広がっています。MLE 予測は不確実性を過小評価しています
- $N = 30$ のとき、予測分布は MLE 予測に近づくが、まだ差がある: パラメータの推定精度は向上しましたが、完全ではないため、予測分布はわずかに広がっています
- $N = 200$ のとき、両者はほぼ一致する: データが十分多いと、パラメータの事後分布が鋭くなり、点推定と予測分布の差は無視できるほど小さくなります。Student-t 分布の自由度が大きくなり、正規分布に近づいています
この結果は、予測分布が「データ量に応じて自動的に不確実性を調整する」という性質を明確に示しています。
続いて、ベイズ線形回帰の予測分布を実装して、入力空間での予測区間がどのように変化するかを可視化しましょう。
Pythonでの実装 — ベイズ線形回帰の予測分布
予測区間の描画
ベイズ線形回帰の予測分布を実装し、予測平均と95%予測区間を描画します。ここでは多項式基底関数を使います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def polynomial_basis(x, degree):
"""多項式基底関数"""
return np.column_stack([x**i for i in range(degree + 1)])
# 真の関数
def true_function(x):
return 0.5 * x**2 - 2 * x + 1
# データ生成
np.random.seed(42)
N = 15
beta = 4.0 # ノイズ精度 (sigma^2 = 0.25)
x_train = np.random.uniform(-1, 4, N)
y_train = true_function(x_train) + np.random.normal(0, 1/np.sqrt(beta), N)
# モデル設定
degree = 3
Phi = polynomial_basis(x_train, degree)
# 事前分布パラメータ
alpha = 0.1 # 事前精度
m_0 = np.zeros(degree + 1)
S_0_inv = alpha * np.eye(degree + 1)
# 事後分布パラメータ
S_N_inv = S_0_inv + beta * Phi.T @ Phi
S_N = np.linalg.inv(S_N_inv)
m_N = S_N @ (S_0_inv @ m_0 + beta * Phi.T @ y_train)
# 予測
x_test = np.linspace(-1.5, 5, 200)
Phi_test = polynomial_basis(x_test, degree)
# 予測平均と分散
y_mean = Phi_test @ m_N
y_var = np.array([phi.T @ S_N @ phi + 1/beta for phi in Phi_test])
y_std = np.sqrt(y_var)
y_epistemic_std = np.sqrt(np.array([phi.T @ S_N @ phi for phi in Phi_test]))
# プロット
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: 予測分布と信頼区間
ax = axes[0]
ax.plot(x_test, y_mean, 'b-', linewidth=2, label='予測平均')
ax.fill_between(x_test, y_mean - 1.96*y_std, y_mean + 1.96*y_std,
alpha=0.15, color='blue', label='95% 予測区間')
ax.fill_between(x_test, y_mean - 1.96*y_epistemic_std, y_mean + 1.96*y_epistemic_std,
alpha=0.3, color='blue', label='95% 信頼区間 (認識的)')
ax.plot(x_test, true_function(x_test), 'k--', linewidth=1.5, label='真の関数')
ax.scatter(x_train, y_train, c='red', s=40, zorder=5, label='訓練データ')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('ベイズ線形回帰の予測分布')
ax.legend(fontsize=9, loc='upper left')
ax.set_ylim(-4, 10)
# 右: データ数による予測区間の変化
ax = axes[1]
for n_sub, color, ls in [(3, 'red', '--'), (8, 'orange', '-.'), (15, 'blue', '-')]:
Phi_sub = polynomial_basis(x_train[:n_sub], degree)
S_sub_inv = S_0_inv + beta * Phi_sub.T @ Phi_sub
S_sub = np.linalg.inv(S_sub_inv)
m_sub = S_sub @ (S_0_inv @ m_0 + beta * Phi_sub.T @ y_train[:n_sub])
y_m = Phi_test @ m_sub
y_v = np.array([phi.T @ S_sub @ phi + 1/beta for phi in Phi_test])
y_s = np.sqrt(y_v)
ax.plot(x_test, y_m, color=color, linestyle=ls, linewidth=2, label=f'N={n_sub}')
ax.fill_between(x_test, y_m - 1.96*y_s, y_m + 1.96*y_s,
alpha=0.1, color=color)
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, alpha=0.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('データ数と予測区間の関係')
ax.legend(fontsize=9)
ax.set_ylim(-4, 10)
plt.tight_layout()
plt.savefig("bayesian_linear_regression_predictive.png", dpi=150, bbox_inches="tight")
plt.show()
左のグラフからは、予測分布の2種類の不確実性が視覚的に確認できます。
- 内側の帯（認識的不確実性のみ）: パラメータ $\bm{w}$ の不確実性のみを反映した区間です。訓練データの近くでは狭く、データから離れると急速に広がっています。これは「モデルがデータのない領域について自信を持てない」ことを正直に表現しています
- 外側の帯（予測区間）: 認識的不確実性に加えて観測ノイズ $\beta^{-1}$ も含めた区間です。この区間は訓練データの近くでも一定の幅を持ちます。これは、たとえモデルが完璧にパラメータを知っていても、ノイズのために予測は不完全であることを示しています
右のグラフからは、データ数が増えるにつれて予測区間が狭くなる様子が確認できます。$N = 3$ では予測が大きく不確実ですが、$N = 15$ ではかなり真の関数に近い予測が得られています。特に、データが存在する領域では予測区間が大幅に狭まっていますが、外挿領域では依然として広い区間が示されています。
ここまでの実装で予測分布の振る舞いを確認しました。次に、予測分布がモデル比較にどのようにつながるかを理論的に整理しましょう。
事前予測分布と周辺尤度 — モデル比較への応用
事前予測分布と事後予測分布
ここまで扱ってきた予測分布は、データを観測した後の事後予測分布（posterior predictive distribution）です。一方、データを観測する前の事前予測分布（prior predictive distribution）も重要な役割を果たします。
事前予測分布は、事後分布の代わりに事前分布を使って積分したものです。
$$ \begin{equation} p(x_{\text{new}}) = \int p(x_{\text{new}} | \theta) \, p(\theta) \, d\theta \end{equation} $$
この式は、パラメータについての事前知識だけで予測を行ったものです。データを見る前に「このモデルはどのようなデータを生成しそうか」を表しています。
周辺尤度との関係
事前予測分布は、周辺尤度（marginal likelihood）と深い関係があります。$N$ 個のデータ $\mathcal{D} = \{x_1, \dots, x_N\}$ に対する周辺尤度は、
$$ \begin{equation} p(\mathcal{D}) = \int p(\mathcal{D} | \theta) \, p(\theta) \, d\theta \end{equation} $$
です。これはデータ全体の事前予測確率であり、モデル（と事前分布）がデータにどの程度適合するかを測る指標です。
周辺尤度は次のように因数分解できます。
$$ p(\mathcal{D}) = p(x_1) \cdot p(x_2 | x_1) \cdot p(x_3 | x_1, x_2) \cdots p(x_N | x_1, \dots, x_{N-1}) $$
ここで各因子 $p(x_{n+1} | x_1, \dots, x_n)$ は、最初の $n$ 個のデータを見た後の事後予測分布を $x_{n+1}$ で評価したものです。つまり、周辺尤度は逐次的な予測分布の積として解釈できます。
この分解が教えてくれるのは、良いモデルとは「各データ点を、それまでのデータから上手く予測できるモデル」だということです。これが周辺尤度によるモデル比較の直感的な基盤です。
ベイズファクターとモデル選択
2つのモデル $M_1$ と $M_2$ のどちらがデータをよく説明するかを比較するには、周辺尤度の比であるベイズファクターを計算します。
$$ \begin{equation} B_{12} = \frac{p(\mathcal{D} | M_1)}{p(\mathcal{D} | M_2)} \end{equation} $$
$B_{12} > 1$ ならモデル $M_1$ がデータをよりよく説明し、$B_{12} < 1$ なら $M_2$ の方が優れています。
ベイズファクターには自動的なオッカムの剃刀（Bayesian Occam’s razor）という重要な性質があります。複雑なモデルは広い範囲のデータパターンに対応できますが、事前分布を広い空間に分散させるため、特定のデータセットに対する周辺尤度は低くなりがちです。一方、単純なモデルは限られたパターンしか説明できませんが、そのパターンに当てはまるデータに対しては高い周辺尤度を与えます。結果として、データの複雑さに見合ったモデルが自然に選ばれるのです。
Pythonでの周辺尤度計算
正規分布モデルにおける周辺尤度を計算し、異なるモデルを比較してみましょう。
import numpy as np
from scipy import stats
from scipy.special import gammaln
import matplotlib.pyplot as plt
def log_marginal_likelihood_normal(data, mu_0, kappa_0, nu_0, sigma2_0):
"""正規-逆ガンマモデルの対数周辺尤度"""
N = len(data)
x_bar = np.mean(data)
s2 = np.var(data)
kappa_N = kappa_0 + N
nu_N = nu_0 + N
nu_N_sigma2_N = (nu_0 * sigma2_0 + N * s2
+ kappa_0 * N / kappa_N * (x_bar - mu_0)**2)
log_ml = (
- N / 2 * np.log(2 * np.pi)
+ 0.5 * np.log(kappa_0 / kappa_N)
+ gammaln(nu_N / 2) - gammaln(nu_0 / 2)
+ nu_0 / 2 * np.log(nu_0 * sigma2_0)
- nu_N / 2 * np.log(nu_N_sigma2_N)
)
return log_ml
# データ生成: 真の分布 N(3, 1.5^2)
np.random.seed(42)
data = np.random.normal(3.0, 1.5, 30)
# 3つの異なる事前分布でモデル比較
models = {
'M1: mu_0=0, 弱い事前': (0.0, 0.1, 1.0, 10.0),
'M2: mu_0=3, 適切な事前': (3.0, 1.0, 2.0, 2.0),
'M3: mu_0=10, 誤った事前': (10.0, 5.0, 5.0, 1.0),
}
print("=== 対数周辺尤度によるモデル比較 ===\n")
log_mls = {}
for name, (mu_0, kappa_0, nu_0, sigma2_0) in models.items():
log_ml = log_marginal_likelihood_normal(data, mu_0, kappa_0, nu_0, sigma2_0)
log_mls[name] = log_ml
print(f"{name}: log p(D|M) = {log_ml:.2f}")
# ベイズファクター
names = list(log_mls.keys())
print(f"\nベイズファクター B(M2/M1) = {np.exp(log_mls[names[1]] - log_mls[names[0]]):.4f}")
print(f"ベイズファクター B(M2/M3) = {np.exp(log_mls[names[1]] - log_mls[names[2]]):.4f}")
# 事前予測分布と事後予測分布の比較
fig, ax = plt.subplots(figsize=(10, 6))
x_range = np.linspace(-5, 12, 300)
# M2の事前・事後予測分布
mu_0, kappa_0, nu_0, sigma2_0 = models['M2: mu_0=3, 適切な事前']
# 事前予測分布
prior_scale = np.sqrt(sigma2_0 * (1 + 1/kappa_0))
prior_pred = stats.t(df=nu_0, loc=mu_0, scale=prior_scale)
# 事後予測分布
N = len(data)
x_bar, s2 = np.mean(data), np.var(data)
kappa_N = kappa_0 + N
mu_N = (kappa_0 * mu_0 + N * x_bar) / kappa_N
nu_N = nu_0 + N
nu_N_sigma2_N = nu_0 * sigma2_0 + N * s2 + kappa_0 * N / kappa_N * (x_bar - mu_0)**2
sigma2_N = nu_N_sigma2_N / nu_N
post_scale = np.sqrt(sigma2_N * (1 + 1/kappa_N))
post_pred = stats.t(df=nu_N, loc=mu_N, scale=post_scale)
# 真の分布
true_dist = stats.norm(loc=3.0, scale=1.5)
ax.plot(x_range, prior_pred.pdf(x_range), 'g--', linewidth=2, label='事前予測分布')
ax.plot(x_range, post_pred.pdf(x_range), 'b-', linewidth=2, label='事後予測分布')
ax.plot(x_range, true_dist.pdf(x_range), 'k:', linewidth=2, label='真の分布')
ax.hist(data, bins=12, density=True, alpha=0.3, color='gray', label='データ')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('密度', fontsize=12)
ax.set_title('事前予測分布と事後予測分布の比較', fontsize=14)
ax.legend(fontsize=11)
plt.tight_layout()
plt.savefig("prior_posterior_predictive.png", dpi=150, bbox_inches="tight")
plt.show()
この実装から、以下のことが読み取れます。
- 適切な事前分布を持つモデル（M2）が最も高い周辺尤度を示す: 真のパラメータに近い事前知識を持つモデルが、データを最もよく「予測」できています
- 弱い事前分布（M1）は、誤った事前分布（M3）よりも良い: 弱い事前分布は広い範囲のデータに対して中程度の尤度を与えますが、誤った事前分布は的外れな予測をするため周辺尤度が低くなります
- 事前予測分布は広く、事後予測分布は狭い: データを観測することで、予測分布が大幅に絞り込まれている様子が視覚的に確認できます。これがベイズ学習の本質です
逐次予測と予測分布の更新
オンライン学習の枠組み
ベイズ予測分布の大きな利点の一つは、データが逐次的に到来する場合に自然に更新できることです。新しいデータ $x_{N+1}$ が観測されると、事後分布が更新され、それに伴い予測分布も自動的に更新されます。
この過程を逐次的に可視化してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def update_posterior(x_new, mu_n, kappa_n, nu_n, nusigma2_n):
"""1データ点で事後分布を更新"""
kappa_new = kappa_n + 1
mu_new = (kappa_n * mu_n + x_new) / kappa_new
nu_new = nu_n + 1
nusigma2_new = (nusigma2_n
+ kappa_n / kappa_new * (x_new - mu_n)**2)
return mu_new, kappa_new, nu_new, nusigma2_new
# 真のパラメータ
true_mu, true_sigma = 5.0, 1.5
np.random.seed(42)
data = np.random.normal(true_mu, true_sigma, 50)
# 弱い事前分布
mu_n, kappa_n, nu_n = 0.0, 0.01, 0.01
nusigma2_n = 0.01 * 4.0
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
checkpoints = [0, 1, 3, 10, 25, 50]
x_range = np.linspace(-2, 12, 300)
for idx, n_obs in enumerate(checkpoints):
ax = axes[idx // 3, idx % 3]
# n_obs個のデータまで更新
mu_c, kappa_c, nu_c, nusigma2_c = 0.0, 0.01, 0.01, 0.01 * 4.0
for i in range(n_obs):
mu_c, kappa_c, nu_c, nusigma2_c = update_posterior(
data[i], mu_c, kappa_c, nu_c, nusigma2_c)
sigma2_c = nusigma2_c / nu_c if nu_c > 0 else 4.0
pred_scale = np.sqrt(sigma2_c * (1 + 1/kappa_c))
if nu_c > 0:
pred_dist = stats.t(df=nu_c, loc=mu_c, scale=pred_scale)
ax.plot(x_range, pred_dist.pdf(x_range), 'b-', linewidth=2, label='予測分布')
ax.fill_between(x_range, pred_dist.pdf(x_range), alpha=0.2, color='blue')
true_dist = stats.norm(loc=true_mu, scale=true_sigma)
ax.plot(x_range, true_dist.pdf(x_range), 'k--', linewidth=1.5, label='真の分布')
if n_obs > 0:
ax.scatter(data[:n_obs], np.zeros(n_obs) - 0.01, c='red',
s=30, zorder=5, marker='|', linewidth=2)
ax.set_title(f'N = {n_obs}', fontsize=14)
ax.set_xlabel('x')
ax.set_ylabel('密度')
ax.set_ylim(-0.03, 0.5)
ax.legend(fontsize=9)
plt.suptitle('逐次更新による予測分布の変化', fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig("sequential_predictive.png", dpi=150, bbox_inches="tight")
plt.show()
この逐次更新のグラフから、ベイズ予測分布の学習過程が明確に見て取れます。
- $N = 0$（データなし）のとき、予測分布は非常に広い: 事前分布が弱いため、あらゆる値に対して小さな確率密度を割り当てています
- $N = 1$ で早くも分布の位置が動く: たった1つのデータ点でも、予測分布の中心が移動し始めます
- $N = 10$ 程度で分布の形が安定してくる: 予測分布が真の分布に近づいていますが、まだやや広めです。これはパラメータの不確実性が残っていることを反映しています
- $N = 50$ で真の分布にほぼ一致: 十分なデータにより、パラメータの不確実性が小さくなり、予測分布は真の分布とほぼ重なっています
特筆すべきは、ベイズ予測分布が常に真の分布よりもやや広い（分散が大きい）ことです。これはパラメータの不確実性を正直に反映した結果であり、過信した予測を避けるという意味で望ましい性質です。
予測分布の応用 — 異常検知
予測分布の実用的な応用の一つが異常検知（anomaly detection）です。予測分布の下での確率が極端に低いデータ点は、モデルが「予想していなかった」データであり、異常値の候補となります。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 正常データの生成
np.random.seed(42)
normal_data = np.random.normal(5.0, 1.5, 50)
# 異常データの追加
anomalies = np.array([12.0, -3.0, 15.0])
all_data = np.concatenate([normal_data, anomalies])
# 正常データから予測分布を構築
mu_0, kappa_0, nu_0, sigma2_0 = 0.0, 0.01, 0.01, 0.01
x_bar = np.mean(normal_data)
s2 = np.var(normal_data)
N = len(normal_data)
kappa_N = kappa_0 + N
mu_N = (kappa_0 * mu_0 + N * x_bar) / kappa_N
nu_N = nu_0 + N
nu_N_sigma2_N = nu_0 * sigma2_0 + N * s2 + kappa_0 * N / kappa_N * (x_bar - mu_0)**2
sigma2_N = nu_N_sigma2_N / nu_N
pred_scale = np.sqrt(sigma2_N * (1 + 1/kappa_N))
pred_dist = stats.t(df=nu_N, loc=mu_N, scale=pred_scale)
# 各データ点の予測確率
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 予測分布と異常値
ax = axes[0]
x_range = np.linspace(-6, 18, 300)
ax.plot(x_range, pred_dist.pdf(x_range), 'b-', linewidth=2, label='予測分布')
ax.fill_between(x_range, pred_dist.pdf(x_range), alpha=0.2, color='blue')
threshold = pred_dist.ppf(0.01)
threshold_upper = pred_dist.ppf(0.99)
ax.axvline(threshold, color='red', linestyle='--', alpha=0.7, label=f'1%閾値 ({threshold:.1f}, {threshold_upper:.1f})')
ax.axvline(threshold_upper, color='red', linestyle='--', alpha=0.7)
ax.scatter(normal_data, np.full(len(normal_data), 0.005), c='green',
s=30, marker='|', linewidth=2, label='正常データ')
ax.scatter(anomalies, np.full(len(anomalies), 0.005), c='red',
s=100, marker='x', linewidth=3, label='異常値')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('密度', fontsize=12)
ax.set_title('予測分布による異常検知', fontsize=14)
ax.legend(fontsize=10)
# 右: 対数予測確率
ax = axes[1]
log_probs_normal = pred_dist.logpdf(normal_data)
log_probs_anomaly = pred_dist.logpdf(anomalies)
ax.bar(range(len(normal_data)), log_probs_normal, color='green', alpha=0.6, label='正常データ')
ax.bar(range(len(normal_data), len(normal_data) + len(anomalies)),
log_probs_anomaly, color='red', alpha=0.8, label='異常値')
ax.axhline(np.log(pred_dist.pdf(threshold)), color='red', linestyle='--',
alpha=0.7, label='1%閾値')
ax.set_xlabel('データインデックス', fontsize=12)
ax.set_ylabel('対数予測確率', fontsize=12)
ax.set_title('各データ点の対数予測確率', fontsize=14)
ax.legend(fontsize=10)
plt.tight_layout()
plt.savefig("anomaly_detection_predictive.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、予測分布を用いた異常検知の仕組みが理解できます。
- 左のグラフ: 予測分布の99%区間の外側に位置するデータ点（赤い×）が異常値として検出されています。Student-t予測分布は正規分布よりも裾が重いため、正規分布ベースの異常検知よりも保守的に（偽陽性を減らす方向で）動作します
- 右のグラフ: 各データ点の対数予測確率を棒グラフで示しています。正常データの対数確率は比較的高い値に集中していますが、異常値は極端に低い値を示しています。閾値を設定することで、異常値を自動的に検出できます
この方法の利点は、予測分布がパラメータの不確実性を考慮しているため、データが少ない場合でも過度に厳しい閾値を設定しないことです。
まとめ
本記事では、ベイズ予測分布の理論と実装について解説しました。
- 点予測の限界: パラメータを一つの値に固定する点予測は、パラメータの不確実性を無視するため、過信した予測を与えます。特にデータが少ない場合に問題が顕著になります
- ベイズ予測分布の定義: 予測分布 $p(x_{\text{new}} | \mathcal{D}) = \int p(x_{\text{new}} | \theta) p(\theta | \mathcal{D}) d\theta$ は、パラメータの不確実性を積分で織り込んだ予測を実現します
- 正規分布モデルでの予測分布: 平均と分散が未知の場合、予測分布はStudent-t分布になります。正規分布よりも裾が重く、パラメータの不確実性を正直に反映しています
- ベイズ線形回帰の予測分布: 予測分散は認識的不確実性（パラメータの不確実性）と偶然的不確実性（観測ノイズ）の和として分解でき、入力空間の位置に応じて変化します
- 周辺尤度との関係: 予測分布は周辺尤度と密接に関連し、ベイズモデル比較の基盤となります
ベイズ予測分布は、不確実性を正直に伝える予測を行うための理論的に美しく、実用的にも強力な道具です。データが少ない段階では広い予測区間を与え、データが増えるにつれて区間を狭めていくという適応的な振る舞いは、多くの実応用で重要な価値を持ちます。
次のステップとして、以下の記事も参考にしてください。
- ベイズモデル平均化(BMA) — 複数モデルの予測を最適に統合する — 予測分布をさらにモデル間で平均する手法
- ベイズモデル比較 — WAIC・LOO-CVの理論と実践 — 予測分布を用いたモデル評価
- ガウス過程回帰を初めから分かりやすく解説 — 非パラメトリックな予測分布