メトロポリス・ヘイスティングス法とは？MCMCの詳細釣り合いから実装まで図解でわかりやすく解説

コインを50回投げて20回表が出た、という観測から「このコインの表が出る確率 $\theta$ はどのあたりか」をベイズ流に推定したいとします。ベイズの定理を書き下すのは簡単です。

$$ p(\theta \mid D) = \frac{p(D \mid \theta)\, p(\theta)}{p(D)} $$

ところが、この式を実際に「使える形」にしようとすると、いきなり壁にぶつかります。分子の $p(D|\theta)p(\theta)$ は、$\theta$ を1つ決めれば電卓でも計算できます。困るのは分母の $p(D) = \int p(D|\theta)p(\theta)\, d\theta$ のほうです。パラメータが1個ならまだしも、階層モデルやニューラルネットのようにパラメータが数十〜数百個あれば、これは数十〜数百重の積分になります。グリッドで刻むのは絶望的です。

ここで問いを立て直してみましょう。分母が計算できないなら、分母を計算しないで済ませる方法はないのか？ 実は、事後分布の「2点での高さの比」だけなら、分母を知らなくても計算できます。$p(\theta_1|D)/p(\theta_2|D)$ を作れば、共通の $p(D)$ は約分されて消えるからです。この当たり前すぎる事実を土台にして、「比の情報だけを頼りに、目的の分布からサンプルを生成する」装置を組み上げたものがマルコフ連鎖モンテカルロ法（MCMC）であり、その中心にあるのがメトロポリス・ヘイスティングス法（Metropolis-Hastings algorithm、MH法）です。

この道具は、ベイズ統計の外にも広く効いています。統計物理学ではイジング模型の平衡状態を調べるのに使われます（そもそもメトロポリス法は1953年に水素爆弾開発の計算物理の現場で生まれました）。ベイズ機械学習では、事後予測分布の計算やベイズニューラルネットの重みのサンプリングに使われます。空間統計、系統樹推定、金融のボラティリティモデル、宇宙論のパラメータ推定 — 「モデルは書けるが正規化定数が計算できない」場面はどこにでもあり、そこには必ずMCMCがいます。

左は1歩の作り方です。現在地から候補を1つ提案し、コイン投げで受け入れるか突き返すかを決める。たったこれだけの操作を何万回も繰り返すと、中央のようにフラフラした軌跡が描かれます。ところがその軌跡の「どこに何回いたか」を数え上げると、右のように目標分布 $\pi(x)$ の形にぴたりと重なります。乱歩の滞在時間分布を設計するのがMCMCだ、と言い換えてもよいでしょう。

本記事の内容

• なぜ独立サンプリングが高次元で破綻し、マルコフ連鎖に頼ることになるのか
• 詳細釣り合い条件の意味と、「詳細釣り合い ⇒ $\pi$ が不変分布」の証明
• MH法の受容確率 $\alpha = \min(1, \text{比})$ の導出と、正規化定数が消える仕組み
• 提案分布のスケール問題と、高次元での最適受容率 0.234
• バーンイン・自己相関・有効サンプルサイズ（ESS）・Gelman-Rubin の $\hat{R}$
• ギブスサンプリングが MH の特殊ケースである理由

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• ベイズの定理とは？わかりやすく解説 — 事後分布と正規化定数の位置づけ
• ベイズ推定とは？仕組みについてわかりやすく解説 — 事前分布から事後分布を得る流れ
• マルコフ連鎖をわかりやすく解説 — 遷移確率と連鎖の基本
• 定常分布と遷移行列 — 本記事の「不変分布」の前提
• モンテカルロ法と重点サンプリング — 独立サンプリング側の代表手法
• 大数の法則をわかりやすく解説 — 標本平均が期待値に収束する根拠

MCMCが解こうとしている問題

欲しいのは積分、手元にあるのは未正規化の密度

ベイズ推論で本当に欲しいものは、事後分布そのものというより、事後分布のもとでの期待値です。事後平均 $\mathbb{E}[\theta|D]$、事後分散、予測分布 $\int p(\tilde{y}|\theta)p(\theta|D)d\theta$、95%信用区間 — どれも積分です。そして積分は、標本さえあればモンテカルロ近似で置き換えられます。

$$ \mathbb{E}_{\pi}[f(X)] \approx \frac{1}{N}\sum_{t=1}^{N} f(x_t), \qquad x_t \sim \pi $$

つまり問題は「積分をどう解くか」ではなく、「目標分布 $\pi$ からどうやって標本を引くか」に翻訳できます。ここで $\pi$ は事後分布 $p(\theta|D)$ のことです。

ところが、その $\pi$ の正体を私たちは完全には知りません。知っているのは正規化されていない形

$$ \tilde{\pi}(\theta) = p(D \mid \theta)\, p(\theta), \qquad \pi(\theta) = \frac{\tilde{\pi}(\theta)}{Z}, \quad Z = p(D) = \int \tilde{\pi}(\theta)\, d\theta $$

だけです。$\tilde{\pi}$ は $\theta$ を代入すれば即座に値が出ます。$Z$ だけが手に負えません。

左が問題の所在です。壊れているのは分母だけで、分子はいつでも評価できます。右がこれから作る仕掛けで、アルゴリズムの中に $\pi$ が「2点での比」の形でしか登場しなければ、$Z$ は分子と分母に同じだけ現れて約分で消えます。MCMCが成立する理由の8割はこの1行に集約されていると言ってよいくらい重要な観察です。

独立サンプリングが高次元で破綻する理由

「比しか分からない分布から標本を引く」古典的な方法として、棄却サンプリング（rejection sampling）があります。目標を包む簡単な分布から候補を引き、目標の高さに応じて受け入れる、という方法です。1次元なら申し分なく動きます。問題は次元です。

包む分布と目標分布のズレは、次元が増えるほど指数的に効いてきます。極端に単純化した例で確かめてみましょう。一辺 $2$ の $d$ 次元立方体から一様に点を打ち、半径 $1$ の球の中に入った点だけを採用する、という「球からのサンプリング」を考えます。受容率は体積比そのもので、

