収束診断と実装の注意｜シミュレーション

← シミュレーション一覧

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

📎 前提：メトロポリス・ヘイスティングス法　|　ベイズ：収束診断

要点（BLUF）

MCMC は「いつ収束したか」が自分では分からないので、収束診断が必須。出力が定常分布に達したかを統計的にチェックします。
R-hat（Gelman–Rubin）：複数連鎖の「連鎖間分散」と「連鎖内分散」を比べる。収束すると 1.0 に近づく（目安 < 1.01）。
有効標本サイズ（ESS）：自己相関を考慮した「独立サンプル換算数」。標準正規を狙う4連鎖で R-hat = 1.0002、ESS を実測します。

1. なぜ診断が要るか

MCMC サンプルには2つの厄介な性質があります。

過渡（バーンイン）：初期は定常分布に達していない（マルコフ連鎖と定常分布）。最初の一定区間は捨てる。
自己相関：連続するサンプルは相関する（メトロポリス・ヘイスティングス法）。 $n$ サンプルあっても独立 $n$ 個ぶんの情報はない。

しかも「収束したか」は連鎖を見ても確定できません（多峰分布の片方の山に嵌って気づかない、など）。そこで複数の診断を併用します。

2. R-hat（Gelman–Rubin 統計量）

異なる初期値から複数の連鎖を走らせ、(a) 連鎖間のばらつき $B$ と (b) 連鎖内のばらつき $W$ を比べます。収束していれば、どの連鎖も同じ定常分布を見ているので両者が一致し、比

\hat{R} = \sqrt{\frac{\widehat{\text{Var}}^+}{W}},\qquad \widehat{\text{Var}}^+ = \frac{L-1}{L}\,W + \frac{1}{L}\,B

が 1に近づく（ $L$ は連鎖長）。未収束だと連鎖がバラバラの場所にいて $B \gg W$ 、 $\hat{R} \gg 1$ になります。目安は $\hat{R} < 1.01$ 。

import numpy as np

# 乱数シードを固定
rng = np.random.default_rng(43)

def mh_chain(seed, n, step=2.0, x0=0.0):
    """標準正規を狙うランダムウォークMH"""
    r = np.random.default_rng(seed)
    x = x0; out = np.empty(n)
    for i in range(n):
        xp = x + r.normal(0, step)
        if np.log(r.random()) < (-0.5*xp**2) - (-0.5*x**2):
            x = xp
        out[i] = x
    return out

M, n = 4, 50_000
# 異なる初期値（散らばった出発点）から4連鎖
chains = np.array([mh_chain(50+m, n, x0=rng.normal(0, 3)) for m in range(M)])
burn = n // 2
ch = chains[:, burn:]                  # 後半だけ使う（バーンイン除去）

m, L = ch.shape
chain_means = ch.mean(axis=1)
chain_vars = ch.var(axis=1, ddof=1)
B = L * chain_means.var(ddof=1)        # 連鎖間分散
W = chain_vars.mean()                  # 連鎖内分散
var_hat = (1 - 1/L)*W + B/L
Rhat = np.sqrt(var_hat / W)
print(f"R-hat = {Rhat:.4f}  (収束なら ~1.0)")

出力：

R-hat = 1.0002  (収束なら ~1.0)

出力の意味：4本の連鎖を散らばった初期値から出したのに R-hat = 1.0002。連鎖間と連鎖内のばらつきがほぼ一致＝全連鎖が同じ標準正規を見ている＝収束した、と判断できます。これが 1.05 や 1.2 なら、連鎖がまだバラバラの場所にいて収束していない（バーンインを伸ばす・連鎖を長くする必要がある）。

3. 有効標本サイズ（ESS）

自己相関があると、 $L$ サンプルの情報量は独立 $L$ 個より少ない。有効標本サイズは

\text{ESS} = \frac{L}{\tau},\qquad \tau = 1 + 2\sum_{k=1}^{\infty}\rho_k

$\tau$ は積分自己相関時間（ $\rho_k$ はラグ $k$ の自己相関）。「実質的に独立サンプルが何個ぶんあるか」を表します。

