8.6 MCMCで事後分布からサンプリングする

ベイズ統計では、観測データを得た後の「事後分布」がすべての答えを握っている。パラメータの推定値も、その不確実性も、すべては事後分布の山の形に書かれている。

しかし困ったことに、事後分布は多くの場合、数式で書き下せても直接サンプリングできない。数十次元、数百次元のパラメータが複雑に絡み合った事後分布の山を、どうやって歩き回り、その「形」を知るのか？

ここで登場するのがMCMC — Markov Chain Monte Carlo。特にもっともシンプルな実装であるGibbsサンプリングは、「全部を一気に動かすのが難しいなら、一つずつ動かせばいい」という素朴で力強いアイデアで、事後分布の山を旅する方法を与えてくれる。

このページでは、8.5節で扱ったEMアルゴリズム — 期待値を「計算」する決定論的な手法 — を、期待値を「サンプリング」する確率的な手法へと書き換えると何が起きるのかを見ていく。EMが「最も可能性の高い1点」を指し示すのに対し、Gibbsサンプラーは「事後分布そのものの形」を浮かび上がらせる。

なぜ事後分布のサンプリングが難しいのか

ベイズ的なモデルを立てたとき、最終的に欲しいのは事後分布 Pr(θ | Z) そのものか、その下での期待値だ。たとえば「パラメータの事後平均」が知りたければ、原理的には事後分布から大量にサンプル θ(1), θ(2), …, θ(M) を引いて、その平均を取ればよい。

問題は、その「サンプルを引く」という操作が一般に極めて難しいことだ。事後分布は、尤度と事前分布の積をデータの周辺尤度で割った形をしている：

\Pr(\theta \mid \mathbf{Z}) = \frac{\Pr(\mathbf{Z} \mid \theta)\,\Pr(\theta)}{\int \Pr(\mathbf{Z} \mid \theta)\,\Pr(\theta)\,d\theta}

分母の積分（周辺尤度とか証拠 (evidence) と呼ばれる量）は、パラメータ空間の次元が上がるとほぼ計算不能になる。仮に分母が分からなくても、分子だけなら計算できる — でも、それで一様にサンプルを取る方法がない。この「分子は分かるが、分母（正規化定数）が分からない」という状況を覚えておいてほしい。Section 6で再び現れる。

事後分布を複雑な地形として表現し、直接サンプリングの困難さを示すアニメーション
事後分布を「地形」として見る。複数の峰がある複雑な地形には、地図（正規化定数）がなければランダムに着地させることができない

事後分布を、高次元空間にある「複雑な地形」だと思ってほしい。山あり谷あり、ときには複数の峰が離れて立っている。私たちは地図を持っていない。持っているのは、「いまいる場所の標高は分かる」というセンサーだけ。この限られた情報から、地形全体の様子を知るには、どんな歩き方をすればいいだろうか？

これがMCMCが解こうとしている問題だ。地図がなくても、賢く歩き回ればその場所を訪れる頻度が標高（確率密度）に比例する — そんな歩き方が存在するのである。

\Pr(\theta \mid \mathbf{Z}) \propto \Pr(\mathbf{Z} \mid \theta)\,\Pr(\theta)

この「比例する」という記号（∝）の意味は深い。正規化定数（分母の積分）が分からなくても、比較だけできれば動ける手法の出発点になっている。

Gibbsサンプリング — 「一つずつ動かす」という発想

事後分布が複数のパラメータ U1, U2, …, UK からなるとき、Gibbsサンプリングの発想はこうだ：

「全部を同時に動かすことができないなら、他を固定して、一つずつ順番に動かそう。」

具体的には、次の手続きをひたすら繰り返す。

  1. 適当な初期値 (U1(0), U2(0), …, UK(0)) を置く
  2. 繰り返し t = 1, 2, … について：
    • U1 を更新：U1(t) ～ Pr(U1 | U2(t-1), …, UK(t-1))
    • U2 を更新：U2(t) ～ Pr(U2 | U1(t), U3(t-1), …, UK(t-1))
    • …（以下同様に K まで続く）
