MCMC入門：Metropolis-HastingsとHMCで事後分布からサンプリングする

はじめに

ベイズ推定では、事後分布の形は「尤度 × 事前」に比例することまでは分かります。ところが分母、つまり周辺尤度の積分が計算できません。正規化定数が手に入らないのに、その分布からサンプルを取りたい。これを可能にするのがMCMC（マルコフ連鎖モンテカルロ）です。

MCMCといっても手法によって効率は桁で違います。今回、相関の強い分布を標的に、基本的なランダムウォーク型のMetropolis-HastingsとHMC（ハミルトニアンモンテカルロ）を比べました。同じ3500サンプルを取っても、そのうち「実質的に独立とみなせる本数」（有効サンプルサイズ、ESS）は、素朴な方法では一桁にとどまり、HMCではほぼ3500全部に届きます。同じ手間から引き出せる情報が、手法しだいで2桁違うということです。この記事では、なぜ正規化定数なしにサンプリングできるのかを押さえ、収束診断まで含めて実装で確かめます。

対象読者:

• ベイズ推定を実装したい方
• 重点サンプリングは知っているが、MCMCは手を動かしたことがない方
• 収束診断（トレース・自己相関・R-hat・ESS）の意味を理解して使いたい方

記事のポイント:

• なぜMCMCが正規化定数なしにサンプリングできるのかを理解します
• Metropolis-Hastingsの受容確率と詳細釣り合いを理解します
• 収束診断を使い、MHとHMCの効率差を数値で確認します

---

正規化定数の壁

事後分布は $p(\theta \mid D) \propto p(D \mid \theta)\, p(\theta)$ と書けますが、これを確率分布にするには周辺尤度 $p(D) = \int p(D\mid\theta)p(\theta)\,d\theta$ で割る必要があります。この積分はパラメータの次元が高いと解けません。

