ハミルトニアンモンテカルロとNUTS

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：ギブスサンプリング　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計）

要点（BLUF）

HMC（ハミルトニアンモンテカルロ）は、事後の勾配と運動量を使い、物理の運動のように「滑空」して遠くへ提案する MCMC。ランダムウォークが強相関・高次元で混合が遅い問題を解消します。
リープフロッグ積分でハミルトン系を数値的に解き、エネルギー保存のずれを受容判定で補正します。
NUTS は軌道長を自動調整する HMC の発展形で、PyMC・Stan の既定サンプラー（API は要最新確認）。実装で、強相関ガウスにおいて HMC の有効サンプルサイズが RW-MH より桁違いに高いことを確かめます。

1. ランダムウォークの限界

メトロポリスヘイスティングスのランダムウォークは、強相関・高次元の事後で少しずつしか動けず、サンプルが互いに強く相関して、独立な情報がなかなか貯まりません。歩幅を大きくすると棄却が増える——進退きわまるのが弱点です。そこで「やみくもに提案する」のをやめ、事後の勾配（どちらが密度の高い方向か）を使って賢く提案します。

2. HMC の考え方：運動量とハミルトニアン

パラメータ $\theta$ に運動量 $p$ （補助変数）を添え、物理系を作ります。位置エネルギーを負の対数事後、運動エネルギーを運動量の二乗で定義します。

U(\theta)=-\log p(\theta\mid D),\qquad K(p)=\tfrac12 p^\top p,\qquad H(\theta,p)=U(\theta)+K(p)

$H$ はハミルトニアン（全エネルギー）。各ステップで運動量を $p\sim\mathcal N(0,I)$ から引き直し、ハミルトンの運動方程式に沿って系を一定時間進めると、エネルギー一定の等高線を滑空して遠くの点に到達します。この軌道を数値的に解くのがリープフロッグ積分（半ステップの運動量更新 → 全ステップの位置更新 → 半ステップの運動量更新）です。

p\leftarrow p+\tfrac{\varepsilon}{2}\nabla\log p(\theta\mid D),\quad \theta\leftarrow\theta+\varepsilon\,p,\quad p\leftarrow p+\tfrac{\varepsilon}{2}\nabla\log p(\theta\mid D)

数値積分の誤差で $H$ がわずかにずれるので、最後に $\min(1,\exp(H_{\text{old}}-H_{\text{new}}))$ で受容判定して補正します（これで正しい事後が定常分布になる）。

3. 実装：HMC は RW-MH より桁違いに効率的

強相関（ $\rho=0.95$ ）の2次元ガウスを目標に、HMC と RW-MH の**有効サンプルサイズ（ESS）**を比べます。ESS は「相関を考慮した実質の独立標本数」（収束診断）。

import numpy as np

rho = 0.95
Sigma = np.array([[1.0, rho], [rho, 1.0]])        # 強相関の目標
Sinv = np.linalg.inv(Sigma)
def log_post(x): return -0.5 * x @ Sinv @ x
def grad_lp(x):  return -Sinv @ x                 # 対数事後の勾配

def ess_1d(x):                                     # 簡易ESS（初期正値までの自己相関和）
    x = x - x.mean(); n = len(x)
    acf = np.correlate(x, x, mode='full')[n-1:]; acf /= acf[0]
    s = 1.0
    for k in range(1, n):
        if acf[k] < 0: break
        s += 2*acf[k]
    return n / s

def hmc(eps=0.18, L=20, N=4000, seed=0):           # リープフロッグ HMC
    rng = np.random.default_rng(seed)
    x = np.zeros(2); out = np.empty((N, 2)); acc = 0
    for i in range(N):
        p = rng.standard_normal(2)
        cur = log_post(x) - 0.5*p@p
        xn, pn = x.copy(), p.copy()
        pn += 0.5*eps*grad_lp(xn)                  # 半ステップ運動量
        for s in range(L):
            xn = xn + eps*pn                        # 全ステップ位置
            if s != L-1: pn += eps*grad_lp(xn)
        pn += 0.5*eps*grad_lp(xn)                  # 半ステップ運動量
        if np.log(rng.uniform()) < (log_post(xn) - 0.5*pn@pn) - cur:
            x = xn; acc += 1
        out[i] = x
    return out, acc/N

