MCMC（マルコフ連鎖モンテカルロ法）とは？仕組み・MH法とHMCの違い・収束判定をPythonで解説

MCMC（Markov Chain Monte Carlo、マルコフ連鎖モンテカルロ法）は、式は書けても直接は計算できない確率分布から、乱数で標本を引いてくる手法です。「目標の分布に落ち着くように設計したマルコフ連鎖」を長く走らせ、その足跡を標本として使います。主な用途はベイズ推定で、事後分布の平均や信用区間を、積分を解かずに標本の集計で求めます。

この記事では、MCMCがなぜ必要になるのか、代表的なアルゴリズム（メトロポリス・ヘイスティングス法、ギブスサンプリング、HMC/NUTS）の違い、収束判定の基準を、実際に動かしたPythonコードの出力とあわせて説明します。なお、日本産婦人科医会の「母と子のメンタルヘルスケア（MCMC）」は同じ略称の別の取り組みで、この記事では扱いません。

まとめ：MCMCの要点と使うときの判断基準

• MCMCは、正規化定数（ベイズ推定なら周辺尤度）が計算できない分布でも、密度の比さえ計算できれば標本を引ける
• 標本は互いに相関するため、引いた個数ではなく有効サンプルサイズ（ESS）で精度を見る
• 収束判定は、初期値をばらした複数チェーンで R-hat が1.01未満かを確認する（Vehtari et al. 2021 の推奨）
• 連続パラメータで勾配を計算できる実務のベイズモデルでは、Stan・PyMC・NumPyroで利用できる NUTS（HMCの自動調整版）から始めるのが基本
• 事後分布が解析的に求まるモデルや、精度より速度が要る大規模データでは、MCMC以外の手段を選ぶ

以下、仕組み、アルゴリズム、Python実装、収束判定、ライブラリの順に説明します。

MCMCの仕組み｜正規化定数が分からない分布から標本を得る考え方

ベイズ推定で事後分布を直接計算できない理由

ベイズの定理では、事後分布は「尤度×事前分布÷周辺尤度」で表されます。分子はモデルを書けば計算できますが、分母の周辺尤度はパラメータ全体にわたる積分です。共役事前分布のように式で解ける組み合わせは限られ、格子点で数値積分する方法は、1次元あたり100点でも10次元で100の10乗点が必要になり、次元が増えるほど計算量が指数的に膨らみます。

MCMCはこの分母を避けて通ります。後述する受容確率は「提案先と現在地の密度の比」だけで決まり、比を取ると分母は約分されて消えるためです。ベイズの定理そのものの意味はベイズの定理の公式と具体例で扱っています。

目標分布を定常分布にする遷移ルールの設計

マルコフ連鎖は、次の状態が現在の状態だけで決まる確率過程です。条件（既約性・非周期性など）を満たす連鎖は、どこから出発しても長く走らせると同じ分布に落ち着きます。この行き着く先が定常分布です。

MCMCはこの性質を逆向きに使います。先に「定常分布が目標分布（事後分布）になる」ように遷移のルールを作り、連鎖を走らせて、落ち着いた後の状態を標本として集めます。遷移ルールの設計に使われる代表的な十分条件が詳細釣り合い条件で、任意の2状態 x、y について「x にいる確率×x から y へ移る確率」と「y にいる確率×y から x へ移る確率」が等しいことを求めます。これを満たせば目標分布は定常分布になります。遷移行列と定常分布の計算例はマルコフ連鎖の仕組みと遷移行列・定常分布にまとめています。

モンテカルロ法・マルコフ連鎖との違い｜MCMCの標本間の相関

モンテカルロ法は、乱数で標本を作り、その平均で期待値や積分を近似する手法全般を指します。円周率を点の打ち込みで求める例が代表で、標本同士は独立です。MCMCはモンテカルロ法の一種ですが、標本をマルコフ連鎖で1つずつ作るため、隣り合う標本が似た値になります。

正の自己相関が強い場合、同じ個数の独立な標本より推定精度が低くなります。その影響を表すのが後述のESSです。ただし、負の自己相関によってESSが実際の標本数を上回る場合もあります。独立な標本を作る側の手法はモンテカルロ法の仕組みと円周率・リスク評価での使い方、棄却サンプリングの基本概念と応用例で扱っています。棄却サンプリングは、高次元で目標密度に近い提案分布と適切な上界を用意できないと受容率が低下します。そのような事後分布ではMCMCが選択肢になります。

代表的なMCMCアルゴリズム｜MH法・ギブスサンプリング・HMC/NUTSの使い分け

