「次数2の多項式と次数9の多項式、どちらのモデルがこのデータをよく説明しているのか？」——機械学習やデータ解析の現場では、こうしたモデル選択の問いに絶えず直面します。素朴に「データへの当てはまりが良い方」を選ぶと、必ずと言っていいほど過適合（オーバーフィッティング）に陥ります。複雑なモデルはどんなデータにもピッタリ合わせられてしまうからです。
この問題に対してベイズ統計が用意する答えが、本記事の主役である 周辺尤度（marginal likelihood）、別名 エビデンス（evidence） です。周辺尤度は、ベイズの定理の分母にひっそりと現れる $p(D)$ という量ですが、これは単なる規格化定数にとどまりません。「そのモデルがデータをどれだけ自然に生成しうるか」を、モデルの複雑さに対するペナルティ込みで測る量なのです。
周辺尤度を理解すると、たとえば次のような場面で強力な武器になります。ひとつは ベイズモデル選択 — 複数の候補モデルから、データを最もよく説明するモデルを過適合を避けつつ選ぶこと。もうひとつは 経験ベイズ法（エビデンス近似） — ハイパーパラメータを周辺尤度の最大化によって最適化すること。ガウス過程回帰のカーネル幅の決定や、ベイズ線形回帰の正則化強度の決定は、まさにこの原理で動いています。さらに統計力学では分配関数、情報理論では符号長と、分野を越えて同じ量が顔を出します。
本記事の内容
- 周辺尤度の直感的な理解（なぜ「分母の積分」が重要なのか）と数学的定義
- ベイズの定理における役割（規格化定数）の導出
- ベイズモデル選択とベイズ因子、そしてオッカムの剃刀効果
- 計算法（解析解・重点サンプリング・ブリッジサンプリング・熱力学積分・変分下界ELBO）
- ベルヌーイ-ベータの解析例とPythonによる実装・検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
周辺尤度の直感 — なぜ「分母の積分」が重要なのか
周辺尤度の定義式をいきなり眺めても、その意味はなかなか伝わってきません。まずはイメージから入りましょう。
あなたが工場の品質管理担当で、ある製造ラインの不良率 $\theta$ を知りたいとします。手元には「データを見る前の見込み」（事前分布 $p(\theta)$）があります。そして実際に製品をいくつか検査して得たデータ $D$ があります。あるパラメータ値 $\theta$ のもとでそのデータが出る確率が尤度 $p(D|\theta)$ です。
ここで素朴な疑問が湧きます。「$\theta$ がいくつなのかは分からないのに、このデータ $D$ そのものは、いったいどれくらい起こりやすいことなのだろう？」。この問いに答えるのが周辺尤度です。$\theta$ の値を一つに決め打ちするのではなく、あらゆる $\theta$ について「その $\theta$ らしさ（事前確率）」で重みづけしながら、尤度を平均するのです。
イメージとしては、無数の専門家がいて、それぞれが「不良率は $\theta$ だ」と異なる主張をしている状況です。各専門家の発言力は事前分布 $p(\theta)$ で決まります。周辺尤度とは、全専門家の予測（尤度）を発言力で加重平均した「総意としてのデータの出やすさ」だと考えられます。
左の図は材料となる事前分布（紫）と尤度（緑、表示用にスケール）です。右の図は両者の積 $p(D|\theta)p(\theta)$ で、このオレンジの曲線が囲む面積こそが周辺尤度 $p(D)$ です。周辺尤度は「尤度と事前が同時に高くなる領域がどれだけあるか」を面積として捉えた量だと読み取れます。事前と尤度の山がうまく重なれば面積は大きく、ずれていれば面積は小さくなります。
この「平均する」という操作が、後でモデル選択における過適合の自動回避につながります。決め打ちした最良の $\theta$ ではなく、$\theta$ の不確かさ全体を考慮するからこそ、複雑すぎるモデルにペナルティがかかるのです。直感がつかめたところで、これを数式できちんと定義しましょう。
ベイズの定理と周辺尤度の定義
周辺尤度に入る前に、ベイズの定理を復習します。得られたデータを $D$（または $\mathcal{D}$）、パラメータのベクトルを $\bm{\theta}$ とすると、ベイズの定理は事後分布 $p(\bm{\theta}|D)$ を次のように与えます。
$$ \begin{equation} p(\bm{\theta}|D) = \frac{p(\bm{\theta})\, p(D|\bm{\theta})}{p(D)} \end{equation} $$
この式を構成する各確率には名前が付いています。左辺は 事後分布、右辺分子の $p(\bm{\theta})$ は 事前分布、$p(D|\bm{\theta})$ は 尤度、そして右辺分母の $p(D)$ が 周辺尤度（エビデンス） です。言葉で書くと次のようになります。
$$ \begin{equation} \text{事後分布} = \frac{\text{事前分布} \times \text{尤度}}{\text{周辺尤度}} \end{equation} $$
さて、この分母 $p(D)$ を具体的に書き下しましょう。$p(D)$ は同時分布 $p(D, \bm{\theta}) = p(D|\bm{\theta})p(\bm{\theta})$ を $\bm{\theta}$ について周辺化（積分消去）したものです。
$$ \begin{equation} p(D) = \int p(D, \bm{\theta})\, d\bm{\theta} = \int p(D|\bm{\theta})\, p(\bm{\theta})\, d\bm{\theta} \end{equation} $$
最初の等号は周辺化（同時分布から一方の変数を積分で消すと、もう一方の周辺分布が得られる）、二番目の等号は条件付き確率の定義 $p(D, \bm{\theta}) = p(D|\bm{\theta})p(\bm{\theta})$ を代入しただけです。(3) 式をよく見ると、尤度 $p(D|\bm{\theta})$ を事前分布 $p(\bm{\theta})$ で重みづけて $\bm{\theta}$ について積分（=周辺化）していることがわかります。これが「周辺尤度」という名前の由来です。
ここで重要な事実を押さえておきましょう。(3) 式は $\bm{\theta}$ で積分消去しているので、結果は $\bm{\theta}$ に依存しません。$p(D)$ はデータ $D$ だけの関数（定数）です。だからこそ、最尤推定やMAP推定で「最適なパラメータ」を一点だけ求めたいときには、$\bm{\theta}$ に依存しない分母 $p(D)$ を無視できるのです。
$$ p(\bm{\theta}|D) \propto p(D|\bm{\theta})\, p(\bm{\theta}) $$
点推定（事後分布の最大値の位置）を求めるだけなら、この比例関係で十分です。では分母 $p(D)$ には何の意味もないのかというと、まったく逆です。次節でその役割を見ます。
周辺尤度の役割① — 規格化定数
事後分布 $p(\bm{\theta}|D)$ は確率分布ですから、全範囲で積分すると1にならなければなりません。
$$ \int p(\bm{\theta}|D)\, d\bm{\theta} = 1 $$
ところが分子 $p(D|\bm{\theta})p(\bm{\theta})$ だけでは、積分しても1にはなりません。これを1にするために割っているのが分母 $p(D)$ です。実際、(1) 式の両辺を $\bm{\theta}$ で積分してみましょう。左辺は1、右辺は
$$ 1 = \int \frac{p(D|\bm{\theta})p(\bm{\theta})}{p(D)}\, d\bm{\theta} = \frac{1}{p(D)}\int p(D|\bm{\theta})p(\bm{\theta})\, d\bm{\theta} $$
ここで $p(D)$ は $\bm{\theta}$ に依存しない定数なので積分の外に出せます。両辺に $p(D)$ を掛ければ、ちょうど (3) 式 $p(D) = \int p(D|\bm{\theta})p(\bm{\theta})d\bm{\theta}$ に戻ります。つまり 周辺尤度とは、分子を確率分布に整える（面積を1にする）ための規格化定数なのです。
左の図は規格化前の分子 $p(D|\theta)p(\theta)$ で、囲む面積は1ではありません。右の図はこれを $p(D)$ で割った事後分布で、面積はちょうど1になっています。$p(D)$ で割るという操作が、形を変えずに高さだけを調整して全体の面積を1に揃えている様子が読み取れます。分子の山の「形」（どの $\theta$ が相対的にもっともらしいか）は規格化の前後で変わらず、変わるのは縦軸のスケールだけです。
規格化定数としての役割だけなら、$p(D)$ は計算の都合上の脇役にすぎません。しかし周辺尤度の真価は、この $p(D)$ という値そのものが「モデルの良さ」を表す指標になる、という点にあります。次節でこの主役としての役割を見ていきましょう。
周辺尤度の役割② — ベイズモデル選択
ここまで、ある一つのモデルの内部での話をしてきました。ここからは視点を上げて、複数のモデルを比べる場面を考えます。
モデルを明示するために記法を拡張します。モデル $M$ のもとでの周辺尤度を $p(D|M)$ と書きます。
$$ \begin{equation} p(D|M) = \int p(D|\bm{\theta}, M)\, p(\bm{\theta}|M)\, d\bm{\theta} \end{equation} $$
これは「モデル $M$ が、自分の事前分布のもとでデータ $D$ をどれだけ生成しやすいか」を表します。さて、モデル自体にも事前確率 $p(M)$ を割り当てれば、データを見た後のモデルの事後確率にベイズの定理を適用できます。
$$ \begin{equation} p(M|D) = \frac{p(D|M)\, p(M)}{p(D)} \end{equation} $$
ここで分母の $p(D) = \sum_M p(D|M)p(M)$ は全モデル共通なので、モデル同士を比較するときには消えます。したがって、モデルの事後確率は
$$ p(M|D) \propto p(D|M)\, p(M) $$
に比例します。もしモデルの事前確率を等しく $p(M_1) = p(M_2)$ と取れば、事後確率の大小はそのまま周辺尤度 $p(D|M)$ の大小で決まります。これがベイズモデル選択の核心です。データへの当てはまり（尤度の最大値）ではなく、周辺尤度（事前で平均した尤度）の大きいモデルを選ぶのです。
左の図は、真の関数が2次のデータに対する1次・2次・7次の多項式フィットです。右の図は次数ごとの対数周辺尤度 $\log p(D|M)$ で、真の次数である2次でちょうどピークを取っています。次数を上げれば訓練データへの当てはまりはいくらでも良くできますが、周辺尤度は2次を超えると下がっていくことが読み取れます。当てはまりの良さだけでは選べない「ちょうど良い複雑さ」を、周辺尤度が自動的に教えてくれているのです。
このように周辺尤度はモデルの絶対評価を与えます。次に、2つのモデルを直接比較する標準的な指標であるベイズ因子を導入しましょう。
ベイズ因子 — モデル比較の物差し
2つのモデル $M_1$ と $M_2$ を比べるとき、最も自然な量はそれぞれの周辺尤度の 比 です。これを ベイズ因子（Bayes factor） と呼びます。
$$ \begin{equation} BF_{12} = \frac{p(D|M_1)}{p(D|M_2)} \end{equation} $$
ベイズ因子の意味は、モデルの事後オッズと事前オッズの関係を見るとはっきりします。(5) 式から2つのモデルの事後確率の比を取ると、共通の分母 $p(D)$ が消えて
$$ \begin{equation} \underbrace{\frac{p(M_1|D)}{p(M_2|D)}}_{\text{事後オッズ}} = \underbrace{\frac{p(D|M_1)}{p(D|M_2)}}_{BF_{12}} \times \underbrace{\frac{p(M_1)}{p(M_2)}}_{\text{事前オッズ}} \end{equation} $$
となります。この式が語っているのは明快です。ベイズ因子とは、データ $D$ を観測したことによって事前オッズが事後オッズへと更新される、その更新の倍率なのです。事前にどちらのモデルも五分五分（事前オッズ1）だと思っていたなら、データを見た後のオッズはそっくりそのままベイズ因子に等しくなります。
ベイズ因子の値をどう解釈するかについては、KassとRafteryによる目安がよく使われます。
この図はベイズ因子の値に応じた「$M_1$ を支持する証拠の強さ」を段階的に示したものです。$BF_{12}$ が1から3程度ならほとんど差はなく、3から20で「それなりに支持」、20から150で「強い支持」、150を超えると「決定的な支持」とされます。重要なのは、これが連続的な強さの尺度であり、頻度論の仮説検定のような「有意/非有意」の二分法ではないことです。データが蓄積されるほどベイズ因子はなめらかに証拠の強さを更新していきます。
ところで、なぜ複雑なモデルは必ずしも有利にならないのでしょうか。当てはまりだけなら複雑な方が常に勝てるはずです。この謎を解く鍵が、次に見る「オッカムの剃刀効果」です。
オッカムの剃刀効果 — 周辺尤度が複雑さを罰する仕組み
「説明は必要以上に複雑にすべきでない」という指針を オッカムの剃刀 と呼びます。驚くべきことに、周辺尤度はこの原理を数学的に内蔵しています。余計なペナルティ項を手で付け加えなくても、$p(D|M)$ を計算するだけで複雑なモデルが自動的に不利になるのです。
そのからくりは、周辺尤度が確率分布であることにあります。周辺尤度 $p(D|M)$ は、起こりうるすべてのデータ $D$ にわたって積分すると1になります。
$$ \int p(D|M)\, dD = 1 $$
これは、各モデルが持つ「確率の総量1」を、説明可能なデータ空間に配分していると見なせます。複雑で柔軟なモデルは、実に多様なデータを説明できます。しかし確率の総量は1で固定されているため、多様なデータに広く薄く確率を配ることになります。一方、単純なモデルは限られた種類のデータしか説明できない代わりに、その狭い範囲に確率を高く集中させます。
この図は横軸に「起こりうるデータ」、縦軸に周辺尤度 $p(D|M)$ を取り、単純・中程度・複雑の3モデルを比べたものです。3本の曲線はいずれも面積が1（規格化されている）ですが、複雑なモデル（赤）は広く薄く、単純なモデル（青）は狭く高い分布になっています。実際に観測されたデータ $D$（破線）の位置で各曲線の高さを比べると、もし観測データがそのモデルの「得意領域」にあるなら、確率を集中させている単純なモデルの方が高い周辺尤度を与えることが読み取れます。
これがオッカムの剃刀効果です。複雑なモデルは「何でも説明できる」がゆえに、特定のデータへの周辺尤度がかえって低くなる。逆に、観測データを過不足なく説明できる「ちょうど良い複雑さ」のモデルが、最も高い周辺尤度を獲得します。先ほどの多項式回帰で2次がピークを取ったのは、まさにこの効果の現れです。
ここまでで周辺尤度の意味と役割は出揃いました。残る課題は「では実際にどうやって $p(D)$ を計算するのか」です。(3) 式の積分は、共役な場合を除いて解析的には解けないことがほとんどです。次節以降で代表的な計算法を順に見ていきます。
計算法① — 解析解（共役モデル）
最も幸運なケースは、事前分布と尤度が 共役 で、積分が閉じた形で解ける場合です。代表例がベルヌーイ尤度とベータ事前の組み合わせです。これは後ほどPython実装でも扱うので、ここで導出しておきます。
不良率 $\theta$ について、事前分布を $\text{Beta}(a_0, b_0)$ とします。
$$ p(\theta) = \frac{\theta^{a_0 – 1}(1-\theta)^{b_0 – 1}}{B(a_0, b_0)} $$
ここで $B(\cdot,\cdot)$ はベータ関数です。$N$ 回の試行で $k$ 回の「表」（不良）が出たとすると、尤度は
$$ p(D|\theta) = \theta^k (1-\theta)^{N-k} $$
です。周辺尤度はこれらの積を $\theta \in [0,1]$ で積分したものです。
$$ \begin{equation} p(D) = \int_0^1 \theta^k (1-\theta)^{N-k} \cdot \frac{\theta^{a_0-1}(1-\theta)^{b_0-1}}{B(a_0, b_0)}\, d\theta \end{equation} $$
指数をまとめると、被積分関数の本体が $\theta^{(a_0+k)-1}(1-\theta)^{(b_0+N-k)-1}$ という形、つまり別のベータ分布の核（規格化前の形）になっていることがわかります。
$$ p(D) = \frac{1}{B(a_0, b_0)} \int_0^1 \theta^{(a_0+k)-1}(1-\theta)^{(b_0+N-k)-1}\, d\theta $$
ベータ関数の定義 $B(\alpha, \beta) = \int_0^1 \theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta$ をそのまま使えば、積分は $B(a_0+k,\, b_0+N-k)$ に等しくなります。よって周辺尤度は
$$ \begin{equation} p(D) = \frac{B(a_0 + k,\; b_0 + N – k)}{B(a_0, b_0)} \end{equation} $$
と閉じた形で求まります。実装上はベータ関数を直接計算すると桁あふれしやすいので、対数ガンマ関数 $\ln\Gamma$（scipy.special.gammaln）を使って対数で評価するのが定石です。ベータ関数は $B(\alpha,\beta)=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$ なので、
$$ \ln p(D) = \big[\ln\Gamma(a_0{+}k) + \ln\Gamma(b_0{+}N{-}k) – \ln\Gamma(a_0{+}b_0{+}N)\big] – \big[\ln\Gamma(a_0) + \ln\Gamma(b_0) – \ln\Gamma(a_0{+}b_0)\big] $$
と書けます。解析解が手に入るのは美しいですが、現実のモデルの多くは共役ではありません。そこで近似計算が必要になります。
計算法② — 重点サンプリング
積分 $p(D) = \int p(D|\bm{\theta})p(\bm{\theta})d\bm{\theta}$ をモンテカルロ法で近似する最も素直な方法は、事前分布から $\bm{\theta}$ をサンプリングして尤度の平均を取ることです。
$$ p(D) = \mathbb{E}_{p(\bm{\theta})}\big[p(D|\bm{\theta})\big] \approx \frac{1}{S}\sum_{s=1}^{S} p(D|\bm{\theta}^{(s)}), \quad \bm{\theta}^{(s)} \sim p(\bm{\theta}) $$
しかしこの素朴な推定は、事前分布が広く尤度のピークが鋭いと、ほとんどのサンプルが尤度のほぼゼロな領域に落ちてしまい、効率が悪くなります。そこでより一般に、扱いやすい 提案分布 $q(\bm{\theta})$ を導入する 重点サンプリング（importance sampling） を使います。
$$ \begin{equation} p(D) = \int p(D|\bm{\theta})p(\bm{\theta})\, d\bm{\theta} = \int \frac{p(D|\bm{\theta})p(\bm{\theta})}{q(\bm{\theta})}\, q(\bm{\theta})\, d\bm{\theta} = \mathbb{E}_{q}\!\left[\frac{p(D|\bm{\theta})p(\bm{\theta})}{q(\bm{\theta})}\right] \end{equation} $$
二番目の等号では分子分母に $q(\bm{\theta})$ を掛けただけです。これにより、$q$ からサンプリングして重み $w^{(s)} = p(D|\bm{\theta}^{(s)})p(\bm{\theta}^{(s)})/q(\bm{\theta}^{(s)})$ の平均を取れば $p(D)$ が推定できます。$q$ を事後分布に近く選ぶほど推定の分散が小さくなります。
この図は、ベルヌーイ-ベータの例で事前分布を提案分布にした重点サンプリング推定が、サンプル数を増やすにつれて解析解（赤の破線）へ収束していく様子です。横軸は対数軸で、初期は推定値が大きく揺れますが、サンプルが数千を超えると解析解にぴったり張り付くことが読み取れます。モンテカルロ法の $1/\sqrt{S}$ の収束則どおり、サンプルを増やすほど推定が安定していきます。
重点サンプリングは手軽ですが、高次元では提案分布の良し悪しに敏感で、推定の分散が爆発しがちです。より頑健に2つの分布の規格化定数比を求める手法として、ブリッジサンプリングがあります。
計算法③ — ブリッジサンプリング
ブリッジサンプリング（bridge sampling）は、規格化定数が既知の提案分布 $g(\bm{\theta})$ と、規格化定数が未知（=求めたい $p(D)$）の事後分布をつなぐ「橋」となる関数 $h(\bm{\theta})$ を導入し、両者の比を安定して推定する方法です。基本となる恒等式は次の形です。
$$ \begin{equation} p(D) = \frac{\mathbb{E}_{g}\big[\tilde p(\bm{\theta})\, h(\bm{\theta})\big]}{\mathbb{E}_{p(\bm{\theta}|D)}\big[g(\bm{\theta})\, h(\bm{\theta})\big]} \end{equation} $$
ここで $\tilde p(\bm{\theta}) = p(D|\bm{\theta})p(\bm{\theta})$ は規格化前の事後（=分子）です。分子は提案分布 $g$ からのサンプルで、分母は事後分布からのサンプル（MCMCで得る）で評価します。橋関数 $h$ を最適に選ぶと、重点サンプリングより格段に分散が小さくなります。
この図は提案分布（紫）と事後分布（青）、そして両者をつなぐ橋関数（オレンジの破線）を描いたものです。重要なのは2つの分布が重なる領域（黄色で強調）の存在です。重点サンプリングは片方の分布だけからサンプリングするのに対し、ブリッジサンプリングは両方からサンプリングして重なり領域を活用するため、2分布のずれが大きくても安定して比を推定できることが読み取れます。橋を架けることで、直接は届かない規格化定数比に手が届くわけです。
ブリッジサンプリングは事後分布と提案分布の重なりに依存します。両者が大きくずれている場合に効果的なのが、温度パラメータを連続的に動かす熱力学積分です。
計算法④ — 熱力学積分
熱力学積分（thermodynamic integration）は、統計力学の自由エネルギー計算からの借用です。事前分布から事後分布へと、逆温度 $t$ を $0$ から $1$ まで連続的に動かす「経路」を考えます。逆温度 $t$ での べき乗事後分布（power posterior） を次で定義します。
$$ p_t(\bm{\theta}|D) \propto p(\bm{\theta})\, p(D|\bm{\theta})^{t} $$
$t=0$ では事前分布、$t=1$ では本来の事後分布になります。各 $t$ での規格化定数を $z(t) = \int p(\bm{\theta})p(D|\bm{\theta})^t d\bm{\theta}$ とすると、$z(0)=1$（事前は規格化済み）、$z(1)=p(D)$ です。$\ln z(t)$ を $t$ で微分すると、期待対数尤度が現れます。
$$ \frac{d}{dt}\ln z(t) = \frac{1}{z(t)}\frac{dz(t)}{dt} = \mathbb{E}_{p_t(\bm{\theta}|D)}\big[\ln p(D|\bm{\theta})\big] $$
二番目の等号は、$z(t)$ の被積分関数を $t$ で微分すると $\ln p(D|\bm{\theta})$ の因子が降りてくることから得られます。両辺を $t=0$ から $1$ まで積分すると、左辺は $\ln z(1) – \ln z(0) = \ln p(D)$ になるので、
$$ \begin{equation} \ln p(D) = \int_0^1 \mathbb{E}_{p_t(\bm{\theta}|D)}\big[\ln p(D|\bm{\theta})\big]\, dt \end{equation} $$
が得られます。これが熱力学積分の公式です。各温度 $t$ で期待対数尤度をMCMCなどで推定し、$t$ について数値積分すれば対数周辺尤度が求まります。
この図は逆温度 $t$（事前の0から事後の1まで）に沿った期待対数尤度の曲線で、この曲線の下の面積がちょうど $\ln p(D)$ になります。$t=0$（事前）では尤度がまだ反映されず低い値、$t=1$（事後）に近づくほどデータに適合して期待対数尤度が上がっていく様子が読み取れます。図のタイトルにあるとおり、数値積分の結果は解析解とよく一致します。温度という1次元のパスに沿って積分するため、高次元でも比較的安定して動くのが長所です。
ここまでの手法はサンプリングに基づくものでした。最後に、最適化に基づく決定論的な近似として、変分推論のELBOを見ます。
計算法⑤ — 変分下界（ELBO）
変分推論（variational inference）では、扱いやすい分布族 $q_\phi(\bm{\theta})$ で事後分布を近似します。このとき対数周辺尤度は次のように分解できます。任意の $q$ に対し、
$$ \begin{equation} \ln p(D) = \underbrace{\mathbb{E}_{q}\!\left[\ln \frac{p(D, \bm{\theta})}{q_\phi(\bm{\theta})}\right]}_{\text{ELBO}\;\mathcal{L}(\phi)} + \underbrace{\mathrm{KL}\big(q_\phi(\bm{\theta}) \,||\, p(\bm{\theta}|D)\big)}_{\geq 0} \end{equation} $$
が成り立ちます。導出は、KLダイバージェンスの定義 $\mathrm{KL}(q||p) = \mathbb{E}_q[\ln q – \ln p(\bm{\theta}|D)]$ に $p(\bm{\theta}|D)=p(D,\bm{\theta})/p(D)$ を代入し、$\ln p(D)$ が $\bm{\theta}$ に依存しないことを使って整理すれば得られます。KLダイバージェンスは常に非負なので、
$$ \ln p(D) \geq \mathcal{L}(\phi) $$
が成り立ちます。この $\mathcal{L}(\phi)$ を ELBO（Evidence Lower BOund、エビデンス下界） と呼びます。名前のとおり、対数周辺尤度（エビデンス）の下界なのです。
この図は変分パラメータ $\phi$ を横軸に取ったELBOの曲線（青）と、一定値である対数周辺尤度 $\ln p(D)$（赤の破線）です。ELBOは常に $\ln p(D)$ 以下にあり、両者の差がちょうど KL$(q||p)$ になっていることが読み取れます。ELBOを最大化することは、この差（KLダイバージェンス）を最小化すること、すなわち近似分布 $q$ を真の事後分布に近づけることと等価です。さらに、最適化されたELBOの値そのものが対数周辺尤度の良い近似（下界）を与えます。これがVAE（変分オートエンコーダ）などの基盤になっている考え方です。
5つの計算法を見てきました。最後に、これらの概念を最も単純なベルヌーイ-ベータモデルで実際にPythonで動かし、過適合との関係まで含めて確認しましょう。
Pythonでの実装と検証
解析解の計算
まず、ベルヌーイ-ベータの周辺尤度を解析解 (9) 式で計算します。対数ガンマ関数で安定に評価するのがポイントです。
import numpy as np
from scipy.special import gammaln
def log_marginal_likelihood_beta_bernoulli(N, k, a0, b0):
"""ベルヌーイ尤度 + ベータ事前 の対数周辺尤度（解析解）
N: 試行回数, k: 表(成功)の回数, a0,b0: 事前 Beta(a0,b0)"""
# log B(a0+k, b0+N-k) - log B(a0, b0)
log_num = gammaln(a0 + k) + gammaln(b0 + N - k) - gammaln(a0 + b0 + N)
log_den = gammaln(a0) + gammaln(b0) - gammaln(a0 + b0)
return log_num - log_den
N, k = 10, 7
a0, b0 = 2.0, 2.0
log_pD = log_marginal_likelihood_beta_bernoulli(N, k, a0, b0)
print(f"対数周辺尤度 log p(D) = {log_pD:.6f}")
print(f"周辺尤度 p(D) = {np.exp(log_pD):.6f}")
このコードを実行すると、$\ln p(D) \approx -6.978$、$p(D) \approx 0.000932$ という値が得られます。これは (9) 式 $p(D)=B(a_0{+}k,\,b_0{+}N{-}k)/B(a_0,b_0)$ をそのまま評価したもので、「事前 $\text{Beta}(2,2)$ のもとで、10回中『ちょうどこの順序で』7回表が出るデータの起こりやすさ」を表す基準値です。以降の近似計算は、すべてこの値に一致するかどうかで検証します。
重点サンプリングとの照合
次に、事前分布を提案分布とした重点サンプリングで $p(D)$ を推定し、解析解と比べます。
import numpy as np
rng = np.random.default_rng(7)
N, k = 10, 7
a0, b0 = 2.0, 2.0
S = 200000
# 事前 Beta(a0,b0) からサンプリングし、尤度を平均
theta = rng.beta(a0, b0, size=S)
lik = theta**k * (1 - theta)**(N - k) # ベルヌーイ尤度（規格化なし）
pD_is = lik.mean() # 重点サンプリング推定（提案=事前）
from scipy.special import gammaln
# 解析解 (9): B(a0+k, b0+N-k) / B(a0, b0)
log_true = (gammaln(a0+k)+gammaln(b0+N-k)-gammaln(a0+b0+N)) \
- (gammaln(a0)+gammaln(b0)-gammaln(a0+b0))
true_pD = np.exp(log_true) # 同じ規格化の解析解
print(f"重点サンプリング推定 p(D) = {pD_is:.6f}")
print(f"解析解 p(D) = {true_pD:.6f}")
print(f"相対誤差 = {abs(pD_is-true_pD)/true_pD*100:.3f} %")
実行すると、重点サンプリング推定値と解析解の相対誤差はおよそ0.04%程度に収まります。提案分布を事前分布に取ったので重み $w=p(D|\theta)p(\theta)/q(\theta)$ が尤度そのものになり、(11) 式が単純な尤度の平均に帰着していることが確認できます。サンプル数 $S$ を増やせば誤差はさらに小さくなり、モンテカルロ近似が解析解に収束することが読み取れます。
モデル選択と過適合のトレードオフ
最後に、本記事の出発点だった「過適合をどう避けるか」を、ベイズ線形回帰の周辺尤度で確認します。多項式の次数を上げながら、訓練誤差と対数周辺尤度がどう変わるかを見ます。
import numpy as np
def log_evidence_linreg(x, y, degree, alpha, beta_noise):
"""ガウス線形回帰の対数周辺尤度（事前 N(0, alpha^-1 I), ノイズ精度 beta_noise）"""
Phi = np.vander(x, degree + 1, increasing=True) # 計画行列
M = Phi.shape[1]
n = len(y)
A = alpha * np.eye(M) + beta_noise * Phi.T @ Phi
mN = beta_noise * np.linalg.solve(A, Phi.T @ y) # 事後平均
resid = y - Phi @ mN
Emn = (beta_noise / 2) * np.sum(resid**2) + (alpha / 2) * mN @ mN
sign, logdetA = np.linalg.slogdet(A)
# Bishop PRML (3.86)
logev = (M/2)*np.log(alpha) + (n/2)*np.log(beta_noise) \
- Emn - 0.5*logdetA - (n/2)*np.log(2*np.pi)
return logev, np.mean(resid**2)
rng = np.random.default_rng(2)
n = 12
x = np.linspace(-1, 1, n)
true = 1.0 + 0.5*x - 1.5*x**2 # 真の関数は2次
sigma = 0.25
y = true + rng.normal(0, sigma, size=n)
alpha, beta_noise = 1.0, 1.0/sigma**2
for d in range(0, 9):
logev, mse = log_evidence_linreg(x, y, d, alpha, beta_noise)
print(f"次数 {d}: 訓練MSE = {mse:.4f}, log p(D|M) = {logev:.3f}")
このコードの出力を見ると、訓練MSEは次数を上げるほど単調に減少し続けます（複雑なモデルほどデータに当てはまる）。ところが対数周辺尤度 $\log p(D|M)$ は、真の次数である2次付近で最大となり、それより複雑にすると減少に転じます。当てはまりの良さと周辺尤度がはっきり別の挙動を示すことが、数値として確認できます。
この図は上のコードの結果を可視化したもので、左軸に訓練誤差（赤）、右軸に対数周辺尤度（オレンジ）を取っています。訓練誤差は次数とともにひたすら下がり続けるのに対し、対数周辺尤度は2次（青い点）でピークを取り、その後ゆるやかに下がっていくことが読み取れます。「訓練誤差の最小化」を信じると際限なく複雑なモデルに引き込まれるのに対し、「周辺尤度の最大化」は自然にちょうど良い複雑さを選んでくれる — オッカムの剃刀効果が、実データで確かに働いていることがわかります。
まとめ
本記事では、周辺尤度（エビデンス）$p(D)$ について、直感から計算法、実装まで解説しました。
- 定義: 周辺尤度は $p(D) = \int p(D|\bm{\theta})p(\bm{\theta})d\bm{\theta}$ で、尤度を事前分布で重みづけて平均した「データの出やすさ」を表す。$\bm{\theta}$ で積分消去されるため $\bm{\theta}$ に依存しない。
- 2つの役割: ベイズの定理の分母として事後分布を1に規格化する定数であると同時に、モデルの良さを測る指標 $p(D|M)$ でもある。
- モデル選択とベイズ因子: 周辺尤度の比 $BF_{12}=p(D|M_1)/p(D|M_2)$ が事前オッズを事後オッズに更新する倍率になる。
- オッカムの剃刀効果: 周辺尤度は確率の総量1を配分するため、複雑なモデルは広く薄く確率を配り、特定のデータへの周辺尤度が下がる。ペナルティを手で付けずとも過適合が自動的に罰される。
- 計算法: 共役なら解析解、それ以外は重点サンプリング・ブリッジサンプリング・熱力学積分・変分下界ELBOで近似する。ELBOは対数周辺尤度の下界 $\ln p(D)\geq\mathcal{L}(\phi)$ を与える。
周辺尤度は、ベイズ統計を「点推定の道具」から「モデルそのものを評価し選ぶ枠組み」へと引き上げる鍵概念です。経験ベイズ法によるハイパーパラメータ最適化、ガウス過程回帰のカーネル選択、VAEの学習など、現代の機械学習の至るところでこの量が中心的な役割を果たしています。
次のステップとして、以下の記事も参考にしてください。