MCMCの発想は、正規化されていない密度 $\pi(x)$ を点ごとに評価できさえすればよい、というものです。目標分布を定常分布に持つマルコフ連鎖を作り、それを長く回せば、連鎖の訪れる頻度が目標分布に一致します。受容確率が密度の比 $\pi(x')/\pi(x)$ だけで決まるため、正規化定数は約分されて消えます。ここがMCMCの鍵です。

Metropolis-Hastings

もっとも基本的なのがMetropolis-Hastingsです。現在地 $x$ から提案分布 $q$ で候補 $x'$ を生成し、次の受容確率で受け入れます。

$$\alpha = \min\left(1, \frac{\pi(x')\,q(x \mid x')}{\pi(x)\,q(x' \mid x)}\right)$$

提案が対称（ランダムウォークなら $q(x'\mid x)=q(x\mid x')$）なら、$\alpha = \min(1, \pi(x')/\pi(x))$ と簡単になります。この受容規則は詳細釣り合い（detailed balance）を満たすように作られており、それが $\pi$ を連鎖の定常分布にすることを保証します。実装はごく短いです。

def metropolis(step, n, seed=0):
    rng = np.random.default_rng(seed)
    x, lp = np.zeros(2), logp(np.zeros(2))
    out, acc = np.zeros((n, 2)), 0
    for i in range(n):
        xp = x + rng.normal(0, step, 2)          # 対称な提案
        lpp = logp(xp)
        if np.log(rng.random()) < lpp - lp:      # 受容判定（対数で比較）
            x, lp = xp, lpp; acc += 1
        out[i] = x
    return out, acc / n

刻み幅の調整が悩ましいところです。小さすぎると受容率は高いですが、少しずつしか動けず自己相関が高くなります。大きすぎると提案がことごとく棄却され、やはり進みません。ランダムウォークMHでは受容率0.234あたりが効率の目安とされています。

より効率的な手法

ギブスサンプリングは、条件付き分布から各変数を順に更新する方法で、共役な場合に強いです。もう一段効率を上げるのがHMC（ハミルトニアンモンテカルロ）です。パラメータに運動量 $p$ を導入し、$H(x,p) = -\log\pi(x) + \frac{1}{2}p^\top p$ というハミルトニアンに従ってしばらく運動させてから次の点を提案します。勾配を使って分布の尾根に沿って遠くまで移動できるため、ランダムウォークにありがちなその場足踏みを避けられます。

運動はリープフロッグ法で数値積分します。実装の核はこの数行です。

p_new = p + 0.5 * step * grad_logp(x_new)        # 半ステップ
for l in range(L):
    x_new = x_new + step * p_new                 # 位置を更新
    if l < L - 1:
        p_new += step * grad_logp(x_new)         # 運動量を更新
p_new += 0.5 * step * grad_logp(x_new)           # 残り半ステップ

軌道長を自動調整するNUTSはHMCの発展版で、PyMCやStanが標準採用しています。今回はHMCの素の挙動を見るため、軌道長は固定しました。

相関ガウスでサンプリングする

相関係数0.95の2次元ガウスを標的に、刻み幅0.05のMH、刻み幅0.6のMH、HMCの3つを各4000サンプル（バーンイン500）回します。勾配は $\nabla\log\pi(x) = -\Sigma^{-1}x$ で解析的に与えました。R-hatは分散した4つの初期値から複数チェーンを走らせて計算します。乱数シードは42に固定しました。

まず、調整したMHとHMCが標的をどう探るかを見ます。MHは局所的にうろつくのに対し、HMCは相関の尾根に沿って全体を素早く覆います。

トレースプロットにその差が端的に出ます。刻み幅0.05のMHはゆっくり漂うだけで全域を回り切れません。HMCはほぼ白色雑音のように見え、毎ステップで大きく動いています。

自己相関関数で定量化すると、HMCは10ラグ程度でほぼ0に落ちるのに対し、MHは長く尾を引きます。

収束診断をまとめます。

サンプラー	受容率	lag-1自己相関	ESS(3500中)	R-hat
MH (刻み0.05)	0.931	0.997	8	1.809
MH (刻み0.6)	0.384	0.958	44	1.008
HMC	0.959	-0.742	3500	1.000

診断結果の読み方

受容率が高いことは、質の高さを意味しません。刻み0.05のMHは受容率0.931と一見良さそうですが、ほとんど動いていないから受容されているだけで、自己相関0.997、ESSはわずか8です。R-hatも1.809と、収束の目安1.1を大きく超えます。3500サンプル取っても、独立な情報としては8個分しかありません。受容率だけを見て安心すると、この罠にはまります。

HMCの効率は際立っています。ESSが3500、つまりサンプルがほぼ独立に取れています。lag-1自己相関が負なのは、HMCが尾根を横切って反対側へ跳ぶために起きる現象で、むしろ効率の良さの表れです。勾配情報を使って分布の形に沿って動けることが、強い相関のある標的でここまで差を生みます。強い相関や高次元の事後で、ランダムウォークが行き詰まったらHMCを検討する、というのが実務の判断になります。

収束診断は必ず複数指標で見ます。R-hatは複数チェーンが同じ分布に収束したかを、ESSは実効的な標本数を測ります。トレースプロットで目視、自己相関で混合を、R-hatで収束を確認します。どれか1つでは見落とします。特にR-hatは、分散した初期値から複数チェーンを回さないと意味をなしません。

手放しで使える手法ではありません。MCMCはバーンイン期間の除去が要り、多峰性の分布では山の間を渡れずに1つの山に閉じ込められることがあります（トンネリング問題）。HMCは勾配が必要なので、微分できない目的や離散パラメータには使えません。PyMCやStanを使うにしても、内部でこうした診断が何を見ているかを理解していないと、収束していない結果を鵜呑みにしかねません。

まとめ

MCMCは、正規化定数を回避してベイズ推定を実用に載せる基盤技術です。今回の相関ガウスでは、同じ3500サンプルからの有効サンプルサイズが、刻みの小さいMHで8、HMCで3500と2桁以上違いました。手法の選択と診断が結果を左右します。PyMCやStanを使う場合でも、まずトレース・自己相関・R-hat・ESSを揃えて確認する習慣をつけ、R-hatが1.1を下回り、ESSが目的の推定精度に足りているかを見てから結果を解釈します。