はじめに
ベイズ推定では、事後分布 \(p(\theta | D)\) からのサンプリングが中心的な課題です。しかし、多くの場合この分布は解析的に扱えず、直接サンプリングもできません。
**マルコフ連鎖モンテカルロ（MCMC）**法は、目標分布を定常分布とするマルコフ連鎖を構築し、そこからサンプルを生成する手法です。十分長い連鎖を走らせることで、事後分布からの近似サンプルが得られます。
マルコフ連鎖の基礎
マルコフ連鎖は、現在の状態のみに依存して次の状態が決まる確率過程です。
\[P(X_{t+1} | X_1, \ldots, X_t) = P(X_{t+1} | X_t) \tag{1}\]MCMCの目的は、目標分布 \(\pi(\theta)\) を定常分布とするマルコフ連鎖を設計することです。定常分布とは、連鎖が十分長く走った後に収束する分布であり、**詳細釣り合い条件（detailed balance）**を満たす遷移確率 \(T(\theta' | \theta)\) を設計すれば保証されます。
\[\pi(\theta) T(\theta' | \theta) = \pi(\theta') T(\theta | \theta') \tag{2}\]なぜ詳細釣り合いが定常性を保証するのか
マルコフ連鎖が分布 \(\pi\) を定常分布として持つとは、次の関係が成り立つことを指します。
\[\int \pi(\theta) T(\theta' | \theta) \, d\theta = \pi(\theta') \tag{2'}\]左辺は「現在の状態が \(\pi\) に従う」と仮定したとき、1ステップ遷移後の状態 \(\theta'\) の周辺分布を表します。これが元の \(\pi(\theta')\) と一致するなら、遷移を繰り返しても分布の形は変わりません。
詳細釣り合い条件 (2) の両辺を \(\theta\) について積分すると、
\[\int \pi(\theta) T(\theta' | \theta) \, d\theta = \int \pi(\theta') T(\theta | \theta') \, d\theta = \pi(\theta') \int T(\theta | \theta') \, d\theta = \pi(\theta')\]が得られ、式(2’)がそのまま成立することが分かります。最後の等号は、\(T(\theta | \theta')\) が \(\theta'\) を出発点とする遷移核であり、到達先 \(\theta\) に関して積分すれば1になる（確率の正規化条件）ことを使っています。したがって、詳細釣り合いを満たす遷移核は自動的に定常性の式を満たします。
逆は成り立たない点に注意してください。定常分布であっても詳細釣り合い（可逆性）を満たすとは限りません。しかし可逆性は定常性を保証する十分条件であり、かつ遷移核を具体的に設計する上で扱いやすいため、MCMCの主要な手法（メトロポリス・ヘイスティングス法、ギブスサンプリング、後述のHMCなど）はいずれもこの可逆性を経由して定常分布を保証する設計になっています。
メトロポリス・ヘイスティングス法
アルゴリズム
- 初期値 \(\theta_0\) を設定
- 提案分布 \(q(\theta' | \theta_t)\) から候補 \(\theta'\) を生成
- 受容確率を計算:
- 確率 \(\alpha\) で \(\theta_{t+1} = \theta'\) 、それ以外は \(\theta_{t+1} = \theta_t\)
- ステップ2-4を繰り返す
提案分布が対称（\(q(\theta' | \theta) = q(\theta | \theta')\) ）の場合、受容確率は次のように簡略化されます（メトロポリス法）。
\[\alpha = \min\left(1, \frac{\pi(\theta')}{\pi(\theta_t)}\right) \tag{4}\]重要な性質として、\(\pi(\theta)\) の正規化定数を知る必要がありません。ベイズ推定では \(p(\theta | D) \propto p(D | \theta) p(\theta)\) の非正規化密度のみで計算可能です。
受容確率が詳細釣り合いを満たすことの証明
式(3)の受容確率がなぜこの特定の形をしているのか、天下り的に受け入れるのではなく証明します。MH連鎖の遷移核は、提案 \(\theta' \neq \theta\) が採用される成分と、棄却されて留まる成分からなり、\(\theta' \neq \theta\) については
\[T(\theta' | \theta) = q(\theta' | \theta) \, \alpha(\theta' | \theta)\]と書けます。したがって、\(\theta' \neq \theta\) の場合に詳細釣り合い \(\pi(\theta) T(\theta' | \theta) = \pi(\theta') T(\theta | \theta')\) が成り立つことを示せば十分で、これは
\[\pi(\theta) \, q(\theta' | \theta) \, \alpha(\theta' | \theta) = \pi(\theta') \, q(\theta | \theta') \, \alpha(\theta | \theta')\]という等式に帰着します。ここで比 \(r = \dfrac{\pi(\theta') q(\theta | \theta')}{\pi(\theta) q(\theta' | \theta)}\) を定義し、式(3)に従って \(\alpha(\theta' | \theta) = \min(1, r)\) 、\(\alpha(\theta | \theta') = \min(1, 1/r)\) とおきます。
\(r \le 1\) の場合：\(\alpha(\theta' | \theta) = r\) かつ \(\alpha(\theta | \theta') = 1\) なので、
\[\text{左辺} = \pi(\theta) q(\theta' | \theta) \cdot r = \pi(\theta) q(\theta' | \theta) \cdot \frac{\pi(\theta') q(\theta | \theta')}{\pi(\theta) q(\theta' | \theta)} = \pi(\theta') q(\theta | \theta') = \text{右辺}\]\(r > 1\) の場合は、\(\alpha(\theta' | \theta) = 1\) 、\(\alpha(\theta | \theta') = 1/r\) となり、対称な計算により両辺とも \(\pi(\theta) q(\theta' | \theta)\) に一致することが確認できます。\(\theta' = \theta\) の場合は両辺とも「留まる確率」の同一の項が現れるため自明に等しくなります。
以上より、式(3)の受容確率はあらゆるケースで詳細釣り合いを満たします。つまりこの式は任意に選ばれた公式ではなく、詳細釣り合いという制約から逆算的に導かれる、可能な受容確率の中で棄却をもっとも減らす（\(\alpha\) を可能な限り1に近づける）設計になっています。
Python実装
正規分布の事後分布からサンプリングする例です。
import numpy as np
import matplotlib.pyplot as plt
def metropolis_hastings(log_target, initial, n_samples, proposal_std=1.0):
"""メトロポリス・ヘイスティングス法"""
samples = [initial]
current = initial
accepted = 0
for _ in range(n_samples):
# 提案分布（ガウスランダムウォーク）
proposal = current + np.random.normal(0, proposal_std)
# 対数受容確率
log_alpha = log_target(proposal) - log_target(current)
# 受容/棄却
if np.log(np.random.random()) < log_alpha:
current = proposal
accepted += 1
samples.append(current)
acceptance_rate = accepted / n_samples
return np.array(samples), acceptance_rate
# --- 目標分布: 混合ガウス分布 ---
def log_target(x):
"""対数目標密度: 0.3*N(-2,1) + 0.7*N(3,0.5)"""
from scipy.special import logsumexp
log_p1 = np.log(0.3) - 0.5 * (x + 2)**2
log_p2 = np.log(0.7) - (x - 3)**2 # sigma=0.5 → 1/(2*0.25)=2
return logsumexp([log_p1, log_p2])
# --- 実行 ---
np.random.seed(42)
samples, acc_rate = metropolis_hastings(log_target, initial=0.0,
n_samples=50000, proposal_std=1.5)
print(f"受容率: {acc_rate:.3f}")
# バーンイン期間を除去
burn_in = 5000
samples = samples[burn_in:]
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# トレースプロット
axes[0].plot(samples[:2000], alpha=0.5, linewidth=0.5)
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Sample value')
axes[0].set_title('Trace Plot')
axes[0].grid(True, alpha=0.3)
# ヒストグラム vs 真の分布
x = np.linspace(-6, 6, 200)
true_pdf = 0.3 * np.exp(-0.5 * (x + 2)**2) / np.sqrt(2*np.pi) + \
0.7 * np.exp(-0.5 * ((x - 3)/0.5)**2) / (0.5 * np.sqrt(2*np.pi))
axes[1].hist(samples, bins=100, density=True, alpha=0.5, label='MCMC samples')
axes[1].plot(x, true_pdf, 'r-', linewidth=2, label='True density')
axes[1].set_xlabel('x')
axes[1].set_ylabel('Density')
axes[1].set_title('Histogram vs True Distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
ギブスサンプリング
多次元の場合、各変数を条件付き分布から順番にサンプリングする手法です。2変数 \((\theta_1, \theta_2)\) の場合:
- \(\theta_1^{(t+1)} \sim p(\theta_1 | \theta_2^{(t)}, D)\)
- \(\theta_2^{(t+1)} \sim p(\theta_2 | \theta_1^{(t+1)}, D)\)
条件付き分布が既知の場合（共役事前分布など）に特に有効です。メトロポリス・ヘイスティングス法の特殊ケースであり、受容率は常に1です。
メトロポリス・ヘイスティングス法の特殊ケースとしての証明
「受容率が常に1」という性質は天下り的な事実ではなく、MHの枠組みから導けます。同時分布 \(\pi(\theta_1, \theta_2)\) に対して、\(\theta_1\) の更新における提案分布を完全条件付き分布そのものに取ります。
\[ q(\theta_1' | \theta_1, \theta_2) = \pi(\theta_1' | \theta_2) \]（\(\theta_2\) を固定したまま、\(\theta_1\) だけを条件付き分布から直接サンプルする提案です。逆向きの提案 \(q(\theta_1 | \theta_1', \theta_2)\) も同じ条件付き分布 \(\pi(\theta_1 | \theta_2)\) に従います。）
同時分布は条件付き分布と周辺分布の積 \(\pi(\theta_1, \theta_2) = \pi(\theta_1 | \theta_2) \pi(\theta_2)\) に分解できます。これをMHの受容確率(3)に代入すると、
\[ \alpha = \min\left(1, \frac{\pi(\theta_1', \theta_2) \, q(\theta_1 | \theta_1', \theta_2)}{\pi(\theta_1, \theta_2) \, q(\theta_1' | \theta_1, \theta_2)}\right) = \min\left(1, \frac{\pi(\theta_1' | \theta_2) \pi(\theta_2) \cdot \pi(\theta_1 | \theta_2)}{\pi(\theta_1 | \theta_2) \pi(\theta_2) \cdot \pi(\theta_1' | \theta_2)}\right) = \min(1, 1) = 1 \]分子・分母に現れる \(\pi(\theta_1' | \theta_2)\) 、\(\pi(\theta_1 | \theta_2)\) 、\(\pi(\theta_2)\) がすべて打ち消し合うため、比は恒等的に1になります。これが「ギブスサンプリングの提案は必ず受容される」ことの証明です。条件付き分布から直接サンプルするというアイデア自体が、すでにMHの受容確率を最大化するように「最適な提案分布」を選んでいることに相当します。
エルゴード性と収束の理論
詳細釣り合いは定常分布の存在を保証しますが、それだけでは「連鎖が実際にその定常分布へ収束する」ことは保証されません。収束を保証するには、連鎖が次の2条件を満たす必要があります。
既約性（irreducibility）
任意の状態 \(\theta\) から、正の確率を持つ任意の状態 \(\theta'\) （\(\pi(\theta') > 0\) ）へ、有限回の遷移で到達できる確率が正であることを言います。式で書けば、ある回数 \(n\) が存在して \(n\) ステップ遷移核 \(T^{(n)}\) について
\[ T^{(n)}(\theta' | \theta) > 0 \]が成り立つことです。既約性が破れると、連鎖は状態空間の一部に永久に閉じ込められ、目標分布の一部の質量に到達できません。
非周期性（aperiodicity）
連鎖がある状態に戻ってくるまでのステップ数の集合の最大公約数が1であることを言います。周期 \(d > 1\) を持つ連鎖は、時刻 \(t\) における分布が \(d\) 種類の分布の間を循環し、単一の極限分布へ収束しません（例えば偶数ステップと奇数ステップで訪れる状態集合が排他的に切り替わるような設計をすると発生します）。ランダムウォーク型の提案分布（連続分布からの加法的ノイズ）を使う限り、通常は非周期性は自動的に満たされます。
エルゴード定理
連鎖が既約かつ非周期的で、\(\pi\) を定常分布として持つなら、任意の初期分布から出発しても、時刻 \(n\) における分布は全変動距離で \(\pi\) へ収束します。
\[ \lim_{n \to \infty} \left\| P(\theta_n \in \cdot) - \pi(\cdot) \right\|_{TV} = 0 \]さらに大数の法則の拡張として、任意の可積分関数 \(f\) について、時間平均が真の期待値へほぼ確実に収束します（エルゴード平均）。
\[ \frac{1}{N} \sum_{t=1}^{N} f(\theta_t) \ \xrightarrow{N \to \infty} \ \mathbb{E}_{\pi}[f(\theta)] \]これがMCMCの理論的な正当化の核心です。ただしこの定理が保証するのは「いつか収束すること」だけであり、「どれだけ速く収束するか」（混合時間）は一切保証しません。以下の2つの具体例で、この違いが実務上いかに重要かを確認します。
具体例1：既約性が厳密に破れるケース
目標分布として、互いに離れた2つの区間上の一様分布 \(\pi(\theta) = 0.5 \cdot \mathrm{Unif}[0, 1] + 0.5 \cdot \mathrm{Unif}[3, 4]\) を考えます。提案分布として、有界な一様ランダムウォーク \(q(\theta' | \theta) = \mathrm{Unif}(\theta - 0.4, \theta + 0.4)\) （1ステップの移動幅が高々0.4）を使うと、区間 \([0, 1]\) 内のどの点からも、提案の到達範囲は高々 \((-0.4, 1.4)\) に収まり、区間 \([3, 4]\) には一度も到達できません。これは近似の問題ではなく、任意の \(n\) について \(T^{(n)}(\theta' | \theta) = 0\) （\(\theta \in [0,1]\) 、\(\theta' \in [3,4]\) ）が厳密に成り立つ、既約性の完全な破れです。連鎖は初期値が属する区間に永久に閉じ込められ、目標分布の残り50%の質量を完全に見逃します。シミュレーションするまでもなく、提案の到達範囲の幾何学から自明に導ける例です。
具体例2：実務で起きる「実質的な」閉じ込め
上記は理論上厳密に既約性が破れる例ですが、実務でより頻繁に問題になるのは、理論上は既約（ガウス提案は \(\mathbb{R}\) 全域に正の密度を持つため、どの状態にも正の確率で到達できる）でありながら、ステップサイズが目標分布の谷の幅に対して小さすぎるために、現実的な計算時間内では実質的に閉じ込められるケースです。
これを実際に確認するため、本記事のPython実装で使った混合ガウス分布 \(0.3 \cdot \mathcal{N}(-2, 1) + 0.7 \cdot \mathcal{N}(3, 0.5)\)
を目標分布とし、初期値 \(\theta_0 = -2\)
（劣モード側）から、ステップサイズ（提案分布の標準偏差）0.1と1.5でそれぞれ2000イテレーション実行しました（乱数シード固定、proposal_std 以外は同一条件）。
| ステップサイズ | 受容率 | モード間の遷移回数 | サンプルの最大値 | \(x > 0.5\) の割合 |
|---|---|---|---|---|
| 0.1 | 0.970 | 0 | 0.127 | 0.0% |
| 1.5 | 0.567 | 48 | 5.442 | 51.5% |
ステップサイズ0.1では、受容率97.0%と一見「効率よく」動いているように見えますが、2000イテレーションの間一度も優勢モード（\(x=3\) 付近）に到達せず、劣モード周辺（サンプル平均 \(-1.97\) ）に完全に閉じ込められました。一方ステップサイズ1.5では48回モード間を行き来し、両モードを探索できています。
さらに、ステップサイズ1.5・初期値 \(-2\) のまま反復回数を伸ばして \(x > 0.5\) の割合の収束を追跡すると（この密度の真の値は後述の通り約0.625）:
| イテレーション数 | \(x > 0.5\) の割合 |
|---|---|
| 2,000 | 0.515 |
| 100,000 | 0.587 |
| 2,000,000 | 0.624 |
「良い」ステップサイズを選んでいても、多峰分布では真の混合比率に収束するまでに数十万〜数百万イテレーション規模の時間を要することが実測できます。これはエルゴード定理が保証する「いつか収束する」ことと、「実務で許容できる時間内に収束する」ことの間にある大きなギャップを示す実例です。
実装上の注意点（正規化定数の落とし穴）：本記事のPythonコードの log_target 関数は、\(\log p_1 = \log(0.3) - 0.5(x+2)^2\)
、\(\log p_2 = \log(0.7) - (x-3)^2\)
という形で各成分の指数部のみを実装しており、各正規分布本来の正規化定数（\(\mathcal{N}(-2,1)\)
の \(1/\sqrt{2\pi}\)
、\(\mathcal{N}(3, 0.5)\)
の \(1/\sqrt{\pi}\)
）を含んでいません。MH法は目標密度全体に共通する正規化定数を知る必要はありませんが、混合分布の各成分ごとの相対的な正規化は結果の混合比率に直接影響します。実際にこの密度を数値積分すると、コード上の「0.3」「0.7」というラベルにもかかわらず、真の混合比率は約 \(0.377 : 0.623\)
になります（\(x > 0.5\)
となる確率は約0.625）。上表の実験値がこの真値へ収束していく様子と整合しているのはこのためです。混合密度を実装する際は、各成分の正規化定数を必ず含めるべき典型的な落とし穴と言えます。
収束診断
エルゴード定理は連鎖がいつか収束することを保証しますが、実際に手元の有限長の連鎖が「十分収束したか」を判定する数学的な手段が別途必要です。ここでは実務で標準的に使われる2つの診断法を、実際にPythonで計算した数値とともに紹介します。
Gelman-Rubin統計量 \(\hat{R}\)
異なる初期値から複数の独立な連鎖（\(m\) 本、各長さ\(n\) ）を走らせ、チェーン間分散とチェーン内分散を比較する診断法です。チェーン \(j\) の平均を \(\bar{\theta}_j\) 、全体平均を \(\bar{\theta}\) として、チェーン内分散の平均
\[ W = \frac{1}{m} \sum_{j=1}^{m} s_j^2, \qquad s_j^2 = \frac{1}{n-1} \sum_{t=1}^{n} (\theta_{j,t} - \bar{\theta}_j)^2 \]と、チェーン間分散（\(n\) 倍してチェーン内分散とスケールを揃える）
\[ B = \frac{n}{m-1} \sum_{j=1}^{m} (\bar{\theta}_j - \bar{\theta})^2 \]を計算し、真の分散のプールされた推定量 \(\widehat{\mathrm{Var}}(\theta) = \frac{n-1}{n} W + \frac{1}{n} B\) を使って、
\[ \hat{R} = \sqrt{\frac{\widehat{\mathrm{Var}}(\theta)}{W}} \]と定義します。すべての連鎖が定常分布に収束していれば、各連鎖の平均 \(\bar{\theta}_j\) はどれも同じ \(\pi\) の平均に近づくため \(B\) は小さくなり、\(\hat{R} \to 1\) となります。逆に連鎖がそれぞれ異なる初期値の影響を引きずっている場合、\(\bar{\theta}_j\) 同士がばらつくため \(B\) が大きくなり \(\hat{R} \gg 1\) になります。経験則として \(\hat{R} < 1.1\) （厳しい基準では\(1.01\) ）を収束の目安とします。
実行検証：混合ガウス分布に対し、極端に離れた初期値 \(\theta_0 \in \{-10, -5, 5, 10\}\) から4本の連鎖（ステップサイズ1.5）を走らせ、連鎖の長さ \(n\) を変えながら（バーンインは除去せず）\(\hat{R}\) を計算しました。
| 連鎖の長さ \(n\) | \(\hat{R}\) |
|---|---|
| 50 | 1.891 |
| 100 | 1.851 |
| 300 | 1.102 |
| 1,000 | 1.016 |
| 5,000 | 1.001 |
| 20,000 | 1.001 |
\(n=50\) 〜\(100\) では初期値の違いがまだ強く残っており \(\hat{R}\) は1.85〜1.89という非常に大きい値を示しますが、\(n=300\) で1.10まで急減し、\(n=1000\) 以降は1.02未満に収まります。極端な初期値からのバーンインが数十〜数百イテレーション程度で解消されるという、後述のバーンイン検証の結果とも整合する挙動です。
さらに、連鎖長を \(n=2000\) に固定し、バーンイン量を変えて \(\hat{R}\) を比較すると:
| バーンイン | \(\hat{R}\) |
|---|---|
| 0 | 1.006 |
| 100 | 1.002 |
| 1,000 | 1.002 |
\(n=2000\) まで伸ばすと、バーンインをまったく除去しなくても \(\hat{R}\) はすでに1.01を下回っており、この設定では十分な連鎖長そのものが初期値依存性を薄める役割を果たしていることが分かります（ただしこれは各連鎖が数十イテレーションで良い領域に到達する今回の例に特有の結果であり、一般には短い連鎖でのバーンイン除去は依然として重要です）。
有効サンプルサイズ（ESS）
MCMCサンプルは自己相関を持つため、\(n\) 個のサンプルは\(n\) 個の独立サンプルと同じ情報量を持ちません。有効サンプルサイズ（ESS）は、実際に得られたサンプル列と同じ分散の推定精度を与える「仮想的な独立サンプル数」です。ラグ\(k\) の自己相関を \(\rho_k\) として、積分自己相関時間
\[ \tau_{\mathrm{int}} = 1 + 2 \sum_{k=1}^{\infty} \rho_k \]を定義すると、\(\mathrm{ESS} = n / \tau_{\mathrm{int}}\) となります。実務では無限和を打ち切る必要があり、本記事ではGeyerの初期正数列規則（\(\rho_{2k} + \rho_{2k+1}\) の対を、正である限り足し合わせて打ち切る方法）を用いました。次節でこのESSをステップサイズ選択の指標として使います。
実用上の注意点
バーンイン
マルコフ連鎖の初期値の影響を除去するため、最初の数千サンプル（バーンイン期間）を捨てます。この必要性を実際に確認するため、混合ガウス分布に対し、極端に離れた5つの初期値 \(\theta_0 \in \{-20, -10, 0, 10, 20\}\) （目標分布の主要な質量はおよそ \([-4, 5]\) の範囲にあります）から連鎖を走らせ、「妥当な領域 \([-4, 5]\) に初めて到達するまでのイテレーション数」を計測しました（ステップサイズ1.5）。
| 初期値 \(\theta_0\) | 妥当な領域に到達するまでの反復数 |
|---|---|
| -20 | 39 |
| -10 | 6 |
| 0 | 0 |
| 10 | 13 |
| 20 | 24 |
いずれの初期値でも40イテレーション以内に妥当な領域へ到達しており、本記事のPython実装で採用したバーンイン期間5000サンプルは、この種の低次元・単純な目標分布に対しては十分すぎるほど保守的な設定であることが確認できます。ただし、高次元・多峰性が強い、あるいは変数間の相関が強い事後分布では、妥当な領域への到達自体には数十イテレーションで済んでも、真の定常分布の統計量（分散や高次モーメントなど）に収束するにはさらに長い時間を要する点に注意してください。前節の \(\hat{R}\) やESSは「妥当な領域に到達したか」ではなく「統計量として収束した度合い」を定量化する診断法であり、両者は区別して考える必要があります。
受容率の調整
提案分布の標準偏差が大きすぎると受容率が低下し、小さすぎると探索が遅くなります。この2つの失敗モードのトレードオフを定量的に確認するため、標準正規分布 \(\mathcal{N}(0,1)\) を目標分布として、ステップサイズ（提案標準偏差）を0.1から10まで変化させ、50000サンプル生成後にバーンイン5000を除いた45000サンプルについて受容率とESSを計算しました。
| ステップサイズ | 受容率 | ESS | ESS比率 | 積分自己相関時間 \(\tau_{\mathrm{int}}\) |
|---|---|---|---|---|
| 0.1 | 96.8% | 95.9 | 0.21% | 469.4 |
| 0.5 | 84.6% | 1,874.5 | 4.17% | 24.0 |
| 1.0 | 70.8% | 5,189.1 | 11.53% | 8.67 |
| 2.38 | 44.8% | 10,508.6 | 23.35% | 4.28 |
| 5.0 | 24.4% | 6,918.7 | 15.37% | 6.50 |
| 10.0 | 12.4% | 3,509.6 | 7.80% | 12.82 |
ESSはステップサイズに対して山型（逆U字型）の関係を示し、ステップサイズ2.38（受容率44.8%）で最大値10,508.6（サンプル数の23.35%）を記録しました。ステップサイズが小さすぎると（0.1）ほぼすべての提案が受容されますが、1歩の移動量が小さいため強い自己相関が残り（\(\tau_{\mathrm{int}} \approx 469\) ）、ESSはサンプル数の0.2%程度にまで低下します。逆に大きすぎると（10.0）提案のほとんどが低確率領域に落ちて棄却され、受容率は12.4%まで落ち込み、ESSも再び低下します。
この最適点における受容率44.8%は、Roberts, Gelman & Gilks (1997) が導出した、1次元の独立同分布ターゲットに対する理論的最適受容率**約44%と非常によく一致しています（多次元・高相関の場合はこの最適値が次元数の増加とともに約23.4%**へ漸近することが知られています）。「受容率を大きくしすぎず小さくしすぎず」という実務上の経験則が、実際に数値実験で再現される結果と言えます。
自己相関と間引き
連続するサンプルは相関を持ちます。独立なサンプルが必要な場合は間引き（thinning、\(k\) 個ごとにサンプルを保持する）を行います。ただし、前節で定義したESSは間引きをしなくても分散などの統計量を正しく推定する上で必要な情報量をすでに定量化しており、実務では間引きよりもESSに基づいて必要なサンプル数を見積もる方が効率的です。間引きはストレージ容量の削減には有効ですが、捨てたサンプルの情報は完全に失われるため、分散推定の効率という観点では間引かずにすべてのサンプルを使う方が（同じ計算コストの下では）優れています。
関連記事
- Hamiltonian Monte Carlo（HMC）の理論とPython実装 - 本記事で触れた自己相関の問題を、勾配情報とハミルトン力学で解決する発展編です。相関ガウス分布で約45倍の効率向上を実測しています。
- 焼きなまし法（Simulated Annealing）の理論とPython実装 - MCMCと同じメトロポリス基準を最適化に応用した手法を解説しています。
- ベイズ最適化の基礎とPython実装 - MCMCで推定した事後分布をベイズ最適化に活用できます。
- 粒子フィルタのPython実装 - 逐次モンテカルロ法（SMC）であり、MCMCと密接に関連しています。
- ベイズ線形回帰の基礎 - MCMCを使わない解析的なベイズ推定を解説しています。
- クロスエントロピー法：モンテカルロ最適化の実践的手法 - モンテカルロ法を最適化に応用したCEMとの関連を理解できます。
- ガウス過程回帰の理論とPython実装 - ベイズ推論の代表例。MCMCの一般的なサンプリング戦略と、GPの解析的事後分布のアプローチを対比できます。
参考文献
- Metropolis, N., et al. (1953). “Equation of State Calculations by Fast Computing Machines”. The Journal of Chemical Physics, 21(6), 1087-1092.
- Hastings, W. K. (1970). “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”. Biometrika, 57(1), 97-109.
- Gelman, A., et al. (2013). Bayesian Data Analysis (3rd ed.). Chapman and Hall/CRC. Chapters 11-12.
- Gelman, A., & Rubin, D. B. (1992). “Inference from Iterative Simulation Using Multiple Sequences”. Statistical Science, 7(4), 457-472.
- Roberts, G. O., Gelman, A., & Gilks, W. R. (1997). “Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms”. The Annals of Applied Probability, 7(1), 110-120.
- Geyer, C. J. (1992). “Practical Markov Chain Monte Carlo”. Statistical Science, 7(4), 473-483.