収束診断｜ベイズ統計

← ベイズ統計一覧

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ハミルトニアンモンテカルロとNUTS　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計）

要点（BLUF）

MCMC は有限回では「まだ定常に達していない（バーンイン）」「強い相関でサンプルが痩せている」「多峰に閉じ込められた」などが起こる。診断でサンプルが事後を代表しているか確かめます。
R-hat（Gelman–Rubin）：分散した初期値から複数チェーンを回し、チェーン間分散とチェーン内分散を比較。 $\hat R\approx1$ で収束、 $1.1$ 以上で未収束のサイン。
ESS（有効サンプルサイズ）：相関を割り引いた「実質の独立標本数」。実装で、よく混合した連鎖（ $\hat R\approx1.00$ ）と二峰に閉じ込められた連鎖（ $\hat R\approx9$ ）を区別します。

1. なぜ診断するのか

MCMC は理論上は無限ステップで事後 $\pi$ に収束しますが、有限回では3つの問題が起こります。

バーンイン：初期値の影響が残る最初の数百〜数千ステップ。捨てる。
自己相関：連続する標本は似ている（前の点から少しずらすだけ）。相関が強いほど実質の情報が少ない。
多峰・閉じ込め：山が複数あると、1本の連鎖が片方に閉じ込められ、事後全体を代表しない。

これらを見抜く2大指標が R-hat と ESS です。

2・3. R-hat と ESS を実装し、良い連鎖と悪い連鎖を区別する

R-hat は、分散した初期値から $m$ 本のチェーン（各長さ $n$ ）を回し、

W=\frac1m\sum_{j}s_j^2\ (\text{チェーン内分散}),\quad B=\frac{n}{m-1}\sum_j(\bar\theta_j-\bar\theta)^2\ (\text{チェーン間分散})

から $\widehat{\mathrm{Var}}=\frac{n-1}{n}W+\frac1n B$ 、 $\hat R=\sqrt{\widehat{\mathrm{Var}}/W}$ 。チェーン間とチェーン内のばらつきが釣り合えば（収束していれば） $\hat R\to1$ 。ESS は自己相関 $\rho_k$ を使って $\mathrm{ESS}\approx N/(1+2\sum_{k\ge1}\rho_k)$ 。両方を実装して、2つのケースで比べます。

import numpy as np

def rhat(chains):                       # chains: 形 (m, n)
    m, n = chains.shape
    W = chains.var(axis=1, ddof=1).mean()           # チェーン内分散
    B = n * chains.mean(axis=1).var(ddof=1)         # チェーン間分散
    var_hat = (n - 1)/n * W + B/n
    return np.sqrt(var_hat / W)

def ess(chains):                        # 自己相関ベースの簡易ESS
    x = (chains - chains.mean()).ravel(); N = len(x)
    acf = np.correlate(x, x, mode='full')[N-1:]; acf /= acf[0]
    s = 1.0
    for k in range(1, N):
        if acf[k] < 0: break                         # 初期正値まで和を取る
        s += 2*acf[k]
    return N / s

def mh_chains(logp, inits, step, n=5000, seed=0):
    rng = np.random.default_rng(seed); out = []
    for x0 in inits:
        x = x0; lp = logp(x); ch = np.empty(n)
        for i in range(n):
            cand = x + rng.normal(0, step); lc = logp(cand)
            if np.log(rng.uniform()) < lc - lp: x, lp = cand, lc
            ch[i] = x
        out.append(ch)
    return np.array(out)

# 良いケース：標準正規・分散した初期・適切な step → よく混合
good = mh_chains(lambda x: -0.5*x**2, inits=[-3,-1,1,3], step=2.0)[:, 1000:]
print(f"良いケース（標準正規）:  R-hat={rhat(good):.3f}  ESS={ess(good):.0f}/{good.size}")

# 悪いケース：二峰 0.5N(-4,.5²)+0.5N(4,.5²)・各チェーンが片方の山に閉じ込め
def logp_bimodal(x):
    return np.logaddexp(-0.5*((x+4)/0.5)**2, -0.5*((x-4)/0.5)**2)
bad = mh_chains(logp_bimodal, inits=[-4,-4,4,4], step=0.5)[:, 1000:]
print(f"悪いケース（二峰・閉じ込め）: R-hat={rhat(bad):.3f}  ESS={ess(bad):.0f}/{bad.size}")
print(f"  各チェーンの平均={bad.mean(axis=1).round(2)}（±4 に割れている＝未収束）")

出力：

良いケース（標準正規）:  R-hat=1.000  ESS=3326/16000
悪いケース（二峰・閉じ込め）: R-hat=9.305  ESS=3/16000
  各チェーンの平均=[-3.94 -3.96  4.01  4.  ]（±4 に割れている＝未収束）