U_k^{(t)} \sim \Pr\!\left(U_k \,\Big|\, U_1^{(t)}, \ldots, U_{k-1}^{(t)}, U_{k+1}^{(t-1)}, \ldots, U_K^{(t-1)}, \mathbf{Z}\right)

更新式で上付きの (t) と (t-1) が混在しているのに気づいただろうか？「すでに更新済みの変数は新しい値、まだ更新していない変数は前回の値」を使う。一つ更新するたびに、最新情報を即座に次の更新に反映する逐次更新が肝だ。

何が変わったか？ 一気に高次元の同時分布からサンプルする難しい仕事が、他のすべてを固定したときの条件付き分布から1つだけサンプルする仕事に分解された。この「条件付き分布」は、高次元地形を1本の線で切り取った「スライス」のような確率分布で、多くの場合、正規分布やガンマ分布のようなよく知られた分布になる。

2次元事後分布上でGibbsサンプラーがカクカクと動く様子
2次元等高線上でのGibbsサンプリング。縦方向（U₁固定で U₂ を動かす）と横方向（U₂固定で U₁ を動かす）を交互に繰り返す。軌跡が階段状になる

イメージを掴もう。2次元の事後分布でハイカーが地形を歩く場面を想像してほしい。まず「南北方向にだけ」進める。現在の東西位置を固定したまま、南北の標高分布に従ってジャンプする。次は「東西方向にだけ」進める。今度は南北位置を固定したまま、東西の標高分布に従ってジャンプする。

これを繰り返すと、軌跡は「カクカク」とした階段状になる。不思議なことに、この階段状の軌跡を十分長く続けると、訪問した場所の頻度が、もとの2次元事後分布の標高に比例するようになる。全方位に自由に動かなくても、各軸方向に「正しい確率で」動くだけで十分なのだ。

マルコフ連鎖と定常分布 — なぜ階段歩きで事後分布が分かるのか

Gibbsサンプラーが生成する系列 (U1(t), …, UK(t)) は、マルコフ連鎖だ。つまり、次の状態の確率分布は現在の状態だけで決まり、それ以前の履歴には依存しない。

マルコフ連鎖の偉大な性質は、緩い条件のもとで定常分布に収束することだ。連鎖を十分長く走らせると、ある時点で現在の状態 X(t) の分布が走らせ続けても変わらなくなる。この「動かない分布」が定常分布である。

Gibbsサンプラーで決定的に重要なのは、次の事実だ：

Gibbsサンプラーの定常分布は、ちょうど私たちがサンプルしたかった同時事後分布 Pr(U1, …, UK | Z) である。

なぜか？ 各更新ステップは「条件付き分布」からサンプルしている。条件付き分布と同時分布は、ベイズの公式を通じて整合している。だから、いまの状態がすでに同時事後分布に従っていたら、条件付きサンプルで更新してもその分布は変わらない — これが定常性の本質だ。

多数のGibbsサンプルがプロットされ、ヒストグラムが事後分布に収束していく様子
サンプルが蓄積されるにつれ、ヒストグラムの形が真の事後分布（等高線）と一致していく。バーンイン期間（灰色）の後、収束が始まる

実用上の意味はこうだ：

\hat{\mathbb{E}}[U_k] = \frac{1}{M-m+1} \sum_{t=m}^M U_k^{(t)}

さらに強力なのは、サンプル列から任意の関数 g(U) の事後期待値も直接推定できることだ：

\hat{\mathbb{E}}[g(U) \mid \mathbf{Z}] = \frac{1}{M-m+1} \sum_{t=m}^{M} g\!\left(U^{(t)}\right) \xrightarrow{M \to \infty} \mathbb{E}[g(U) \mid \mathbf{Z}]