import numpy as np

def mh_chain(seed, n, step=2.0, x0=0.0):
    r = np.random.default_rng(seed)
    x = x0; out = np.empty(n)
    for i in range(n):
        xp = x + r.normal(0, step)
        if np.log(r.random()) < (-0.5*xp**2) - (-0.5*x**2):
            x = xp
        out[i] = x
    return out

def ess(x):
    x = x - x.mean(); n = len(x)
    acf = np.correlate(x, x, mode='full')[n-1:] / (x.var() * np.arange(n, 0, -1))
    tau = 1.0
    for k in range(1, n):
        if acf[k] < 0: break            # 自己相関が負になったら打ち切り
        tau += 2*acf[k]
    return n/tau, tau

chain = mh_chain(60, 50_000)[25_000:]   # バーンイン後の25000サンプル
e, tau = ess(chain)
print(f"積分自己相関時間 tau = {tau:.1f}")
print(f"有効標本サイズ ESS (L=25000) = {e:.0f}")

出力：

積分自己相関時間 tau = 4.7
有効標本サイズ ESS (L=25000) = 5292

出力の意味：25000サンプルあっても、自己相関時間 $\tau=4.7$ のせいで実質的に独立なのは約5300個ぶん。つまり連続する約4.7サンプルが「ほぼ1個ぶん」の情報。歩幅 step を最適化すれば $\tau$ が下がり ESS が増えます。推定の精度（標準誤差）は $L$ ではなく ESS で決まる（ $\sigma/\sqrt{\text{ESS}}$ ）ので、ESS こそが MCMC の「実質サンプル数」です。

4. 実装の作法

複数連鎖を散らばった初期値から：R-hat には必須。1連鎖では収束を確認できません。
バーンイン除去：前半（または十分な区間）を捨てる。R-hat や trace を見て決める。
trace プロット：サンプルの時系列。複数連鎖が同じ帯に重なって「毛虫」状なら収束のサイン。
シン（間引き）は基本不要：自己相関を減らすための間引きは情報を捨てるだけで、現在は非推奨（ESS で評価すれば足りる）。
多峰分布に注意：モード間を行き来できないと、R-hat が良くても一部しか探索していないことがあります。

数式の直観的意味

R-hat は「みんな同じ景色を見ているか」のチェックです。違う場所から出発した登山者（連鎖）が、十分歩いた後に全員が同じ高度分布を報告すれば、その分布が真の地形（定常分布）。報告がバラついていれば、まだ各自が出発点近くの局所にいる証拠（ $B \gg W$ ）。ESS は「サンプルの水増し度」の補正です。自己相関した連続サンプルは「同じ情報の繰り返し」を含むので、額面の $L$ を $\tau$ で割って実質の独立数に直す。 $\tau$ は「前の状態を忘れるのにかかる歩数」で、これが MCMC の効率そのもの。R-hat で「収束したか」、ESS で「どれだけ精度があるか」——この2つで MCMC の信頼性が定量化されます。

⚠️ よくある誤解・落とし穴

「1連鎖で R-hat が計算できる」ではない：複数連鎖（散らばった初期値）が必須。
「R-hat ≈ 1 なら完全に収束」ではない：多峰でモードを見落としても R-hat が良く見えることがある。trace・複数初期値で補強。
「サンプル数が多ければ精度が高い」ではない：効くのは ESS。自己相関が強いと $L$ が大きくても ESS は小さい。
「間引き（thinning）で精度が上がる」ではない：情報を捨てるだけ。メモリ節約以外の利点は基本なし。
「バーンインは固定長でいい」ではない：収束の速さは問題依存。診断を見て決めます。
「収束診断は推論固有」ではない：これはサンプリング技法の数値挙動の問題。推論への応用はベイズへ（境界）。

対応シミュレーション参照

本文の R-hat（4連鎖・default_rng(43) と default_rng(50+m)）と ESS（default_rng(60)、 $\tau=4.7$ , ESS≈5300）。応用はベイズへ。