def rwmh(step=0.3, N=4000, seed=0):                # 比較用ランダムウォーク MH
    rng = np.random.default_rng(seed)
    x = np.zeros(2); lp = log_post(x); out = np.empty((N, 2)); acc = 0
    for i in range(N):
        cand = x + rng.normal(0, step, 2); lc = log_post(cand)
        if np.log(rng.uniform()) < lc - lp:
            x, lp = cand, lc; acc += 1
        out[i] = x
    return out, acc/N

h, ah = hmc(); m, am = rwmh()
print(f"HMC  : 受容率={ah:.2f}  分散=({h[:,0].var():.2f},{h[:,1].var():.2f})  ESS(x1)={ess_1d(h[:,0]):.0f}/{len(h)}")
print(f"RW-MH: 受容率={am:.2f}  分散=({m[:,0].var():.2f},{m[:,1].var():.2f})  ESS(x1)={ess_1d(m[:,0]):.0f}/{len(m)}")

出力：

HMC  : 受容率=0.95  分散=(1.00,1.00)  ESS(x1)=4000/4000
RW-MH: 受容率=0.60  分散=(1.32,1.30)  ESS(x1)=16/4000

出力の意味：同じ $4000$ ステップでも、HMC の有効サンプルは約 $4000$ （ほぼ独立）、RW-MH はわずか $16$ ——HMC の効率は実に2桁以上上です。RW-MH は強相関の谷に沿ってちまちま動くため分散も $1.32$ と過大（まだ広がりきっていないサイン）ですが、HMC は勾配に沿って一気に対角線方向へ滑空し、分散 $1.00$ ・相関 $0.95$ と目標を正確に再現します。勾配情報を使う効果が、強相関・高次元ほど効いてきます。

4. NUTS と実務

HMC の弱点は軌道長 $L$ とステップ幅 $\varepsilon$ の調整です。 $L$ が短いとランダムウォーク化し、長すぎると軌道が U ターンして無駄足になります。

NUTS（No-U-Turn Sampler）は、軌道が引き返し始めた（U ターンした）ところで自動的に打ち切ることで $L$ の手調整を不要にした HMC の発展形です。実務では自分で HMC を書くより、PyMC・Stan・NumPyro が既定で積む NUTS を使うのが普通です（モデルを宣言的に書けば勾配は自動微分、サンプラーは NUTS が自動調整。具体的な API・既定値は更新が速いので要最新確認。PPL の使い方は確率的プログラミング概観で扱います）。本ノートの自作 HMC は「中で何が起きているか」を見るための最小例です。

flowchart LR
  RW["ランダムウォーク MH<br/>強相関・高次元で混合が遅い"] --> HMC["HMC：勾配＋運動量で滑空<br/>ESS が桁違いに向上"]
  HMC --> T["でも軌道長 L の手調整が必要"]
  T --> NUTS["NUTS：U ターンで自動打ち切り<br/>PyMC/Stan の既定（要最新確認）"]

⚠️ よくある誤解

「HMC は常に最強」ではない：勾配が必要なので、離散パラメータには直接使えません（周辺化やギブスと併用）。勾配計算のコストも要考慮。
「ESS が標本数なら多いほどよい」だけではない：ESS は相関を割り引いた実質独立数。 $4000$ 標本で ESS $16$ なら、実質 $16$ 個ぶんの情報しかない（収束診断）。
「NUTS なら調整不要で必ず収束」ではない：強い相関・漏斗型の事後では NUTS でも発散（divergence）が出ます。再パラメータ化が要る（階層モデルの実例と再パラメータ化）。
「自作 HMC を実務で使う」べきではない：教育目的の最小例です。実務は PyMC/Stan の NUTS（自動微分・自動調整）を使います（要最新確認）。