事後分散・分位点・信用区間、すべてがサンプルから直接推定できる。「複雑な事後分布の山の形を、サンプルの分布として写し取る」という壮大な仕事が、条件付き分布から順番にサンプルするという単純な操作に還元されたのだ。

ガウス混合モデルへのGibbsサンプラー

ここまで一般論を見てきた。では、具体例で「条件付きで一つずつ」というアイデアが実際にどう動くのかを見てみよう。

8.5節のEMアルゴリズムが解いた問題 — 2成分ガウス混合モデル — を、今度はGibbsサンプラーで解いてみる。設定は同じだ：データ y1, …, yN は、混合比 π で2つの正規分布から生成され、各 yi がどちらから来たかを示す潜在変数 Δi ∈ {0, 1} がある。

EMアルゴリズムでは、Eステップで Δi の期待値（責任度 γi）を計算した。「データ点 yi が第2成分から来ている確率」だ。

\gamma_i(\theta) = \frac{\pi\, \phi_{\theta_2}(y_i)}{(1-\pi)\,\phi_{\theta_1}(y_i) + \pi\,\phi_{\theta_2}(y_i)}

Gibbsサンプラーが行うのは、この期待値を確率として実際にコイン投げをすることだ。各データ点 yi について、確率 γi で Δi = 1、確率 1-γi で Δi = 0 をランダムに決める。

\Delta_i^{(t)} \sim \text{Bernoulli}\!\left(\gamma_i(\theta^{(t-1)})\right)

そして、サンプルされた Δi を使って、今度はパラメータ μ1, μ2 の条件付き分布からサンプルする。第1成分に割り当てられたデータ（Δi = 0 のもの）の平均が μ1 の事後平均となり、その周りの正規分布から μ1 をサンプルする。

ガウス混合のGibbsサンプリング。データ点の所属がランダムに切り替わりながらパラメータが更新される様子
上段：各データ点が確率に応じて青（成分1）か赤（成分2）に色分けされ、2つの正規分布曲線の位置が更新される。下段：μ₁, μ₂ の時系列プロット

Algorithm 8.4 をまとめると：

  1. 初期値 θ(0) = (μ1(0), μ2(0)) を設定
  2. 各反復 t = 1, 2, … について：
    • (a) 潜在変数のサンプル：各 i について、確率 γi(θ(t-1)) で Δi(t) = 1 を、確率 1 - γi(θ(t-1)) で Δi(t) = 0 を生成する（ベルヌーイ試行）
    • (b) パラメータのサンプル：サンプルされた Δ(t) を使って、μ1, μ2 の条件付き分布からそれぞれをサンプルする

これを200回繰り返すと、μ1, μ2, π の事後分布からのサンプル集合が得られる。EMアルゴリズムが「期待値を計算して確定させる」のに対し、Gibbsサンプラーは「確率としてランダムに決める」— このわずかな違いが、出力の性質を根本的に変える。

EMとGibbsの違い — 期待値か、サンプルか

EMアルゴリズムとGibbsサンプラーは、構造的にとてもよく似ている。どちらも「潜在変数を扱うステップ」と「パラメータを更新するステップ」を交互に繰り返す。にもかかわらず、結果が表すものが根本的に違う。

EMアルゴリズムGibbsサンプラー
潜在変数のステップ期待値 γi = E[Δi | θ] を計算同じ γi を確率としてサンプル
パラメータのステップ尤度を最大化する点を決定論的に取る条件付き分布からサンプルする
軌跡の挙動単調に尤度が増え、ピークに収束事後分布をランダムウォーク
出力パラメータの点推定 θ̂MAP事後分布からのサンプル列
同じ尤度面の上で、EM軌跡（単調登山）とGibbs軌跡（ランダムウォーク）を並べて比較するアニメーション
左：EMは単調にピークへ向かう（赤線）。右：Gibbsはピーク周辺をランダムウォークし、分布全体の「形」を写し取る（黄線）

EMが「事後分布の最も高い場所はどこか」を問うのに対し、Gibbsサンプラーは「事後分布の山全体はどんな形か」を問う。