$$ \frac{V_d}{2^d} = \frac{\pi^{d/2}}{\Gamma(d/2+1)\, 2^d} $$

です。この値を実際に計算すると、$d=2$ で 0.785、$d=5$ で 0.164、$d=10$ で $2.49 \times 10^{-3}$、$d=20$ では $2.46 \times 10^{-8}$ まで落ちます。20万点を一様に打つ数値実験でも、$d=10$ では 0.252% しか採用されず、$d=20$ では1点も球に入りませんでした。$d=50$ にいたっては受容率 $1.5 \times 10^{-28}$ で、宇宙の年齢をかけても1点も採れません。

これは「包み方が下手」なのではなく、高次元ではどんな素朴な包み方をしても、包む分布と目標分布の重なりが指数的にゼロへ向かうという構造的な事情です。重点サンプリングでも事情は同じで、重みが一部の標本に集中して有効サンプル数が崩壊します（詳しくはモンテカルロ法と重点サンプリングを参照）。

発想の転換 — 独立性をあきらめる

破綻の原因をたどると、「毎回ゼロから独立に候補を引き直している」ことに行き当たります。前回どこにいたかを完全に忘れて、広い空間に賭け直しているから当たらないのです。

そこで発想を変えます。独立性をあきらめる。前回いた場所の近くから次の候補を出せば、少なくとも「密度がそこそこ高い領域」から出発できます。もちろん、そうやって作った標本列 $x_1, x_2, \dots$ は互いに独立ではありません。しかし独立でなくても、標本平均が正しい期待値に収束してくれるなら、モンテカルロ近似としては十分に役に立ちます。

こうして目指すものが定まりました。目標分布 $\pi$ を不変分布に持つマルコフ連鎖を設計し、それを長く走らせる。 次の節では、そのために連鎖が満たすべき条件を確認します。

マルコフ連鎖が目標分布を「覚える」条件

不変分布 — 一度たどり着いたら出ていかない分布

マルコフ連鎖は、遷移核 $P(x \to x’)$（現在 $x$ にいるとき次に $x’$ へ移る確率密度）で決まります。ある分布 $\pi$ が不変分布（invariant distribution）であるとは、$\pi$ に従って現在地を配ってから1歩進めても、分布が $\pi$ のまま変わらないことを言います。式で書けば

$$ \pi(x’) = \int \pi(x)\, P(x \to x’)\, dx $$

です。左辺は「1歩後に $x’$ にいる確率」、右辺は「どこかにいた確率 × そこから $x’$ に移る確率」を全経路について足したもの。これが釣り合っていれば、連鎖は $\pi$ を保ったまま回り続けます。イメージとしては、水槽の中で水が対流していても水位が変わらない状態です。

既約性と非周期性 — どこへでも行けて、周期に囚われない

不変分布があるだけでは足りません。連鎖が変な場所に閉じこもってしまっては意味がないので、次の2つを課します。

• 既約性（irreducibility）: 状態空間のどこからでも、有限のステップで（正の確率で）どこへでも到達できる。目標分布が2つの山に分かれていても、山と山を行き来できる必要があります。
• 非周期性（aperiodicity）: 「必ず偶数ステップでしか戻ってこない」といった周期構造を持たない。MH法は棄却によって同じ場所に留まる可能性があるので、この条件は自動的に満たされます。棄却が非周期性を保証してくれる、というのは少し面白い副作用です。

エルゴード定理 — 時間平均が空間平均に化ける

不変分布 $\pi$ を持ち、既約かつ非周期な連鎖では、エルゴード定理が成り立ちます。

$$ \frac{1}{N}\sum_{t=1}^{N} f(x_t) \xrightarrow[N \to \infty]{} \mathbb{E}_{\pi}[f(X)] = \int f(x)\, \pi(x)\, dx $$

左辺は1本の連鎖に沿った時間平均、右辺は分布に関する空間平均です。標本が独立でなくても、連鎖を長く走らせれば時間平均が期待値に収束する — これが、独立性をあきらめても構わない理論的な根拠です。独立標本に対する大数の法則の、相関がある場合への一般化だと思えばよいでしょう。

さて、条件は分かりました。しかし「$\pi$ を不変分布に持つ遷移核を作れ」と言われても、上の積分方程式をにらんで直接設計するのは大変です。次の節で、はるかに扱いやすい十分条件を導入します。

詳細釣り合い条件

2つの状態のあいだで、行きと帰りの流れを釣り合わせる

不変分布の条件は「全体として辻褄が合っていればよい」という大域的な要求です。これを、もっと厳しいが確認しやすい局所的な要求に置き換えます。それが詳細釣り合い条件（detailed balance condition）です。

$$ \begin{equation} \pi(x)\, P(x \to x’) = \pi(x’)\, P(x’ \to x) \qquad \text{（すべての } x, x’ \text{ で）} \end{equation} $$

（日本語の教科書では「詳細つり合い条件」と平仮名で書かれることも多く、同じものを指します。）

左辺は「状態 $x$ に $\pi(x)$ だけの確率が溜まっていて、そのうち $P(x\to x’)$ の割合が $x’$ へ流れ出す量」、つまり $x \to x’$ の確率の流量です。右辺は逆向きの流量。詳細釣り合いとは、どの2状態を取り出しても、行きと帰りの流量が完全に打ち消し合うという要求です。

左では $0.6 \times 0.5 = 0.30$ と $0.3 \times 1.0 = 0.30$ が一致し、正味の流れがゼロになっています。右では行きが 0.48、帰りが 0.15 で、差し引き 0.33 が右向きに流れ続けます。この状態を放置すると $x’$ に確率が溜まり続け、分布は $\pi$ から離れていきます。詳細釣り合いが「分布を動かさないための条件」だという直感が、この図から掴めるはずです。

詳細釣り合いならば $\pi$ は不変分布である

では、なぜ詳細釣り合いを満たせば不変分布が保証されるのでしょうか。証明はごく短いので、1行ずつ追ってみましょう。目標は、不変分布の定義式 $\pi(x’) = \int \pi(x) P(x\to x’)\, dx$ を示すことです。