選び方の軸は「勾配が計算できるか」と「条件付き分布が既知か」の2つです。パラメータが連続で自動微分が使えるならNUTS、離散パラメータを含むならギブスやMHを組み合わせる、という順で考えます。

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

MH法は、現在地 x から提案分布 q で候補 x’ を作り、確率 min(1, π(x’)q(x|x’) / π(x)q(x’|x)) で移動し、却下なら x に留まります。π は目標分布の密度です。提案分布が対称（正規分布によるランダムウォークなど）なら q が約分され、密度の比だけで決まるメトロポリス法になります。

効率を左右するのは提案分布の幅です。幅が狭いと受容率は高くても少しずつしか動かず、広いと候補がほとんど却下されます。Roberts, Gelman, Gilks（1997）は、一定の滑らかさなどの条件を満たす同一の1次元密度の積を目標分布とし、正規ランダムウォーク提案を用いる高次元極限で、漸近的に最適な受容率が約0.234になることを示しました。次章の実測では、この幅の違いがESSに直接表れます。

ギブスサンプリングが使える条件

ギブスサンプリングは、パラメータを1つずつ、他を固定した条件付き分布から直接引く方法です。提案の却下が起きないため調整の手間がなく、共役事前分布を使った階層モデルや混合ガウスモデルのように条件付き分布が標準的な分布になる場合に向きます。一方、パラメータ同士の相関が強いと、1軸ずつしか動けないため連鎖がほとんど進まなくなります。BUGSやJAGSがこの方式を中心にしたソフトウェアです。

HMCとNUTSが既定になっている理由

HMC（ハミルトニアンモンテカルロ法）は、密度の勾配を使って物理の運動のように遠くの候補を作るため、ランダムウォークより少ない反復で分布全体を動けます。ただしステップ幅と軌道の長さという2つの調整値に敏感です。NUTS（No-U-Turn Sampler）は軌道の長さを自動で決める改良版で、CmdStanの既定サンプラーもNUTSです。CmdStanの既定値は、ウォームアップ1,000回・本サンプリング1,000回・目標受容率（adapt delta）0.8・木の最大深さ10です。

PythonでMH法を実装する｜コイン投げの事後分布を解析解と照合

答えが分かっている問題で動かすと、MCMCの挙動を確かめられます。コインを10回投げて7回表が出たとき、表の確率 p の事後分布は、一様な事前分布 Beta(1,1) のもとで Beta(8,4) になり、事後平均は 8/12＝0.667 です。これをMH法で求めます。NumPy以外のライブラリは使いません（NumPyの基本はNumPyの読み方と基本的な使い方を参照）。

import numpy as np

# 観測: コインを10回投げて表が7回。事前分布は一様分布 Beta(1,1)
n, k = 10, 7

def log_post(p):
    if p <= 0 or p >= 1:
        return -np.inf                             # 範囲外は確率0
    return k * np.log(p) + (n - k) * np.log(1 - p)   # 正規化定数は不要

def metropolis(step, n_iter=20000, init=0.5, seed=0):
    rng = np.random.default_rng(seed)
    p, lp = init, log_post(init)
    chain, accepted = np.empty(n_iter), 0
    for i in range(n_iter):
        prop = p + rng.normal(0, step)             # 対称な提案分布（ランダムウォーク）
        lp_prop = log_post(prop)
        if np.log(rng.uniform()) < lp_prop - lp:   # 受容確率 min(1, π(prop)/π(p))
            p, lp = prop, lp_prop
            accepted += 1
        chain[i] = p
    return chain, accepted / n_iter

def split_rhat(chains):
    # 各チェーンを前半・後半に割り、チェーン間とチェーン内の分散を比べる
    half = chains.shape[1] // 2
    c = np.vstack([chains[:, :half], chains[:, half:2 * half]])
    m, L = c.shape
    W = c.var(axis=1, ddof=1).mean()
    B = L * c.mean(axis=1).var(ddof=1)
    return np.sqrt(((L - 1) / L * W + B / L) / W)

def ess(x):
    # 自己相関が負に転じる手前までを足し合わせる単純な推定
    x = x - x.mean()
    n = len(x)
    acf = np.correlate(x, x, mode="full")[n - 1:] / (x.var() * n)
    s = 0.0
    for t in range(1, n):
        if acf[t] < 0:
            break
        s += acf[t]
    return n / (1 + 2 * s)

