メトロポリス・ヘイスティングス法(M-H法)の仕組みを解説して実装する

メトロポリス・ヘイスティングス法(Metropolis-Hastings algorithm; M-H法)は、MCMC(マルコフ連鎖モンテカルロ法)における最も汎用的なアルゴリズムの1つです。ギブスサンプリングの一般化と位置づけることもできます。

今回は、メトロポリス・ヘイスティングス法の仕組みについてわかりやすくまとめていきます。

本記事の内容

  • M-H法の基本的な考え方
  • 提案分布と受理確率
  • アルゴリズムの定式化
  • Python での実装と可視化

前提知識

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

M-H法の基本的な考え方

目的の分布 $p(\theta)$ から直接サンプリングできない場合、M-H法は以下の戦略を取ります。

  1. 現在の状態 $\theta^{(t)}$ から、提案分布 $q(\theta’|\theta^{(t)})$ に従って新しい候補 $\theta’$ を生成する
  2. 受理確率 $\alpha$ を計算して、候補を受理するか棄却するかを確率的に決定する
  3. 受理すれば $\theta^{(t+1)} = \theta’$、棄却すれば $\theta^{(t+1)} = \theta^{(t)}$

この手続きを繰り返すことで、マルコフ連鎖の定常分布が目的の分布 $p(\theta)$ に一致するようになります。

受理確率

M-H法の受理確率は次のように定義されます。

$$ \begin{equation} \alpha = \min\left(1, \frac{p(\theta’) q(\theta^{(t)}|\theta’)}{p(\theta^{(t)}) q(\theta’|\theta^{(t)})}\right) \end{equation} $$

詳細つり合い条件

受理確率がこの形になる理由は、詳細つり合い条件(detailed balance condition)を満たすためです。

$$ p(\theta) T(\theta’|\theta) = p(\theta’) T(\theta|\theta’) $$

ここで $T(\theta’|\theta) = q(\theta’|\theta) \cdot \alpha(\theta, \theta’)$ は遷移確率です。この条件を満たすマルコフ連鎖は、$p(\theta)$ を定常分布として持ちます。

対称な提案分布の場合(メトロポリス法)

提案分布が対称 $q(\theta’|\theta) = q(\theta|\theta’)$ の場合、受理確率は簡略化されます。

$$ \alpha = \min\left(1, \frac{p(\theta’)}{p(\theta^{(t)})}\right) $$

これが元のメトロポリス法です。確率が高い方向への移動は常に受理され、確率が低い方向への移動は確率的に受理されます。

アルゴリズム

  1. 初期値 $\theta^{(0)}$ を設定
  2. $t = 0, 1, 2, \dots$ について繰り返し: – 提案分布 $q(\theta’|\theta^{(t)})$ から候補 $\theta’$ を生成 – 受理確率 $\alpha$ を (1) 式で計算 – 一様乱数 $u \sim \text{Uniform}(0, 1)$ を生成 – $u \leq \alpha$ なら $\theta^{(t+1)} = \theta’$(受理) – $u > \alpha$ なら $\theta^{(t+1)} = \theta^{(t)}$(棄却)

Python での実装

多峰性の分布からサンプリングする例を実装しましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 目的の分布: 2つのガウス分布の混合
def target_pdf(x):
    return 0.4 * norm.pdf(x, -2, 0.8) + 0.6 * norm.pdf(x, 3, 1.2)

def metropolis_hastings(target, n_samples, proposal_std=1.0, x0=0.0, burn_in=1000):
    """メトロポリス・ヘイスティングス法"""
    samples = np.zeros(n_samples + burn_in)
    x = x0
    n_accepted = 0

    for t in range(n_samples + burn_in):
        # 提案分布(対称な正規分布)から候補生成
        x_proposal = x + np.random.normal(0, proposal_std)

        # 受理確率の計算
        alpha = min(1, target(x_proposal) / (target(x) + 1e-300))

        # 受理/棄却
        if np.random.uniform() <= alpha:
            x = x_proposal
            if t >= burn_in:
                n_accepted += 1

        samples[t] = x

    acceptance_rate = n_accepted / n_samples
    return samples[burn_in:], acceptance_rate

np.random.seed(42)
n_samples = 20000

# 異なる提案分布の幅で比較
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

proposal_stds = [0.1, 1.0, 10.0]
x_range = np.linspace(-6, 8, 500)
target_values = target_pdf(x_range)

for idx, std in enumerate(proposal_stds):
    samples, acc_rate = metropolis_hastings(target_pdf, n_samples, proposal_std=std)

    # 上段: サンプリングの軌跡
    ax_trace = axes[0, idx]
    ax_trace.plot(samples[:500], linewidth=0.5, alpha=0.7)
    ax_trace.set_xlabel('Iteration', fontsize=10)
    ax_trace.set_ylabel('x', fontsize=10)
    ax_trace.set_title(f'Trace (proposal_std={std}, acc={acc_rate:.2f})', fontsize=11)
    ax_trace.grid(True, alpha=0.3)

    # 下段: ヒストグラムと理論分布
    ax_hist = axes[1, idx]
    ax_hist.hist(samples, bins=80, density=True, alpha=0.5, color='blue', label='M-H samples')
    ax_hist.plot(x_range, target_values, 'r-', linewidth=2, label='Target')
    ax_hist.set_xlabel('x', fontsize=10)
    ax_hist.set_ylabel('Density', fontsize=10)
    ax_hist.set_title(f'proposal_std={std}', fontsize=11)
    ax_hist.legend(fontsize=9)
    ax_hist.grid(True, alpha=0.3)

plt.suptitle('Metropolis-Hastings Algorithm: Effect of Proposal Distribution Width', fontsize=14)
plt.tight_layout()
plt.show()

提案分布の幅が小さすぎると(左)探索が遅く受理率は高いが収束が遅い、適切な場合(中央)効率的にサンプリングでき、大きすぎると(右)多くの候補が棄却され効率が悪くなることが確認できます。

ギブスサンプリングとの比較

特性 M-H法 ギブスサンプリング
適用範囲 任意の分布 完全条件付き分布が既知
棄却 あり なし(常に受理)
提案分布 設計が必要 完全条件付き分布を使用
チューニング 提案分布の調整が必要 不要

まとめ

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

  • M-H法は提案分布から候補を生成し、受理確率に基づいて受理/棄却を行うMCMC法
  • 受理確率は詳細つり合い条件を満たすように設計されている
  • 提案分布の幅の調整が性能に大きく影響する(受理率20-50%が目安)
  • ギブスサンプリングよりも汎用的だが、チューニングが必要

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