出発点は詳細釣り合いの式 (1) です。この両辺を、$x$ について全空間で積分します。$x’$ は固定したままにしておくのがポイントです。

$$ \int \pi(x)\, P(x \to x’)\, dx = \int \pi(x’)\, P(x’ \to x)\, dx $$

右辺の $\pi(x’)$ は積分変数 $x$ を含まない定数なので、積分の外に出せます。

$$ \int \pi(x)\, P(x \to x’)\, dx = \pi(x’) \int P(x’ \to x)\, dx $$

ここで、残った積分 $\int P(x’ \to x)\, dx$ に注目します。これは「状態 $x’$ から1歩進んだとき、どこかしらに到達する確率」の総和なので、確率の規格化条件からちょうど $1$ です。したがって

$$ \int \pi(x)\, P(x \to x’)\, dx = \pi(x’) \cdot 1 = \pi(x’) $$

となり、不変分布の定義式そのものが得られました。証明完了です。

大事なのは、この含意が一方通行だということです。詳細釣り合いは不変分布であるための十分条件であって、必要条件ではありません（詳細釣り合いを満たさないのに正しい不変分布を持つMCMCも存在します。ギブスサンプリングを決まった順序で回す方式や、Hamiltonian Monte Carlo の一部の変種がそれにあたります）。それでも詳細釣り合いが使われ続けるのは、設計が圧倒的に楽だからです。大域的な積分方程式を解く代わりに、2点ずつの単純な等式を満たせばよいのですから。

この十分条件を武器に、いよいよ具体的な遷移核を組み立てます。

メトロポリス・ヘイスティングス法の受容確率

遷移を「提案」と「受容判定」に分解する

遷移核をゼロから設計するのは難しいので、2段構えにします。

1. 提案分布 $q(x’|x)$ から候補 $x’$ を引く。$q$ は自分で好きに選んでよい、扱いやすい分布（たとえば $\mathcal{N}(x, s^2)$）。
2. 確率 $\alpha(x \to x’)$ で候補を受容し、$1-\alpha$ で棄却して現在地に留まる。

この $\alpha$ は「受容確率」のほか「受理確率」と訳されることもあります（acceptance probability）。本記事では受容で統一しますが、文献によって表記が揺れる点は頭の片隅に置いておいてください。

このとき、$x’ \neq x$ に対する遷移核は2つの因子の積になります。

$$ P(x \to x’) = q(x’ \mid x)\, \alpha(x \to x’) $$

自由に選べる $q$ はいったん脇に置き、詳細釣り合いを満たすように $\alpha$ を決める、というのが戦略です。

受容確率の導出

上の分解を詳細釣り合いの式 (1) に代入します。

$$ \pi(x)\, q(x’ \mid x)\, \alpha(x \to x’) = \pi(x’)\, q(x \mid x’)\, \alpha(x’ \to x) $$

$\alpha$ について整理したいので、両辺を $\pi(x)q(x’|x)\,\alpha(x’ \to x)$ で割ります（いずれも正であるとします）。すると、2つの受容確率の比についての条件式が得られます。

$$ \begin{equation} \frac{\alpha(x \to x’)}{\alpha(x’ \to x)} = \frac{\pi(x’)\, q(x \mid x’)}{\pi(x)\, q(x’ \mid x)} \;\equiv\; r(x, x’) \end{equation} $$

右辺の $r(x,x’)$ をMH比と呼ぶことにします。ここで気づいてほしいのは、式 (2) は $\alpha$ を一意には決めていないということです。比さえ合っていればよいので、たとえば両方を1000分の1にしても条件は満たされます。しかし $\alpha$ を小さくすると棄却ばかりになって連鎖が動かなくなるので、条件を満たす範囲でできるだけ大きい $\alpha$ を選びたい。確率である以上 $\alpha \le 1$ という上限があるので、この制約下で比 (2) を満たす最大の選び方を探すと、答えは次の形になります。

$$ \begin{equation} \alpha(x \to x’) = \min\left(1, \; \frac{\pi(x’)\, q(x \mid x’)}{\pi(x)\, q(x’ \mid x)}\right) = \min\big(1, \, r(x,x’)\big) \end{equation} $$

これがメトロポリス・ヘイスティングスの受容確率です。本当に条件 (2) を満たしているか、$r$ の大小で場合分けして確かめておきましょう。逆向きのMH比が $r(x’,x) = 1/r(x,x’)$ である点を使います。

$r \le 1$ のときは $\alpha(x\to x’) = r$、逆向きは $\alpha(x’\to x) = \min(1, 1/r) = 1$ なので、比は $r/1 = r$。$r > 1$ のときは $\alpha(x\to x’) = 1$、逆向きは $\alpha(x’\to x) = 1/r$ なので、比は $1/(1/r) = r$。どちらの場合も式 (2) が成立しています。

正規化定数が消える瞬間

ここが本記事の山場です。式 (3) に $\pi = \tilde{\pi}/Z$ を代入してみます。

$$ r(x,x’) = \frac{\pi(x’)\, q(x \mid x’)}{\pi(x)\, q(x’ \mid x)} = \frac{\big(\tilde{\pi}(x’)/Z\big)\, q(x \mid x’)}{\big(\tilde{\pi}(x)/Z\big)\, q(x’ \mid x)} = \frac{\tilde{\pi}(x’)\, q(x \mid x’)}{\tilde{\pi}(x)\, q(x’ \mid x)} $$

分子と分母に同じ $1/Z$ が現れて、きれいに約分されました。つまり、アルゴリズムを回すのに $Z = p(D)$ の値は一度も必要ありません。冒頭で立てた「分母を計算しないで済ませられないか」という問いに、これで答えが出ました。尤度と事前分布さえ書ければ、事後分布から標本が引ける。ベイズ統計が実務で使える道具になった決定的な理由が、この約分ひとつなのです。

なお、実装では $\tilde{\pi}$ をそのまま計算すると簡単に桁あふれするので、対数で扱うのが鉄則です。$\log r = \log\tilde{\pi}(x’) – \log\tilde{\pi}(x) + \log q(x|x’) – \log q(x’|x)$ を計算し、$\log u < \log r$（$u \sim \mathrm{Uniform}(0,1)$）で判定します。