for step in (0.01, 0.2, 2.0):
    chains, rates = [], []
    for seed, init in enumerate((0.1, 0.4, 0.6, 0.9)):   # 初期値をばらした4本
        ch, r = metropolis(step, init=init, seed=seed)
        chains.append(ch[2000:])                          # 先頭2,000回はバーンイン
        rates.append(r)
    chains = np.array(chains)
    print(f"step={step:<5} 受容率={np.mean(rates):.2f} "
          f"事後平均={chains.mean():.3f} R-hat={split_rhat(chains):.3f} "
          f"ESS={sum(ess(c) for c in chains):.0f}")

Python 3.9・NumPy 2.0.2 で実行した出力です（初期値0.1／0.4／0.6／0.9の4本、各2万回からバーンイン2,000回を捨てた計72,000標本）。

step=0.01  受容率=0.98 事後平均=0.672 R-hat=1.042 ESS=133
step=0.2   受容率=0.59 事後平均=0.667 R-hat=1.000 ESS=13601
step=2.0   受容率=0.08 事後平均=0.667 R-hat=1.002 ESS=3365

提案分布の幅による受容率とESSの変化

幅0.01は受容率98%で、一見うまく動いているように見えます。しかし、この簡易計算でのESS推定値は72,000標本に対して133にとどまり、R-hatも1.042で基準の1.01を超えています。受容率が高いことは、よく混ざっていることを意味しません。幅0.2では事後平均が解析解の0.667と一致し、ESSは幅0.01の約100倍になりました。同じ関数で幅0.2・初期値0.5・seed=42として102,000回走らせ、先頭2,000回を除いた10万標本の2.5・97.5パーセンタイルは [0.385, 0.890] でした。Beta(8,4) の理論上の95%信用区間（等裾）は [0.390, 0.891] で、ほぼ重なります。

このコードの R-hat と ESS は仕組みを示すための簡易版です。実務では、順位正規化した R-hat と bulk/tail ESS を計算する ArviZ や Stan の関数を使ってください。

MCMCの収束判定｜R-hat・ESS・トレースプロットの見方

R-hat（Gelman-Rubin統計量）の目安｜順位正規化版で1.01未満

R-hat は、初期値を変えた複数チェーンについて「チェーン間のばらつき」と「チェーン内のばらつき」を比べる指標です。全チェーンが同じ分布に落ち着いていれば1に近づきます。かつては1.1未満が目安とされましたが、Vehtari et al.（2021）は順位正規化した split R-hat を提案し、一般的な閾値として1.01を推奨しています（R-hatとESSの推奨基準の原論文）。従来の R-hat は、収束前でも1.1を下回る例があると同論文は指摘しています。Stanのリファレンスマニュアルもこの基準を紹介しています。1本でも前半・後半に分割すれば split R-hat は計算できますが、異なる初期値からの探索を比較するため、独立した4本以上のチェーンを走らせることが推奨されます。

ESS（有効サンプルサイズ）による推定精度の評価

ESSは、相関のある標本列が独立な標本何個分の情報を持つかを表します。事後平均のモンテカルロ誤差はおよそ「事後標準偏差÷√ESS」なので、ESSが100なら誤差は標準偏差の1割です。ArviZや Stan は分布の中心部を見る bulk-ESS と、裾（信用区間の端）を見る tail-ESS を分けて出力します。信用区間を報告するなら tail-ESS も確認します。

トレースプロットの診断｜帯の安定とチェーン間の重なり

トレースプロットは、横軸に反復回数、縦軸に標本の値を描いた図です。収束して混ざりのよいチェーンは、一定の帯の中を細かく上下する毛虫のような「ギザギザ」になり、複数チェーンの線が重なります。逆に、ゆっくり波打つ、途中で水準が変わる、チェーンごとに別の高さに張り付く、といった形は収束していないサインです。上の実測で幅0.01のチェーンは、100回離れた標本同士でも自己相関が0.71残っており（幅0.2では10回離れると0.03）、ゆっくり波打つ形になります。

初期値の影響を減らすため、先頭部分をバーンインとして除外します。Stan・PyMCのウォームアップやチューニングでは、これに加えてステップ幅などのサンプラー設定も調整します。何回捨てるかに固定の正解はなく、捨てた後の R-hat と ESS で足りたかを判断します。

MCMCライブラリの選び方｜PyMC・Stan・NumPyroの違い

版番号は2026年9月26日時点のPyMCのリリース履歴、CmdStan 2.40の公式リリース、NumPyroのリリース履歴、ArviZのリリース履歴で確認した値です。PyMCの pm.sample() は、既定で draws＝1,000、tune＝1,000、引数を省略した場合の並列数はCPU数と4の小さい方、チェーン数はその並列数と2の大きい方です。NUTSのサンプラーは、nutpie（0.16.10以上）がインストールされ、全変数をNUTSで引けるなどの条件がそろえば nutpie を自動で選び、それ以外は PyMC 自身の実装を使います。古い書籍やブログにある import pymc3 は、2022年6月のv4.0で PyMC に改称される前の書き方で、PyMC3 の最終版は3.11.6です。

