マルコフ連鎖をわかりやすく解説【MCMC法への応用】
1. マルコフ連鎖の基本原理
1.1 マルコフ連鎖
マルコフ連鎖 (Markov Chain) は、確率過程の一種です。離散時間マルコフ連鎖を例にとると、時刻 ${t = 0, 1, 2, \dots}$ に観測される状態 ${X_t}$がとりうる値を状態空間${S}$と呼び、状態遷移が以下の「マルコフ性（メモリレス性）」を満たします。
$${P(X_{t+1} = j \mid X_t = i, X_{t-1} = i_{t-1}, \dots, X_0 = i_0) = P(X_{t+1} = j \mid X_t = i) = p_{ij}}$$
${p_{ij}}$ は「状態 ${i}$ から状態 ${j}$ に遷移する確率」です。
つまり、一度現在の状態がわかれば、将来の状態の確率分布は過去の履歴には依存しない、というのがマルコフ連鎖の大きな特徴です。
1.2 定常分布とエルゴード性
定常分布
マルコフ連鎖の長期的なふるまいを理解する上でだいじな概念が「定常分布」です。
定常分布 ${\pi}$ は
$${ \pi = \pi P}$$
を満たす確率分布です。
${P は推移確率行列 (${p_{ij}}$) を並べた行列です。
もしマルコフ連鎖が何らかの条件でこの分布に収束すれば、長期的に見ると状態は ${\pi}$ に従って分布するということになります。
この辺は、統計検定準一級で出ていますね。
エルゴード性
エルゴード性は「あるマルコフ連鎖が、初期状態に依存しない一意の定常分布に収束する性質」を指します。
エルゴード的なマルコフ連鎖であれば、十分に時間を進めたあとに得られるサンプルは、初期値の影響をほぼ受けなくなるため、定常分布を“代表”する標本として利用できます。
2. MCMCとマルコフ連鎖
2.1 なぜマルコフ連鎖がサンプリングに使われるのか
ベイズ統計や各種推論タスクでは「複雑な確率分布（事後分布など）からランダムサンプリングしたい」という場面が頻繁に登場します。
しかし、高次元であったり、正規化定数が不明であったりする分布から直接サンプルを得るのは難しい場合が多いです。
そこで着目されるのが、エルゴード的なマルコフ連鎖の性質です。
適切に設計したマルコフ連鎖が、ある確率分布 ${\pi}$（これを目標分布と呼びます）を定常分布として持つならば、
- 十分にたくさん連鎖を進めると、理論的には初期状態の影響が無視できるレベルになる。
- その後に得られるサンプルは、定常分布 ${\pi}$（つまり目的の分布）からのサンプルとみなせる。
という流れでサンプリングが実現できます。
これをうまく応用した手法群がMCMC法 (Markov Chain Monte Carlo methods) です。
MCMC でのマルコフ連鎖の使い方
サンプリングしたい分布に従うマルコフ連鎖を構築し、そこからのサンプルを採集する
3. 代表的な MCMC アルゴリズムと数式
MCMC 法にはさまざまなアルゴリズムが提案されていますが、代表的な2つの手法が メトロポリス・ヘイスティングス (Metropolis-Hastings) と ギブスサンプリング (Gibbs Sampling) です。
3.1 メトロポリス・ヘイスティングス法
メトロポリス・ヘイスティングス法は、提案分布 ${q(x’ \mid x)}$ を用いて次の状態 ${x’}$ をサンプリングし、それが「どれだけ妥当か」を受容・棄却のステップで調整しながらマルコフ連鎖を作る手法です。
これによって、最終的に目標分布 ${p(x)}$ に定常分布をもつ連鎖が作れます。
- 現在の状態を ${x}$ とする。
- 提案分布 ${q(x’ \mid x)}$ から次の状態候補 ${x’}$ をサンプリングする。
- 採択確率 ${\alpha(x’ \mid x)}$ を以下で計算する。
$${\alpha(x’ \mid x) = \min\!\Bigl\{\, 1,\; \frac{p(x’)\, q(x \mid x’)}{\,p(x)\, q(x’ \mid x)\!} \Bigr\}}$$
- 確率 ${\alpha(x’ \mid x)}$ で提案された状態 ${x’}$ を採択し、そうでなければ ${x}$ に留まる。
- 次のステップへ進む。
このように「提案分布」と「採択確率」の2段階で連鎖を作るため、高次元空間でも柔軟に適用できるメリットがあります。
一方で、採択率を調整しないと効率が下がる場合があるため、実務ではチューニングがしばしば大きな課題です。
3.2 ギブスサンプリング
ギブスサンプリングは、多変量確率分布を扱うときに特に有効な手法です。
複数の変数 ${(X_1, \dots, X_d)}$があるとき、それぞれの条件付き分布から順番にサンプリングを行います。
- 初期値 ${(x_1, \dots, x_d)}$をなんらかの方法で与える。
- 各ステップで、たとえば${X_1}$ から順に次のように更新する
$${X_1 \sim P(X_1 \mid X_2 = x_2, \dots, X_d = x_d), ⋮ \vdots⋮ X_d \sim P(X_d \mid X_1 = x_1′, \dots, X_{d-1} = x_{d-1}’)}$$
変数が一斉にではなく、1つずつ更新される点が重要です
- この手順を繰り返すと、やがて連鎖は目標とする多変量分布に従うようになります（エルゴード的である場合）。
ギブスサンプリングは「条件付き分布をサンプリングしやすい」という場面では実装が簡単で効率も高いですが、条件付き分布が複雑でサンプリングが難しいときは適用しにくい場合もあります。
4. ベイズ推論への応用
マルコフ連鎖が MCMC によって頻繁に使われる場面のひとつがベイズ推論です。
例えば、ある観測データ ${D}$ に対して、
- 事前分布 ${p(\theta)}$
- 尤度関数 ${p(D \mid \theta)}$
から計算される事後分布${p(\theta \mid D)}$ を得たい場面を考えましょう。
通常、
$${p(\theta \mid D) = \frac{p(D \mid \theta)\, p(\theta)}{\int p(D \mid \theta)\, p(\theta)\, d\theta}}$$
の分母（正規化定数）を厳密に計算するのは困難です。
しかし、メトロポリス・ヘイスティングス法などであれば、定数を必要とせず、比 だけで受容確率が決まるため、
$${\alpha(\theta’ \mid \theta) = \min\!\Bigl\{\,1,\; \frac{p(D \mid \theta’)\,p(\theta’)}{p(D \mid \theta)\,p(\theta)} \times \frac{q(\theta \mid \theta’)}{q(\theta’ \mid \theta)} \Bigr\}}$$
の形を直接計算できます。
そして、長期的に見ればサンプリングされた ${\theta}$ が事後分布に従うため、事後分布の期待値や分散、信用区間（credible interval）の推定などが可能になります。
6. サンプルコード: MH法による2次元正規分布サンプリング
ここでは、2次元の正規分布からのサンプリングをメトロポリス・ヘイスティングス法で実装する例を示します。コードは Python を想定し、NumPy を使います。
注意: 実務ではさらに多次元かつ複雑な分布に対して行うケースが多いため、チューニング（提案分布の選び方など）やバーンイン期間の設定が重要です。
import numpy as np
import matplotlib.pyplot as plt
# 目標分布: 2次元正規分布 N(mu, Sigma)
# muは2次元ベクトル、Sigmaは2x2正定値行列
mu = np.array([0.0, 0.0])
Sigma = np.array([[1.0, 0.5],
[0.5, 1.0]])
inv_Sigma = np.linalg.inv(Sigma)
det_Sigma = np.linalg.det(Sigma)
def target_pdf(x):
"""
目標分布の密度 p(x) を計算
x: np.array([x1, x2])
"""
diff = x - mu
exponent = -0.5 * diff @ inv_Sigma @ diff
coeff = 1.0 / (2.0 * np.pi * np.sqrt(det_Sigma))
return coeff * np.exp(exponent)
def proposal_sampler(x, scale=0.5):
"""
提案分布 q(x'|x) = N(x, scale^2 * I)
"""
return x + scale * np.random.randn(len(x))
# ここで関数
def metropolis_hastings(n_samples=5000, scale=0.5):
# 初期値を適当にセット
x_current = np.array([0.0, 0.0])
samples = []
for _ in range(n_samples):
# 提案分布からx_proposalをサンプリング
x_proposal = proposal_sampler(x_current, scale=scale)
# 採択確率 alpha(x'|x)
p_current = target_pdf(x_current)
p_proposal = target_pdf(x_proposal)
# 今回の提案分布は対称分布 (q(x'|x) = q(x|x'))
# なので ratio は p_proposal / p_current
ratio = p_proposal / p_current
# 受容or棄却
if np.random.rand() < min(1.0, ratio):
x_current = x_proposal
samples.append(x_current)
return np.array(samples)
# 実行
samples = metropolis_hastings(n_samples=10000, scale=0.5)
# 可視化
plt.figure(figsize=(6, 6))
plt.scatter(samples[:, 0], samples[:, 1], s=5, alpha=0.3)
plt.xlabel("x1")
plt.ylabel("x2")
plt.title("Metropolis-Hastings sampling from 2D Normal")
plt.axis("equal")
plt.show()
この例では、提案分布 ${q(x’ \mid x)}$ を「現在の点 ${x}$ を平均とする等方的ガウス」に設定し、対称性から採択確率の式が単純化されています。
実際には問題に応じて、提案分布の分散や形状を慎重に選び、バランスのとれたサンプリング（受容率が極端に低くも高くもならないように）を目指します。
- 多次元化や非対称な提案分布への拡張
- バーンインの除去
- 事後サンプルの依存度（サンプリング間の相関）を下げるための「サブサンプリング」
などなど、実務では更に工夫が行われます。