メトロポリス法 — 提案が対称な場合

提案分布が対称、すなわち $q(x’|x) = q(x|x’)$ を満たす場合を考えます。現在地を中心とした正規分布 $\mathcal{N}(x, s^2)$ や一様分布 $\mathrm{Uniform}(x-s, x+s)$ はこの性質を持ちます。このとき $q$ の因子が約分されて、受容確率は劇的に単純になります。

$$ \alpha(x \to x’) = \min\left(1, \frac{\pi(x’)}{\pi(x)}\right) $$

これが1953年の原型、メトロポリス法です（$q$ を非対称にまで一般化したのが1970年のヘイスティングスの仕事で、MH法という名前はこの二人に由来します）。この形は言葉にすると非常に分かりやすくなります。

• $\pi(x’) \ge \pi(x)$ のとき（上り坂）: $\alpha = 1$ で、必ず受容する。
• $\pi(x’) < \pi(x)$ のとき（下り坂）: 高さの比 $\pi(x’)/\pi(x)$ の確率で受容する。

左は $x=-3.5$ から $x’=-2.2$ への提案で、密度比は 5.62 なので $\alpha = \min(1, 5.62) = 1.00$、必ず受け入れられます。右は山頂 $x=2.5$ から裾 $x’=5.2$ への提案で、密度比は 0.05 しかないため受容率は約5%、つまり100回提案しても95回は突き返されて同じ場所に留まります。「下り坂も確率的には受け入れる」ことがMCMCの生命線です。上り坂しか登らなければ単なる山登り法になり、モードに貼りついて分布を再現できません。低い場所も低いなりの頻度で訪れるからこそ、訪問頻度が $\pi$ に比例するのです。

アルゴリズムと最小実装

ここまでを手順としてまとめます。

1. 初期値 $x_0$ を決める。
2. $t = 0, 1, 2, \dots$ について繰り返す。
3. 提案 $x’ \sim q(\cdot \mid x_t)$ を引く。
4. MH比 $r$ を（対数で）計算する。
5. $u \sim \mathrm{Uniform}(0,1)$ を引き、$u < \min(1, r)$ なら $x_{t+1} = x'$、そうでなければ $x_{t+1} = x_t$。

コードにすると、本体は20行ほどにしかなりません。冒頭のコイン投げの例（50回中20回成功、事前分布 Beta(2,2)）で動かしてみます。この設定では事後分布が Beta(22, 32) と解析的に分かるので、答え合わせができます。

import numpy as np

# 観測: 50回中20回成功。事前分布は Beta(2, 2)
N_TRIAL, N_SUCC, A0, B0 = 50, 20, 2.0, 2.0

def log_unnorm_posterior(theta):
    """未正規化事後 log[p(D|θ)p(θ)]。正規化定数 p(D) は含まない。"""
    t = float(theta[0])
    if t <= 0.0 or t >= 1.0:
        return -np.inf                                   # 台の外は密度ゼロ
    return (N_SUCC + A0 - 1) * np.log(t) + (N_TRIAL - N_SUCC + B0 - 1) * np.log(1 - t)

def metropolis(log_target, x0, n_iter, scale, rng):
    """対称な正規提案によるランダムウォーク・メトロポリス法"""
    x = np.atleast_1d(np.asarray(x0, dtype=float))
    lp = log_target(x)
    chain = np.empty((n_iter, x.size))
    n_acc = 0
    for t in range(n_iter):
        prop = x + scale * rng.standard_normal(x.size)   # 提案（対称なので q は約分で消える）
        lp_prop = log_target(prop)
        if np.log(rng.random()) < lp_prop - lp:          # u < min(1, r) と同値
            x, lp = prop, lp_prop
            n_acc += 1
        chain[t] = x                                     # 棄却時は同じ値をもう一度記録する
    return chain[:, 0], n_acc / n_iter

rng = np.random.default_rng(22)
chain, acc = metropolis(log_unnorm_posterior, 0.5, 200_000, 0.18, rng)
samples = chain[5000:]                                   # バーンイン5000反復を破棄

print(f"受容率           : {acc:.3f}")
print(f"事後平均   推定  : {samples.mean():.4f}   解析解: {(A0 + N_SUCC) / (A0 + B0 + N_TRIAL):.4f}")
print(f"事後標準偏差 推定: {samples.std(ddof=1):.4f}   解析解: 0.0663")

実行結果は受容率 0.408、事後平均の推定値 0.4078（解析解 0.4074）、事後標準偏差 0.0661（解析解 0.0663）でした。小数点以下3桁まで一致しています。注目してほしいのは、このコードのどこにも $p(D)$ が現れないことです。log_unnorm_posterior は尤度と事前分布の対数を足しただけで、積分は1回も実行していません。それでも正しい事後分布が再現されました。

なお、実装で見落としがちなのが chain[t] = x の位置です。棄却したときも同じ値をもう一度記録する必要があります。棄却を「何も起きなかった」として記録を飛ばしてしまうと、滞在時間の情報が失われ、分布が歪みます。棄却は失敗ではなく、「その場に留まる」という立派な1歩なのです。

左は二峰性の混合ガウス（真の平均 0.9250、真の標準偏差 2.3701）を目標にした場合で、19万5千サンプルからの推定値は 0.9169 と 2.3739。右が上のコードのベータ事後分布で、推定 0.4078 / 0.0661 に対し解析解が 0.4074 / 0.0663 です。左のずれ 0.0081 は、後述するモンテカルロ標準誤差 0.0172 の範囲内に収まっており、統計的な揺らぎとして説明がつきます。目標分布が二峰でも、正規化定数が未知でも、ヒストグラムは正しい形に収束することが確認できました。

左の二峰性の例も、先ほど定義した metropolis 関数をそのまま流用できます。書き換えるのは目標分布だけです。

# 二峰性の目標分布: 0.35*N(-2, 0.8^2) + 0.65*N(2.5, 1.1^2)
W = np.array([0.35, 0.65]); MU = np.array([-2.0, 2.5]); SD = np.array([0.8, 1.1])