RからStanを使う場合の環境はR言語の言語仕様と環境構築もあわせて参照してください。

MCMCを使うべきでない場面と、よくある失敗

次の場面では、MCMCより先に別の手段を検討します。

• 事後分布が解析的に求まる：ベータ二項モデルのような共役な組み合わせなら、上の例のように式で答えが出ます。MCMCは検算用にとどめます。
• データ件数が多く、精度より速度が優先：各反復で全データの尤度を計算するため時間がかかります。変分推論（ADVI など）は近似誤差の検証が必要ですが、モデルによっては計算時間を短縮でき、PyMC・NumPyro とも実装を持っています。
• パラメータが識別できないモデル：ラベルの入れ替えで同じ尤度になる混合モデルなどでは、チェーンごとに別の山へ張り付き、R-hat が下がりません。制約や再パラメータ化が先です。

実務で多い失敗は、NUTSの発散（divergent transitions）の警告を無視することです。PyMCのドキュメントは、ウォームアップ後に発散が残る場合は結果が信頼できない可能性があるとし、再パラメータ化か target_accept を上げることを勧めています。PyMC 6.3.2 の既定サンプラー（nutpie）では既定値が0.9で、公式チュートリアルは発散が出た例で0.99に上げています。階層モデルで分散パラメータが0付近に寄ると発生しやすく、非中心化パラメータ化で解消することが多い問題です。もう1つは、チェーン1本・R-hat未確認のまま平均だけを報告することで、上の幅0.01の例のように、値がもっともらしく見えても収束していない場合があります。

よくある質問

MCMCは何回回せばよいですか？

回数に一律の基準はありません。まず4チェーンそれぞれでウォームアップ1,000回と本サンプリング1,000回を実行し、R-hat が1.01未満か、ESS が十分か（Vehtari et al. 2021 は順位正規化した ESS が400を超えることを推奨）で判断します。足りなければ回数を増やすより先に、モデルの再パラメータ化を検討します。

MCMCと変分ベイズ法の違いは何ですか？

MCMCは十分長く回せば事後分布そのものに近づく標本を得る方法、変分ベイズ法は事後分布を扱いやすい分布族で近似し、その最適化問題を解く方法です。変分ベイズ法は速い代わりに、近似分布の形で決まる偏り（分散を小さく見積もりがちなど）が残ります。精度を優先するならMCMC、速度を優先するなら変分ベイズ法、という使い分けが基本です。

MCMCはGPUで速くなりますか？

NumPyro は JAX による自動微分と JIT コンパイルで GPU・TPU・CPU 上で動き、NUTS の木構築まで含めてコンパイルします。PyMC でも pm.sample(nuts_sampler="numpyro") のように JAX 系の実装へ切り替えられます。効果が出やすいのは、勾配計算が重い大規模モデルや、多数のチェーンをまとめて回す場合です。

隠れマルコフモデルや状態空間モデルとMCMCの関係は？

隠れマルコフモデルや状態空間モデルは「観測の裏に、マルコフ連鎖に従う見えない状態がある」というモデルで、MCMCは「そのモデルのパラメータをどう推定するか」という計算手法です。線形・正規の状態空間モデルならカルマンフィルタで尤度を計算し、その尤度を使ってパラメータの事後分布をMCMCで求める、という組み合わせがよく使われます。カルマンフィルタ自体はカルマンフィルタの仕組みと使い方で解説しています。

PyMC3とPyMCの違いは何ですか？

同じライブラリの旧名と現行名です。2022年6月のv4.0で PyMC に改称され、計算の土台も Theano から Aesara、現在は PyTensor に移りました。PyMC3 の最終版は3.11.6で、新規に使うなら pip install pymc で入る PyMC（2026年9月時点で6.3.2）を選びます。

関連記事

• マルコフ連鎖とは？仕組み・遷移行列・定常分布とAI実装での使いどころ
• モンテカルロ法とは？乱数で解く仕組みと円周率・リスク評価・強化学習での使い方【2026年版】
• ベイズの定理とは？公式の意味・具体例からナイーブベイズ実装まで
• エンジニア向け入門：棄却サンプリングとは何か？基本概念と応用例をステップバイステップで徹底解説
• カルマンフィルタとは？仕組みと使い方をわかりやすく解説【線形の基礎から応用まで】