\text{EM}: \quad \hat{\theta} = \arg\max_\theta \Pr(\theta \mid \mathbf{Z})
\text{Gibbs}: \quad \theta^{(1)}, \theta^{(2)}, \ldots, \theta^{(M)} \sim \Pr(\theta \mid \mathbf{Z})

実用上の意味も大きい。EMで得られる答えは1点だけ — その点の周りの不確実性（信用区間、分散など）を知るには別途計算が必要だ。一方Gibbsサンプラーは、サンプル列そのものから事後平均、事後分散、信用区間、任意の関数の事後期待値、すべてが直接得られる。

その代わり、Gibbsサンプラーには計算コストという代償がある。バーンインと十分なサンプル取得のために数千〜数万ステップ走らせる必要があり、各ステップで条件付き分布からのサンプリングが必要だ。EMはずっと速い。

「速さと点推定 vs. 計算と不確実性の全体像」 — このトレードオフが、EMとMCMCを使い分ける鍵になる。

Gibbsの先へ — Metropolis-Hastingsと現代のMCMC

Gibbsサンプラーには、ある重要な前提があった — それは「各変数の条件付き分布から、簡単にサンプリングできる」ということだ。ガウス混合モデルのような綺麗な構造では、条件付き分布が正規分布やベータ分布になり、標準的なサンプラーで扱える。しかし、もっと複雑なモデルでは、条件付き分布が「形は分かるけど、サンプリングは簡単じゃない」ものになりうる。

そのときに登場するのが、より一般的なMetropolis-Hastingsアルゴリズムだ。発想はこうである：

  1. 現在の状態 θ(t) から、適当な「提案分布」 q(θ* | θ(t)) を使って候補 θ* を提案
  2. 受理確率 α を計算する（下の式）
  3. 確率 α で受理（θ(t+1) = θ*）、確率 1-α で棄却（θ(t+1) = θ(t)）
\alpha(\theta^* \mid \theta^{(t)}) = \min\!\left(1, \frac{\Pr(\theta^* \mid \mathbf{Z})\, q(\theta^{(t)} \mid \theta^*)}{\Pr(\theta^{(t)} \mid \mathbf{Z})\, q(\theta^* \mid \theta^{(t)})}\right)

受理確率の式を分解するとこう読める：

要するに「事後密度が高い場所への提案は受け入れやすく、低い場所への提案は確率的に棄却する」という単純なルールだ。

Metropolis-Hastingsの提案と受理/棄却の動作を見せるアニメーション
1次元事後分布の上でのMH法。高い場所への提案は受理（黄色）、低い場所への提案は確率的に棄却（赤でフェード）。ランダムな意思決定が入る

驚くべきことに、Section 1で困った正規化定数（分母の積分）が完全にキャンセルされる — 分子と分母の比を取るからだ。事後密度の比だけ計算できれば、絶対値が分からなくても動く。これがMH法の強力なところだ。

実は、GibbsサンプラーはMetropolis-Hastingsの特殊ケースとして理解できる — 提案分布が条件付き分布で、受理確率が常に1になるケースだ。

現代の機械学習で使われるMCMCは、もっと洗練されている：

しかし、これらすべての出発点にGibbsサンプリングがある。「全部いっぺんが無理なら、条件付きで一つずつ」という素朴で力強い発想が、複雑な事後分布の山を旅する道を切り拓いたのだ。

まとめ

  • MCMCは、ベイズ統計の「事後分布を知る」という核心的な問題に、ランダムウォークで答える手法だ
  • Gibbsサンプリングは、条件付き分布から一つずつ変数を更新する素朴な実装。定常分布として、欲しい同時事後分布が得られる
  • EMとの違いは、「期待値を計算する」か「サンプルする」かの一点だが、出力は点推定 vs. 分布サンプル列という大きな違いを生む
  • Metropolis-HastingsはGibbsを一般化し、条件付き分布が簡単にサンプルできない場合にも対応する