def log_mixture(x):
    """混合ガウスの対数密度。裾で桁落ちしないよう logsumexp で計算する"""
    a = np.log(W) - 0.5 * ((float(x[0]) - MU) / SD) ** 2 - np.log(SD * np.sqrt(2 * np.pi))
    m = a.max()
    return float(m + np.log(np.exp(a - m).sum()))

rng = np.random.default_rng(21)
chain, acc = metropolis(log_mixture, 0.0, 200_000, 2.5, rng)
s = chain[5000:]
print(f"受容率 {acc:.3f}  平均 {s.mean():.4f} (真値 0.9250)  標準偏差 {s.std(ddof=1):.4f} (真値 2.3701)")

出力は「受容率 0.537 平均 0.9169 標準偏差 2.3739」で、上の図の左パネルと同じ数字です。目標分布のコードを1つ差し替えるだけで、まったく同じサンプラーが使い回せるのがMH法の汎用性で、山が2つあっても3つあっても、密度を評価する関数さえ書ければ手続きは変わりません。ただし多峰性には落とし穴もあります。ここでは提案幅 2.5 が山と山の距離（4.5）と同じオーダーなので行き来できていますが、山が遠く離れていると連鎖が片方の山に閉じ込められます。この失敗の実例は次の節で見ます。

ランダムウォークMH と 独立MH

提案分布 $q$ の選び方には、大きく2つの流儀があります。

ランダムウォークMH は $q(x’|x) = \mathcal{N}(x’ \mid x, s^2 I)$ のように、現在地を中心に置く方式です。局所探索なので提案が「そこそこ良い」場所に落ちやすく、目標分布の形をまったく知らなくても動きます。上のコードもこれです。最も広く使われる標準形ですが、局所的にしか動けないので、遠く離れたモードへ移るのは苦手です。

独立MH は $q(x’|x) = q(x’)$ と、現在地に依存しない固定の提案を使う方式です。この場合MH比は $r = \big(\tilde{\pi}(x’)/q(x’)\big) \big/ \big(\tilde{\pi}(x)/q(x)\big)$ となり、重点サンプリングの重み $w = \tilde{\pi}/q$ の比という形になります。重点サンプリングとMCMCが地続きであることが、ここに透けて見えます。$q$ が目標分布をうまく覆っていれば非常に高速で、実際に先ほどの混合ガウスに $q = \mathcal{N}(0.9, 3^2)$ を使うと、受容率 58.3%、3万8千反復あたりのESSは 13844 という良い成績が出ました。ただしこれは提案が当たっている場合の話で、$q$ の裾が目標より薄いと連鎖が一点に貼りついて動かなくなります。高次元では当たりを引くのがほぼ不可能なので、実務ではランダムウォーク型が基本になります。

提案分布の「置き方」は決まりました。しかし、まだ肝心のパラメータが残っています — 提案の幅です。

提案分布のスケール問題

小さすぎても、大きすぎてもいけない

ランダムウォークMHの提案幅 $s$ は、性能を決定づける唯一にして最重要のハイパーパラメータです。両極端を考えると、なぜトレードオフになるのかがすぐに分かります。

• $s$ が小さすぎる: 提案先の密度は現在地とほとんど変わらないので、比はほぼ1、受容率はほぼ100%。全部受け入れられますが、1歩が極小なので連鎖はその場でにじむだけです。受容率が高いこと自体はまったく良い知らせではないという点に注意してください。
• $s$ が大きすぎる: 提案先は密度がほぼゼロの荒野に落ちるので、比はほぼ0、ほとんど棄却されます。連鎖は同じ場所に何百歩も貼りついたままになります。

二峰性の混合ガウス（真の平均 0.925、真の標準偏差 2.370）を目標に、3万8千サンプルで比較した結果です。$s=0.15$（左）は受容率 95.4% と一見良さそうですが、軌跡はごく狭い範囲をにじむように動くだけで、1500反復のあいだに山を1回移るのがやっとです。3万8千反復を通しても山の行き来がほとんど起こらず、ヒストグラムは左の山を大きく取りこぼしました。推定平均も 1.974 と、真値 0.925 から大きく外れています。$s=40$（右）は受容率 5.9% で、軌跡が階段状 — 棄却が続いて水平に伸びている区間 — になっています。それでもヒストグラムの形は概ね正しく、推定平均 0.894 は悪くありません。$s=2.5$（中央）は受容率 53.9% で、軌跡が全域を細かく往復し、ヒストグラムも目標に密着しています。受容率が高い左のほうが、受容率の低い右より結果は悪いという逆転が、この図の教訓です。

最適受容率 0.234

では、受容率はどのくらいを狙えばよいのでしょうか。ここで有名な理論的指針が登場します。Roberts, Gelman & Gilks (1997) は、目標分布が独立同分布な成分の積 $\pi(x) = \prod_{i=1}^d f(x_i)$ である場合について、次元 $d \to \infty$ の極限で解析しました。結論はこうです。

• 最適な提案スケールは $s^\star = 2.38 / \sqrt{d}$（成分の標準偏差が1のとき）
• そのときの受容率は $0.234$ に収束する

$\sqrt{d}$ で割るところが重要です。次元が上がるほど提案幅は縮めなければならない。$d$ 個の成分すべてを同時に「ハズさない」場所に提案する必要があるので、各成分の余裕がどんどん小さくなるためです。

10次元標準正規分布で実際に測ってみた結果です。左は受容率がスケールとともに単調に落ちる様子、中央はESS（有効サンプルサイズ、次節で説明）が山型になる様子。ESSのピークは $s = 0.75$ に出ており、理論値 $2.38/\sqrt{10} = 0.753$ と小数点以下2桁まで一致しました。右はこの2つを合成した「受容率 vs ESS」で、ピークが 0.234 の破線のすぐ上に来ています（ピーク時の実測受容率は 0.262）。10次元という決して高くない次元でも、漸近理論の予測が実用的な精度で当たることが分かります。