出力の意味：良いケースは $\hat R=1.000$ （チェーン間と内のばらつきが一致＝みな同じ分布を見ている）、ESS も $16000$ 中 $3326$ と健全。悪いケースは $\hat R=9.305$ で一目で未収束——4本のチェーンが $-4$ と $+4$ の山に2本ずつ閉じ込められ、チェーン間分散 $B$ が爆発したためです。ESS はわずか $3$ （実質3個ぶんの情報しかない）。分散した初期値で複数チェーンを回し $\hat R$ を見るだけで、この閉じ込めを検出できます。1本だけ見ていたら「片方の山できれいに収束した」と誤認したでしょう。

4. 実務の目安と対処

指標	目安	外れたときの対処
$\hat R$	$<1.01$ （厳しめ）。 $>1.1$ は明確に未収束	チェーン数・反復数を増やす、初期値を分散、再パラメータ化
ESS	関心量ごとに数百以上（区間端なら多めに）	反復を増やす、提案・サンプラー改善（HMC/NUTS）、間引きは基本不要
受容率	RW-MH で 0.2〜0.5、HMC で 0.6〜0.9	step（提案幅・ $\varepsilon$ ）を調整（メトロポリスヘイスティングス）
発散（divergence）	0 が理想（NUTS）	再パラメータ化（階層モデルの実例と再パラメータ化）、step 小さく

トレースプロットは最初の目視点検——よく混合した連鎖は水平な帯（毛虫状）、未収束はドリフトや帯の分離が見えます。良いケースと悪いケースのトレースを並べます。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用

def mh_chains(logp, inits, step, n=5000, seed=0):
    rng = np.random.default_rng(seed); out = []
    for x0 in inits:
        x = x0; lp = logp(x); ch = np.empty(n)
        for i in range(n):
            cand = x + rng.normal(0, step); lc = logp(cand)
            if np.log(rng.uniform()) < lc - lp: x, lp = cand, lc
            ch[i] = x
        out.append(ch)
    return np.array(out)

good = mh_chains(lambda x: -0.5*x**2, [-3,-1,1,3], 2.0)
bad  = mh_chains(lambda x: np.logaddexp(-0.5*((x+4)/0.5)**2, -0.5*((x-4)/0.5)**2), [-4,-4,4,4], 0.5)
fig, ax = plt.subplots(1, 2, figsize=(11, 4), sharey=False)
for c in good: ax[0].plot(c[:1500], lw=0.6)
ax[0].set_title("良い：4本が重なり帯状（R-hat≈1.00）"); ax[0].set_xlabel("ステップ")
for c in bad: ax[1].plot(c[:1500], lw=0.6)
ax[1].set_title("悪い：±4 の山に分離（R-hat≈9.3）"); ax[1].set_xlabel("ステップ")
plt.tight_layout(); plt.show()

グラフの意味：左は4本のチェーンが同じ帯に重なり合い、よく混合（収束）。右は2本が下の山（ $-4$ ）、2本が上の山（ $+4$ ）に貼り付いたまま行き来できず、帯がくっきり分離します。これが $\hat R=9.3$ の正体で、トレースを見れば一目で分かります。

⚠️ よくある誤解

「1本のチェーンが帯状なら収束」ではない：多峰では片方に閉じ込められても帯状に見えます。複数チェーン＋ $\hat R$ が必須（§3）。
「標本数が多ければ ESS も多い」ではない：相関が強いと $N$ が大きくても ESS は小さい（悪いケースは $16000$ 標本で ESS $3$ ）。
「 $\hat R=1$ なら絶対に正しい」ではない：未踏の山があれば全チェーンがそこを見逃して $\hat R\approx1$ になりえます。事前知識・事後予測チェック（事後予測チェック）と併用。
「間引き（thinning）で相関を消すべき」ではない：原則は全標本を使うほうが推定は良い。間引きは保存容量の都合のときだけ。

まとめ（Phase 4）

第4章では、共役が崩れた事後からサンプルを生成する道具を揃えました——なぜサンプリングか（なぜサンプリングか）、規格化なしで引く MH（メトロポリスヘイスティングス）、条件付きで引くギブス（ギブスサンプリング）、勾配で滑空する HMC/NUTS（ハミルトニアンモンテカルロとNUTS）、そして出力を信じてよいか確かめる収束診断（本ノート）。得たサンプルは第3章の要約（事後分布の要約）にそのまま渡せます。次章では、この計算力を前提に、グループ構造を持つデータをまとめて扱う階層ベイズモデルへ進みます。