MCMC入門と実践:マルコフ連鎖モンテカルロの理論・アルゴリズム・実装ガイド

概要 — MCMCとは何か

マルコフ連鎖モンテカルロ(MCMC: Markov chain Monte Carlo)は、複雑な確率分布からのサンプリングを行うための一連の手法です。ベイズ推論や高次元パラメータ空間の探索、統計物理学、機械学習などで広く用いられます。MCMCは、目標分布π(x)(通常は事後分布)を直接サンプリングできない場合に、漸近的にπを標本化するマルコフ連鎖を構築することで解決します。

基礎理論

MCMCのコアはマルコフ連鎖の性質です。状態空間上で移動する連鎖が持つ遷移確率P(x->x')を設計し、
長く走らせることで連鎖の分布が一意の不変分布(stationary distribution)に収束すれば、その不変分布からのサンプルを得られます。重要な概念は以下です。

  • 不変分布: πはπP=πを満たす分布。
  • 詳細釣り合い(detailed balance): π(x)P(x->x')=π(x')P(x'->x)。多くのMCMCアルゴリズムはこれを満たすように設計される。
  • エルゴード性: 初期分布に関係なく時間平均が期待値へ収束する性質。漸近的にπからの代表的な標本が得られるために必要。

代表的アルゴリズム

ここでは主要なMCMCアルゴリズムの仕組みと特徴を解説します。

1. メトロポリス法(Metropolis)とメトロポリス・ヘイスティングス(Metropolis-Hastings)

メトロポリス法(Metropolis, 1953)は、提案分布q(x'|x)に従って候補x'を生成し、受容確率α=min(1, π(x')/π(x))で受容する単純な仕組みです。提案分布が対称(例: 正規分布)ならこの式になります。これを一般化したのがMetropolis-Hastings(Hastings, 1970)で、非対称提案を扱うため受容確率は
α=min(1, [π(x') q(x|x')] / [π(x) q(x'|x)])
となります。利点は単純さと汎用性。欠点は高次元では混合が遅くなることです。

2. ギブスサンプリング(Gibbs Sampling)

ギブスは複数次元パラメータを変数ごとの条件付き分布から順にサンプリングする手法です。各成分のフル条件付き分布が解析的にサンプル可能であれば効率的に動作します。データ拡張や階層モデルで非常に有用ですが、各条件付きが得られない場合は適用が難しい。

3. ハミルトニアン(Hamiltonian)モンテカルロ(HMC)

HMCは確率微分方程式に基づき、勾配情報を利用して高次元空間を効率的に探索します。潜在的エネルギーU(x)=-log π(x)と運動量pを導入し、ハミルトン方程式を近似解く(通常はリープフロッグ積分)ことで提案点を生成します。受容確率は数値誤差のためにMetropolis判定で補正します。HMCは連続かつ微分可能なモデルで非常に効果を発揮し、サンプルの自己相関を大幅に減らします。ステップサイズεとステップ数Lの調整が必要で、NUTS(No-U-Turn Sampler)はこのLを自動化する手法です(Hoffman & Gelman, 2014)。

4. アンサンブル手法(例: emcee のアフィン不変サンプラー)

複数のウォーカー(並列チェーン)を使って相互に提案を生成することで、スケーリングや形状の歪みに強い手法です。計算の面で並列化しやすく、分布が細長い場合などに有利です。

実装上の重要事項とチューニング

  • バーンイン(burn-in): 初期値の影響を取り除くため、初期サンプルを破棄します。何サンプル破棄するかは診断で判断します。
  • 自動調整(adaptation): ウォームアップ期間中にステップサイズや質量行列を適応的に調整することが一般的(例: Stanのadapt_deltaや自動質量行列推定)。
  • 受容率の目安: Metropolis系では受容率が極端に高い(小刻み)・低い(拒否多)と効率が悪くなる。経験則として次元に依存するが、簡易な目安は0.2〜0.5程度。ただしHMC系では受容率を高めに保つ傾向がある。
  • 初期値と複数チェーン: 複数の初期化で独立にチェーンを走らせ、収束性を評価します。局所モードに捕捉されるリスクを減らします。

収束診断とサンプルの質

MCMCの正当性を担保するためにいくつかの診断が重要です。

  • Gelman-Rubin診断(R-hat): 複数チェーンの間の分散とチェーン内分散を比較し、R-hat≈1に近ければ収束の目安。一般に1.1未満が推奨される。
  • 有効サンプルサイズ(ESS: Effective Sample Size): 自己相関を考慮した実質的な独立サンプル数。ESSが十分でないと推定の誤差が大きくなる。
  • トレースプロットと自己相関プロット: チェーンの混合や周期性、収束の有無を可視化する基本手段。
  • 事後予測チェック(Posterior Predictive Check): モデルが観測データをどれだけ再現するかを確認し、モデル適合性を評価する。

計算資源と並列化

MCMCは計算負荷が高いことが多く、実運用では効率化が重要です。並列化の典型は独立チェーンの同時計算(複数コア利用)です。HMCのように勾配計算が多いアルゴリズムでは、GPUや自動微分(AD)を利用すると大幅に高速化できます(StanはAD、TensorFlow ProbabilityやPyTorchベースの実装はGPUを活用)。また、分散MCMCや並列化可能なサブサンプリング手法も研究されていますが、バイアス管理には注意が必要です。

ソフトウェアと実務での選択肢

  • Stan(https://mc-stan.org): HMC(NUTS)を実装し、AD・自動適応機構を備えた高性能ソフト。統計コミュニティで広く使用。
  • PyMC(https://www.pymc.io): Pythonベースの確率的プログラミングツール。NUTSや変分推論をサポート。
  • emcee(https://emcee.readthedocs.io): 天文学でよく使われるアフィン不変アンサンブルサンプラー。
  • JAGS(https://mcmc-jags.sourceforge.net)/BUGS: 階層ベイズモデル向けの古典的ツール。

現場で気をつける点(IT視点)

  • スケーリング: 高次元では単純なランダムウォークは実用的でない。勾配法(HMC)や事前のパラメータ変換が有効。
  • 数値安定性: 対数確率を扱い、オーバーフロー・アンダーフローを防ぐ。対数和exp_log_sumなどのテクニックを使う。
  • デバッグ: モデルが適切に実装されているかを検証するため、既知解や簡易ケースでテストする。分布のモードや形状を確認。
  • 監査と再現性: シード管理、チェーンの保存、ウォームアップ設定の記録を行い再現可能にする。

MCMCの限界と代替手法

MCMCは汎用性が高い反面、計算コストや収束診断の難しさがあります。大規模データやリアルタイム推定では変分推論(VI: Variational Inference)や近似ベイズ法が現実的な選択肢になることがあります。ただしVIは近似誤差の評価が難しく、MCMCはゴールドスタンダードとして残る場面が多いです。

まとめと実践チェックリスト

  • まずモデルと目標分布を明確化する(解析的に条件付きが取れるか、勾配が使えるか)。
  • 適切なアルゴリズム選択(Metropolis系、Gibbs、HMC/NUTS、アンサンブル)を行う。
  • 複数チェーン、ウォームアップ、自動調整を行い、R-hat・ESS・トレースプロットで診断する。
  • 実装は既存ライブラリ(Stan、PyMC、emcee等)活用を検討。必要ならGPU/ADで高速化。
  • 事後予測チェックでモデルの妥当性を検証する。

参考文献