実務では 0.234 という数字を厳密に守る必要はなく、受容率20〜50%くらいに収まっていれば概ね良好と考えます。ESSの山は頂上付近で平たく、$s=0.75$ でのESS 877 に対し、理論値近傍の $s=0.80$ でも 859 と5%も違いません。むしろ危険なのは、受容率が 90% を超えているとき（幅が狭すぎる）と 5% を切っているとき（幅が広すぎる）です。適応的MCMC（adaptive MCMC）は、この受容率をバーンイン中に監視して $s$ を自動調整する仕組みで、Stan や PyMC のようなライブラリには標準で入っています。

スケールだけでなく「形」も合わせる

提案幅は1つの数字ですが、多次元の提案分布には本来「形」もあります。パラメータ間に強い相関があるとき、等方的な（真円の）提案は非常に効率が悪くなります。

相関 $\rho = 0.95$ の2次元ガウス分布での実測です。目標分布は細長い谷の形をしています。等方提案では、谷の幅（短軸方向）を飛び出さないように幅を小さく抑えるしかなく、その結果、長軸方向にも小さくしか進めません。$s=0.05$ ではESSがわずか27、$s=0.40$ に広げても600止まりです。一方、提案の共分散を目標の共分散に合わせた右のケース（橙の破線が目標の等高線と同じ向きに寝ている）ではESSが 7827 に跳ね上がりました。同じ反復回数で13倍の情報量です。

このため、実用的なMCMCではバーンイン期間中に標本共分散を推定して提案共分散に流用する、という手が広く使われます。あるいは、そもそも変数を事前に無相関化（再パラメータ化）しておくのも有効です。「MCMCが遅い」と感じたとき、原因が幅ではなく形にあることは非常に多い、と覚えておいてください。

チューニングの指針は見えてきました。しかし、そもそも「良い連鎖」をどう測るのか。次の節で、標本の質を定量化する道具を揃えます。

サンプルの質を測る

バーンイン — 初期値の記憶を捨てる

エルゴード定理が保証するのは $N \to \infty$ の挙動です。走り始めの数百歩は、初期値の記憶を色濃く引きずっていて、まだ $\pi$ に従っているとは言えません。この過渡状態をバーンイン期間（burn-in）と呼び、対応する標本は捨てるのが慣例です。

左は、真の分布から遠く離れた $x_0 = 60$ から出発させたときのトレースプロットです。61反復かけて目標領域まで滑り降り、そのあとは一定の帯の中で振動しています。この「毛虫が横たわっているような」見た目が、健全に混合している連鎖のサインです。右は $x$ の累積平均で、初期値が違う3本ともやがて真の平均 0.925 に収束していきます。

バーンインを捨てないとどうなるか、数字で確認しておきましょう。$x_0 = 60$ から4000反復回した場合、全反復を使った平均は 1.403 で真値 0.925 から大きく外れます。先頭500反復を捨てると 0.835 まで改善しました。過渡状態の標本はたった61個（全体の1.5%）ですが、値が真値から遠く離れているため、平均に無視できないバイアスを持ち込むのです。

捨てる長さに厳密な決まりはありません。トレースプロットを目で見て過渡が消えたところ、というのが実務的な判断で、迷ったら全体の10〜50%を捨てておけば安全側です。

自己相関と有効サンプルサイズ（ESS）

MCMCの標本は独立ではありません。連続する標本は互いによく似ているので、「1万サンプル得た」と言っても、独立標本1万個ぶんの情報量はありません。この目減りを定量化するのが自己相関時間と有効サンプルサイズ（Effective Sample Size, ESS）です。

ラグ $k$ の自己相関を $\rho_k$ とすると、積分自己相関時間 $\tau$ とESSは次で定義されます。

$$ \tau = 1 + 2\sum_{k=1}^{\infty} \rho_k, \qquad \mathrm{ESS} = \frac{N}{\tau} $$

$\tau$ は「独立な情報1個を得るのに何反復かかるか」を表します。すべての標本が独立なら $\rho_k = 0$ で $\tau = 1$、ESS $= N$。強く相関していれば $\tau$ が大きくなり、ESSは激減します。実際の推定では無限和を打ち切る必要があり、隣り合う2ラグの和が正である間だけ足す Geyer の初期正値列法などが使われます。

先ほどの3水準について、自己相関関数を並べたものです。$s=0.15$（赤）はラグ300でもまだ相関0.6を保っており、$\tau = 1826$、3万8千反復あたりのESSはわずか21。$s=2.5$（緑）はラグ20ほどでゼロに落ち、$\tau = 9.5$、ESSは 3994。$s=40$（紫）は $\tau = 34$、ESSは 1117 でした。反復回数を3万8千と報告しても、実質的な標本数は21個から3994個まで190倍の開きがあるわけです。MCMCの結果を報告するときに反復回数だけを書くのは、ほとんど意味がありません。

ESSが分かると、推定値の誤差も評価できます。モンテカルロ標準誤差（MCSE）は

$$ \mathrm{MCSE} = \frac{\hat{\sigma}}{\sqrt{\mathrm{ESS}}} $$

で計算します。先ほどの混合ガウスの例では、19万5千反復でESS 18991、標準偏差の推定 2.374 だったので MCSE は 0.0172。推定平均 0.9169 と真値 0.9250 の差 0.0081 は、この MCSE の半分以下です。「事後平均は 0.917 ± 0.017（MC誤差）」のように、推定値と一緒にMCSEを報告するのが望ましい作法です。

間引き（thinning）はしたほうがよいのか

「相関が邪魔なら、10反復に1つだけ残して独立に近づければよいのでは」と考えたくなります。これを間引き（thinning）と呼びます。結論から言うと、メモリが足りない場合を除いて、間引きは損です。

先ほどの $s=2.5$ の連鎖（3万8千サンプル、ESS 3994）で確かめてみました。5反復ごとに間引くと標本数は7600個に減り、ESSは 3665 に下がります。20反復ごとにすると標本数1900個、ESSも 1900 まで落ちました。確かに残った標本はほぼ独立になりましたが、捨てた標本が持っていた情報は戻ってきません。相関があっても、標本はゼロよりは情報を持っているからです。推定精度を決めるのは標本数ではなくESSなので、全部残すのが最善です。

間引きが正当化されるのは、標本1個あたりのメモリや後処理コストが大きい場合（高次元パラメータを何百万個も保存できない、など）に限られます。

収束診断 — トレースプロットと Gelman-Rubin の $\hat{R}$

最後に、そもそも連鎖が定常状態に達したのかをどう判定するかです。残酷な事実として、MCMCの収束を証明する方法は存在しません。診断できるのは「収束していないこと」だけで、どのテストも「まだ収束していない証拠は見つからなかった」としか言えません。それでも、実務で使える強力な診断はあります。

第一がトレースプロットです。1本の連鎖の時系列を目で見て、太い帯の中で一様に振動していれば良好、ゆっくりしたトレンドや長い階段が見えれば異常。手軽ですが、驚くほどよく効きます。

第二がGelman-Rubin 統計量 $\hat{R}$ です。アイデアは「ばらばらの初期値から $m$ 本の連鎖を走らせ、連鎖どうしが区別できなくなったら収束とみなす」というもの。$n$ 反復ぶんについて、連鎖内分散 $W$ と連鎖間分散 $B$ を

$$ W = \frac{1}{m}\sum_{j=1}^{m} s_j^2, \qquad B = \frac{n}{m-1}\sum_{j=1}^{m}\left(\bar{\theta}_j – \bar{\theta}\right)^2 $$

と定義します（$s_j^2$ は $j$ 番目の連鎖の標本分散、$\bar\theta_j$ はその平均）。これを混ぜた分散推定量

$$ \widehat{\mathrm{var}}^{+} = \frac{n-1}{n}W + \frac{1}{n}B $$

を作り、その比の平方根を取ったものが $\hat{R}$ です。

$$ \hat{R} = \sqrt{\frac{\widehat{\mathrm{var}}^{+}}{W}} $$

収束していなければ、連鎖どうしがまだ別の場所にいるので $B$ が大きく、$\hat{R} > 1$ になります。収束すれば連鎖間の差は連鎖内の揺らぎに埋もれ、$B/n \to 0$ となって $\hat{R} \to 1$ に近づきます。

提案幅をわざと狭く（0.7）して山の行き来を遅くした、混合が悪い設定での実測です。左を見ると、4本の連鎖がそれぞれ数百反復にわたって同じ山に居座っているのが分かります。右の $\hat{R}$ は $n=100$ で 1.456、$n=1000$ でもまだ 1.283 と高く、$n=3918$ でようやく実務の目安 1.01 を下回りました。3万反復まで回すと 1.0005 です。1000反復で止めていたら、間違いなく偏った事後分布を報告していたことになります。

実務上の目安は $\hat{R} < 1.01$ です（かつては 1.1 が使われていましたが、Vehtari ら (2021) の指摘を受けて基準は厳しくなりました）。連鎖は最低4本、初期値は意図的に散らして走らせるのが定石です。1本だけ走らせて $\hat{R}$ が計算できないのは、診断を放棄しているのと同じだと考えてください。

これで品質管理の道具が揃いました。最後に、MH法と並ぶもう一つの代表的MCMCとの関係を確認しましょう。

ギブスサンプリングは MH の特殊ケース

受容確率が常に1になる提案

ギブスサンプリングは、多次元パラメータを一度に動かすのではなく、成分を1つずつ順番に「他を固定した条件付き分布」から引き直す方法です。

$$ x_i’ \sim \pi(x_i \mid x_{-i}), \qquad x_{-i} = (x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d) $$

一見すると受容判定がないので、MH法とは別のアルゴリズムに見えます。しかし実は、ギブスサンプリングは提案分布に完全条件付き分布を選んだMH法であり、そのときMH比が恒等的に1になるのです。確かめてみましょう。

提案は成分 $i$ だけを動かすので、$x’ = (x_i’, x_{-i})$ で $x’_{-i} = x_{-i}$ です。$q(x’|x) = \pi(x_i’|x_{-i})$ を式 (3) のMH比に代入すると

$$ r = \frac{\pi(x’)\, q(x \mid x’)}{\pi(x)\, q(x’ \mid x)} = \frac{\pi(x_i’, x_{-i})\, \pi(x_i \mid x’_{-i})}{\pi(x_i, x_{-i})\, \pi(x_i’ \mid x_{-i})} $$

となります。ここで $x’_{-i} = x_{-i}$ なので、$\pi(x_i|x’_{-i}) = \pi(x_i|x_{-i})$ と書き換えられます。さらに同時分布を条件付き分布と周辺分布の積 $\pi(x_i, x_{-i}) = \pi(x_i|x_{-i})\,\pi(x_{-i})$ に分解して代入すると

$$ r = \frac{\pi(x_i’ \mid x_{-i})\, \pi(x_{-i})\, \pi(x_i \mid x_{-i})}{\pi(x_i \mid x_{-i})\, \pi(x_{-i})\, \pi(x_i’ \mid x_{-i})} = 1 $$

分子と分母がまったく同じ3つの因子の積になり、$r = 1$、したがって $\alpha = \min(1,1) = 1$ です。ギブスサンプリングで棄却が起きないのは偶然ではなく、提案の選び方から必然的にそうなるわけです。

相関 $\rho=0.95$ の2次元ガウスでの比較です。ギブス（左）は一度に1成分しか更新しないので、軌跡が必ず座標軸に平行な階段状になります。MH（中央）は斜めにも動けますが、棄却されるとその場に留まります。1万9千反復あたりのESSはギブスが 1023、等方提案のMH（提案幅0.4、受容率51%）が 195 で、この目標分布ではギブスの圧勝でした。棄却がない分、そして条件付き分布という「完璧な提案」を使える分、有利に働いています。

ただしギブスにも弱点があります。第一に、完全条件付き分布が解析的に求まる必要があります。共役事前分布を使った線形モデルや混合モデルなら綺麗に出ますが、そうでなければ使えません。第二に、軸に平行にしか動けないので、上の図のような細長い谷では階段を1段ずつ登ることになり、相関がさらに強い（$\rho \to 1$）と急速に効率が落ちます。MH法は提案の設計が必要な代わりに、どんな分布にも適用できる汎用性があります。

両者の性格の違いを表にまとめておきます。

要するに、「常に受容できてチューニングも要らないが、条件付き分布が導出できなければ使えないギブス」と「どんな分布にも使えるがチューニングが要るMH」という住み分けです。実務では両者を組み合わせることも多く、条件付き分布が求まる成分はギブスで、求まらない成分だけMHで更新する Metropolis-within-Gibbs が定番の折衷案になっています。

ギブスサンプリングを混合ポアソン分布のパラメータ推定に適用した具体例は、次の記事で詳しく扱っています。

混合ポアソン分布のパラメータをギブスサンプリングで求める

潜在変数とモデルパラメータを分けてサンプリングする、混合モデルへのギブスサンプリングの具体的な適用例。完全条件付き分布の導出まで追える。

その先へ — 勾配を使うMCMC

ランダムウォークMHの弱点は、提案が完全に無情報だという点にあります。目標分布の形を一切見ずに、ただ現在地の周りに等方的に候補をばらまいているので、次元が上がるほど当たりにくくなる。$2.38/\sqrt{d}$ というスケーリング則は、まさにその代償です。

そこで、目標分布の勾配 $\nabla \log \pi(x)$ を提案に組み込むのが現代の主流です。

• MALA（Metropolis-adjusted Langevin algorithm）: 提案を $x’ \sim \mathcal{N}\big(x + \tfrac{\epsilon^2}{2}\nabla \log \pi(x),\, \epsilon^2 I\big)$ とし、密度が高い方向へ少し押し出します。提案が非対称になるので、MH比では $q$ の因子を省略できません。最適スケーリングは $\epsilon \propto d^{-1/6}$ で、$d^{-1/2}$ のランダムウォークより次元に強くなります。
• HMC（Hamiltonian Monte Carlo）: 補助的な運動量変数を導入し、$-\log\pi$ をポテンシャルとするハミルトン力学の軌道を数値積分して遠くまで一気に移動します。離散化誤差をMH判定で補正するので、これも本質的にはMH法の一種です。最適スケーリングは $d^{-1/4}$。
• NUTS（No-U-Turn Sampler）: HMCの軌道長を自動で決める拡張。Stan や PyMC のデフォルトサンプラーで、実務でMCMCを使うなら実際に動かすのはこれになることが多いはずです。

面白いのは、これらがすべて「提案分布 $q$ を賢くしたMH法」として同じ枠組みに収まることです。詳細釣り合いと $\alpha = \min(1, \text{比})$ という骨格は変わりません。本記事で導いた式が、現代のベイズ計算ライブラリの心臓部でもそのまま動いているわけです。実際にPyMCでこれらを使う方法はPyMCによる確率的プログラミングで扱っています。

まとめ

本記事では、MCMCの一般論とメトロポリス・ヘイスティングス法について解説しました。

• 動機: ベイズ事後分布の正規化定数 $p(D)$ は高次元積分で計算できない。しかし2点での密度の比なら計算できる。棄却サンプリングの受容率は次元とともに指数的に落ちる（10次元で 0.25%、20次元で $2.5\times10^{-8}$）ため、独立サンプリングは高次元で使えない。
• 発想: 独立性をあきらめ、目標分布 $\pi$ を不変分布に持つマルコフ連鎖を作って走らせる。既約性・非周期性が揃えば、エルゴード定理により時間平均が空間平均に収束する。
• 詳細釣り合い: $\pi(x)P(x\to x’) = \pi(x’)P(x’\to x)$。両辺を $x$ で積分すれば $\pi$ が不変分布であることが数行で示せる。大域的な積分方程式を、2点ずつの局所条件に置き換えるのが設計上の急所。
• 受容確率: 遷移を「提案 $q$ ＋受容判定 $\alpha$」に分解し、詳細釣り合いから $\alpha = \min\big(1, \pi(x’)q(x|x’) / \pi(x)q(x’|x)\big)$ を導いた。$\pi$ が比の形でしか現れないので、正規化定数が約分で消える。提案が対称なら $\alpha = \min(1, \pi(x’)/\pi(x))$ で、「上り坂は必ず受容、下り坂は高さの比で受容」。
• チューニング: 提案幅は小さすぎても大きすぎても情報を稼げない。高次元の漸近理論では最適スケール $2.38/\sqrt{d}$、最適受容率 0.234。10次元での実測でもピークは $s=0.75$（理論値 0.753）、受容率 0.262 だった。スケールだけでなく提案の形（共分散）を目標に合わせると、相関 0.95 の例でESSが 600 から 7827 へ13倍に改善する。
• 品質管理: バーンインで初期値の記憶を捨てる。反復回数ではなくESS $= N/\tau$ で実質的な標本数を測り、MCSE $= \hat\sigma/\sqrt{\mathrm{ESS}}$ で誤差を報告する。間引きはESSを減らすだけなので基本的に不要。収束は複数チェーンの $\hat{R} < 1.01$ で診断する（収束の証明はできず、非収束の検出しかできない点に注意）。
• ギブスとの関係: 提案に完全条件付き分布を使うとMH比が恒等的に1になり、常に受容される。ギブスはMH法の特殊ケース。

MCMCの美しさは、「正規化定数が計算できない」という致命的に見えた障害を、「比を取れば消える」という一手でまるごと迂回してしまうところにあります。そしてその一手を正当化しているのが、たった数行で証明できる詳細釣り合い条件です。ベイズ統計が理論の枠を超えて実務の道具になったのは、この単純な仕掛けのおかげだと言っても大げさではありません。

次のステップとして、以下の記事も参考にしてください。

• ギブスサンプリングの仕組みをわかりやすく解説して実装する — MH法の受容確率が常に1になる特殊ケース
• PyMCによる確率的プログラミング — NUTSを実際に動かしてベイズモデルを推定する
• ベイズ混合ガウスモデルとサンプリング — MCMCを適用する典型的なモデル
• 粒子フィルタの仕組みを理解してPythonで実装する — 逐次データに対するモンテカルロ法
• 予測分布をわかりやすく解説 — MCMC標本から予測分布を作る