統計誤差の解析：自己相関・バッチ平均・ジャックナイフとブートストラップ

0. この記事の要点

乱数を使った計算で平均を求めるとき、結果には「どれくらい信用できるか」を表す誤差棒を付けます。ところが、前の結果から次の結果を作る方法（Markov 連鎖）では、隣り合う標本が似ているので、独立な標本の公式をそのまま使うと誤差棒が小さすぎます。この記事では、似ている標本があっても誤差棒を正しく付ける方法と、平均以外の量（比や分散など）の誤差をデータの作り直しで見積もる方法（ジャックナイフとブートストラップ）を、どこまで信用できてどこで壊れるかまで含めて調べます。読み終えると、誤差棒の付け方を選び、その偏りと限界を説明できるようになります。

• 相関のある標本 $N$ 個の平均の分散は、独立な標本 $N/(2\tau_{\mathrm{int}})$ 個ぶんです。相関を無視した標準誤差は $\sqrt{2\tau_{\mathrm{int}}}$ 倍だけ小さく出ます。この章の仕事は、$\tau_{\mathrm{int}}$（または長期分散 $\sigma_\infty^2$）をデータから推定する方法の偏りと分散を評価することです。
• 有限・可逆な連鎖では、自己共分散は $\Gamma(t) = \sum_j c_j\lambda_j^{t}$（$c_j \ge 0$）と書け、隣り合う 2 項の和 $\Gamma(2k) + \Gamma(2k+1)$ は非負・単調非増加・凸です（定理 3.1）。推定した対の和が負になったら雑音しか残っていないと判断して和を打ち切れます。
• 窓 $W$ で打ち切ると、長期分散の推定量の期待値を $2\Gamma(0)$ で割った量の偏りは $-\sum_{t > W}\rho(t) + O(W/N)$ です（命題 4.1）。分散は $W/N$ に比例して増えます（Madras–Sokal の近似公式）。バッチ平均の偏りはバッチ長 $b$ について $O(1/b)$ です。バッチ数を固定したままでは、バッチ長をいくら大きくしても推定量は一致するとは限りません（独立な正規標本をバッチ 2 個に分けた場合に、一致しないことを示します。定理 5.2）。
• ジャックナイフ（データを 1 個ずつ除いて同じ計算を繰り返し、ばらつきから分散を見積もる方法）の分散推定量は、標本平均の滑らかな関数については delta 法の分散に概収束します（定理 6.3）。中央値のような滑らかでない統計量では一致しません。
• ブートストラップ（手元のデータから重複を許して取り直し、分布を再現する方法）は標本平均の分布を $O(N^{-1/2})$ の精度で再現しますが（定理 7.3）、一様分布の最大値では Kolmogorov 距離が $1 - 1/e$ より小さくなりません（命題 7.5）。
• 初期分布による偏りは $O(1/(\gamma_* N))$ で統計誤差 $O(N^{-1/2})$ より高次です（命題 8.2）。複数の連鎖から作る $\hat R$ は、連鎖がエルゴード的なら $1$ に収束しますが、$\hat R \approx 1$ は収束の十分条件ではありません（命題 8.4）。

この記事で使う事前知識

• 標準誤差: 独立な標本 $N$ 個の平均がどれくらいばらつくかを表す量で、$\sigma/\sqrt N$ です。→ 標準誤差
• 共分散: 2 つの確率変数が一緒に大きく（小さく）なりやすいかを表す量です。→ 分散・共分散
• 定常分布と可逆性: 1 ステップ進めても変わらない分布が定常分布で、行って戻る流れが釣り合っている性質が可逆性です。→ 定常分布・可逆性
• スペクトルギャップ: 連鎖がどれだけ速く定常分布に近づくかを決める、遷移行列の固有値の隙間です。→ スペクトルギャップ
• 概収束と分布収束: 確率 1 で値が近づく収束と、分布の形が近づく収束です。→ 概収束、分布収束

1. 動機

Monte Carlo 積分の章では、独立な標本 $N$ 個の平均の標準誤差が $\sigma/\sqrt N$ で、標本分散から推定できることを示しました。Markov 連鎖 Monte Carlo の標本は独立ではありません。隣り合う標本が似ているので、$N$ 個の標本がもつ情報は独立な $N$ 個より少なく、独立と見なして計算した標準誤差は小さすぎます。誤差棒が小さすぎる結果は、誤差棒のない結果より害があります。読者が「理論値と $5\sigma$ ずれている」と誤って結論するからです。

もう 1 つの問題は、報告したい量が標本平均そのものでないことです。比熱は分散、質量は相関関数の対数の差、再重みづけした期待値は 2 つの平均の比です。平均の非線形な関数には $O(1/N)$ の偏りがあり、誤差の伝播の式を量ごとに手で導くのは間違いのもとです。ジャックナイフとブートストラップは、データを作り直して同じ計算を繰り返すだけで誤差を出す方法です。便利な反面、どの統計量で正しく、どこで壊れるかを知らずに使うと危険です。

この章では、相関の扱い（§3〜§5）、非線形な量の誤差（§6・§7）、連鎖の初期の扱いと収束の診断（§8）の順に、各推定量の偏りと分散を定理として評価します。自己相関時間の定義、相関のある標本平均の分散、ビニング、ジャックナイフによる偏りの除去は格子上の Monte Carlo 法の章が正本です。この章では必要な定義と主張を述べ直したうえで、そこに無い結果を証明します。

2. 準備

2.1. 定常な列と積分自己相関時間

定義 2.1（弱定常列・自己共分散・長期分散）↗

2 乗可積分な実確率変数の列 $A_1, A_2, \ldots$ が弱定常であるとは、$\mathbb{E}[A_s] = \mu$ が $s$ によらず、$\mathrm{Cov}(A_s, A_{s+t})$ が $s$ によらないことをいう。このとき

$$\Gamma(t) = \mathrm{Cov}(A_s, A_{s+|t|}) \quad (t \in \mathbb{Z}), \qquad \rho(t) = \frac{\Gamma(t)}{\Gamma(0)} \quad (\Gamma(0) > 0)$$

を自己共分散、自己相関関数という。$\sum_{t \in \mathbb{Z}}|\Gamma(t)| < \infty$ のとき

$$\sigma_\infty^2 = \sum_{t \in \mathbb{Z}}\Gamma(t), \qquad \tau_{\mathrm{int}} = \frac12 + \sum_{t=1}^{\infty}\rho(t) = \frac{\sigma_\infty^2}{2\Gamma(0)}$$

を長期分散、積分自己相関時間という。

$\tau_{\mathrm{int}}$ の定義は 定義 6.1［格子場のモンテカルロ法］ と同じです。独立な列では $\rho(t) = 0$（$t \ne 0$）なので $\tau_{\mathrm{int}} = 1/2$ です。標本平均 $\bar A_N = N^{-1}\sum_{s=1}^N A_s$ の分散は、$\sum_t|\Gamma(t)| < \infty$ のもとで

$$\mathrm{Var}(\bar A_N) = \frac1N\sum_{|t| < N}\Bigl(1 - \frac{|t|}{N}\Bigr)\Gamma(t), \qquad \lim_{N\to\infty}N\,\mathrm{Var}(\bar A_N) = \sigma_\infty^2 = 2\tau_{\mathrm{int}}\Gamma(0)$$

です（定理 6.2［格子場のモンテカルロ法］）。前半は $\mathrm{Var}(\sum_s A_s) = \sum_{s,u}\Gamma(s-u)$ で差 $t = s - u$ をもつ組が $N - |t|$ 個あることから出ます。後半は $(1 - |t|/N)\Gamma(t) \to \Gamma(t)$ と絶対総和可能性による項別の極限です。

定義 2.2（有効標本数）↗

$\sum_t|\Gamma(t)| < \infty$、$\sigma_\infty^2 > 0$ を満たす弱定常列の $N$ 個の標本について

$$N_{\mathrm{eff}} = \frac{N}{2\tau_{\mathrm{int}}} = \frac{N\,\Gamma(0)}{\sigma_\infty^2}$$

を有効標本数という。

上の式から $\mathrm{Var}(\bar A_N) \approx \Gamma(0)/N_{\mathrm{eff}}$ で、$N_{\mathrm{eff}}$ は「同じ分散を与える独立標本の個数」です。標準誤差は $\sqrt{\sigma_\infty^2/N}$ で、相関を無視した $\sqrt{\Gamma(0)/N}$ の $\sqrt{2\tau_{\mathrm{int}}}$ 倍です。したがって誤差を付ける問題は、$\sigma_\infty^2$ を推定する問題です。

2.2. 有限・可逆な連鎖

§3・§4・§8 では、次の設定の Markov 連鎖（定義 2.1［Markov 連鎖と Poisson 過程］）を考えます。

• $S$ は $n$ 点の有限集合、$P$ は $S$ 上の既約（どの状態からどの状態へも行ける）な遷移行列。
• $\pi$ は $S$ 上の確率分布で、すべての $x$ で $\pi(x) > 0$、$P$ について可逆、すなわち $\pi(x)P(x,y) = \pi(y)P(y,x)$（定義 2.4［Markov 連鎖と Poisson 過程］）。このとき $\pi$ は定常分布です（命題 2.5［Markov 連鎖と Poisson 過程］）。

$S$ 上の実関数に内積 $\langle f, g\rangle_\pi = \sum_x \pi(x)f(x)g(x)$ とノルム $\|f\|_\pi = \sqrt{\langle f, f\rangle_\pi}$ を入れ、$\pi(f) = \sum_x\pi(x)f(x)$、$\mathrm{Var}_\pi(f) = \|f - \pi(f)\|_\pi^2$ と書きます。$(Pf)(x) = \sum_y P(x,y)f(y)$ です。

補題 2.3（可逆な遷移行列の固有分解）✓

上の設定で、実数 $1 = \lambda_1 > \lambda_2 \ge \cdots \ge \lambda_n \ge -1$ と、$\langle\cdot,\cdot\rangle_\pi$ について正規直交な実関数 $\varphi_1, \ldots, \varphi_n$ があって、$P\varphi_j = \lambda_j\varphi_j$、$\varphi_1 \equiv 1$ となる。

証明（補題 2.3）

$D = \mathrm{diag}(\pi(x))$ とし、$K = D^{1/2}PD^{-1/2}$、すなわち $K(x,y) = \pi(x)^{1/2}P(x,y)\pi(y)^{-1/2}$ とおきます。可逆性から

$$K(x,y) = \frac{\pi(x)P(x,y)}{\sqrt{\pi(x)\pi(y)}} = \frac{\pi(y)P(y,x)}{\sqrt{\pi(x)\pi(y)}} = K(y,x)$$

で、$K$ は実対称行列です。スペクトル定理（定理 4.2［スペクトル定理］）から、実固有値 $\lambda_1 \ge \cdots \ge \lambda_n$ と正規直交な固有ベクトル $u_1, \ldots, u_n$ があります。$K$ が実行列で固有値が実数なので、$u_j$ は実ベクトルに取れます（$u_j$ の実部と虚部はどちらも同じ固有値の固有ベクトルで、各固有空間の実ベクトルからなる基底に Gram–Schmidt の直交化をすれば実の正規直交基底になります）。$\varphi_j = D^{-1/2}u_j$ とおくと $P\varphi_j = D^{-1/2}Ku_j = \lambda_j\varphi_j$、$\langle\varphi_i, \varphi_j\rangle_\pi = \sum_x u_i(x)u_j(x) = \delta_{ij}$ です。

固有値の範囲を調べます。$P\varphi = \lambda\varphi$、$\varphi \ne 0$ とし、$|\varphi|$ が最大になる点を $x_0$ とします。$|\lambda||\varphi(x_0)| = |\sum_y P(x_0,y)\varphi(y)| \le \sum_y P(x_0, y)|\varphi(x_0)| = |\varphi(x_0)|$ から $|\lambda| \le 1$ です。

定数関数 $1$ は $P1 = 1$ を満たし $\|1\|_\pi = 1$ です。固有値 $1$ の固有関数が定数に限ることを示します。$P\varphi = \varphi$ とし、$\varphi$ が最大値 $m$ を取る点を $x_0$ とします。$m = \sum_y P(x_0,y)\varphi(y)$ は $m$ 以下の数の重みつき平均なので、$P(x_0,y) > 0$ であるすべての $y$ で $\varphi(y) = m$ です。これを繰り返すと、$x_0$ から正の確率でたどり着けるすべての点で $\varphi = m$ となり、既約性から $\varphi \equiv m$ です。したがって固有値 $1$ の重複度は $1$ で、$\varphi_1 \equiv 1$、$\lambda_2 < 1$ と取れます。

∎

定義 2.4（絶対スペクトルギャップ）↗

補題 2.3 の固有値について $\gamma_* = 1 - \max(|\lambda_2|, |\lambda_n|)$ を絶対スペクトルギャップという（$n \ge 2$）。

通常のスペクトルギャップ $\gamma = 1 - \lambda_2$（命題 6.2［Markov 連鎖と Poisson 過程］）との違いは、$-1$ に近い固有値も数えることです。確率 $1/2$ でその場に留まる「怠惰な」連鎖 $(I + P)/2$ では固有値がすべて $0$ 以上になり、$\gamma_* = \gamma$ です。

3. 可逆な連鎖の自己共分散

3.1. スペクトル表示と対の和の正値性

ここでの狙いは、自己共分散 $\Gamma(t)$ を連鎖の固有値で書き表し、推定値のどこまでを信用してよいかの判定規則を得ることです。推定値は大きなラグで雑音だけになるので、足し上げをどこかで止める必要があります。隣り合う 2 項 $\Gamma(2k) + \Gamma(2k+1)$ をまとめると、奇数ラグの負の成分が偶数ラグの正の成分と打ち消し合い、必ず非負になります。この性質が「和が負になったら止める」という規則の根拠になります。

連鎖を定常分布から始め（$X_0 \sim \pi$）、関数 $f : S \to \mathbb{R}$ について $A_t = f(X_t)$ とおきます。

定理 3.1（可逆な連鎖の自己共分散（Geyer））✓

§2.2 の設定で $X_0 \sim \pi$、$A_t = f(X_t)$ とし、$\tilde f = f - \pi(f)$、$c_j = \langle\tilde f, \varphi_j\rangle_\pi^2$（$j = 2, \ldots, n$）とおく。

1. $(A_t)$ は弱定常で、$t \ge 0$ について $$\Gamma(t) = \langle\tilde f, P^t\tilde f\rangle_\pi = \sum_{j=2}^{n}c_j\lambda_j^{t}, \qquad c_j \ge 0, \qquad \sum_{j=2}^n c_j = \Gamma(0) = \mathrm{Var}_\pi(f).$$
2. $G(k) = \Gamma(2k) + \Gamma(2k+1)$（$k = 0, 1, 2, \ldots$）は $$G(k) \ge 0, \qquad G(k) - G(k+1) \ge 0, \qquad G(k) - 2G(k+1) + G(k+2) \ge 0$$ を満たす。$c_j > 0$ となる $j$ のうちに $\lambda_j \notin \{0, -1\}$ のものがあれば、すべての $k$ で $G(k) > 0$ である。

証明（定理 3.1）

1. $\pi$ は定常分布なので、すべての $s$ で $X_s \sim \pi$ です。したがって $\mathbb{E}[A_s] = \pi(f)$ です。$n$ ステップの遷移確率が $P^t$ で与えられること（命題 2.3［Markov 連鎖と Poisson 過程］）から $\mathbb{P}(X_s = x, X_{s+t} = y) = \pi(x)P^t(x,y)$ で、

$$\mathrm{Cov}(A_s, A_{s+t}) = \sum_{x,y}\pi(x)P^t(x,y)\tilde f(x)\tilde f(y) = \langle\tilde f, P^t\tilde f\rangle_\pi$$

は $s$ によりません。補題 2.3 の正規直交基底で $\tilde f = \sum_{j=1}^n\langle\tilde f, \varphi_j\rangle_\pi\varphi_j$ と展開します。$\langle\tilde f, \varphi_1\rangle_\pi = \pi(\tilde f) = 0$ なので $j = 1$ の項はありません。$P^t\varphi_j = \lambda_j^t\varphi_j$ と正規直交性から

$$\langle\tilde f, P^t\tilde f\rangle_\pi = \sum_{j=2}^n\langle\tilde f, \varphi_j\rangle_\pi^2\lambda_j^t = \sum_{j=2}^n c_j\lambda_j^t$$

です。$t = 0$ とすれば $\sum_j c_j = \|\tilde f\|_\pi^2 = \mathrm{Var}_\pi(f)$ です。

2. 1 から

$$G(k) = \sum_{j=2}^n c_j\lambda_j^{2k}(1 + \lambda_j)$$

です。$c_j \ge 0$、$\lambda_j^{2k} \ge 0$、$1 + \lambda_j \ge 0$（$\lambda_j \ge -1$）なので各項は非負です。同様に

$$G(k) - G(k+1) = \sum_{j=2}^n c_j\lambda_j^{2k}(1 + \lambda_j)(1 - \lambda_j^2), \qquad G(k) - 2G(k+1) + G(k+2) = \sum_{j=2}^n c_j\lambda_j^{2k}(1 + \lambda_j)(1 - \lambda_j^2)^2$$

で、$|\lambda_j| \le 1$ から $1 - \lambda_j^2 \ge 0$ なので、どちらも非負です。$c_j > 0$ かつ $\lambda_j \notin \{0, -1\}$ の $j$ があれば、$G(k)$ の和のその項が $c_j\lambda_j^{2k}(1 + \lambda_j) > 0$ なので $G(k) > 0$ です。

∎

$\Gamma(t)$ そのものは正とは限りません。負の固有値の成分があると、奇数のラグで負になります。

例 3.2（奇数ラグの負の相関）

$S = \{1, 2\}$、$P = \begin{pmatrix}1-a & a\\ a & 1-a\end{pmatrix}$（$0 < a < 1$）は $\pi = (1/2, 1/2)$ について可逆です。固有値は $1$ と $\lambda_2 = 1 - 2a$、$\varphi_2 = (1, -1)$ です。$f = \varphi_2$ とすると $\Gamma(t) = (1 - 2a)^t$ で、$a > 1/2$ なら奇数の $t$ で負です。一方 $G(k) = (1-2a)^{2k}(2 - 2a) > 0$ です。このとき

$$\tau_{\mathrm{int}} = \frac12 + \sum_{t\ge1}(1-2a)^t = \frac12 + \frac{1 - 2a}{2a} = \frac{1 - a}{2a}$$

で、$a > 1/2$ では $\tau_{\mathrm{int}} < 1/2$、つまり独立な標本より分散が小さくなります。$a = 1$ は交互に入れ替わる連鎖（$\lambda_2 = -1$）で、$\tau_{\mathrm{int}}$ を定義する級数は収束しませんが、$G(k) = 0$ です。

3.2. 初期正値列による打ち切り

$\tau_{\mathrm{int}}$ をデータから求めるには、自己共分散の推定値を足し上げます。ところが $\Gamma(t)$ の推定値の雑音はラグ $t$ が大きくなっても小さくならないので（真の $\Gamma(t)$ は指数的に減衰します）、全部足すと雑音の和になります。どこかで打ち切る必要があり、定理 3.1 がその規則を与えます。

定義 3.3（自己共分散の推定量・窓つき推定量）↗

データ $A_1, \ldots, A_N$ と $0 \le t < N$ について

$$\hat\Gamma(t) = \frac{1}{N - t}\sum_{s=1}^{N-t}(A_s - \bar A_N)(A_{s+t} - \bar A_N)$$

とおく。$0 \le W < N$ について

$$\hat\sigma_W^2 = \hat\Gamma(0) + 2\sum_{t=1}^{W}\hat\Gamma(t), \qquad \hat\tau_W = \frac{\hat\sigma_W^2}{2\hat\Gamma(0)} = \frac12 + \sum_{t=1}^W\frac{\hat\Gamma(t)}{\hat\Gamma(0)}$$

を、窓 $W$ の長期分散・積分自己相関時間の推定量という。

定義 3.4（初期正値列）↗

$\hat G(k) = \hat\Gamma(2k) + \hat\Gamma(2k+1)$ とし、$\hat G(k) \le 0$ となる最小の $k$ を $K$ とする。

$$\hat\tau_{\mathrm{IPS}} = -\frac12 + \frac{1}{\hat\Gamma(0)}\sum_{k=0}^{K-1}\hat G(k)$$

を初期正値列（initial positive sequence）による推定量という。窓 $W = 2K - 1$ の $\hat\tau_W$ に等しい。

真の値では $\sum_{k\ge0}G(k) = \sum_{t\ge0}\Gamma(t) = \Gamma(0)(\tau_{\mathrm{int}} + 1/2)$ なので、$-1/2$ を足せば $\tau_{\mathrm{int}}$ です。定理 3.1 から真の $G(k)$ は負になりません。したがって推定値 $\hat G(k)$ が $0$ 以下になったら、そのラグでは信号が雑音に埋もれています。そこで打ち切るのが初期正値列です。単調性と凸性も使い、$\hat G(k)$ を単調非増加・凸になるように直してから足す方法（initial monotone / convex sequence）もあります。どれも、可逆な連鎖であることを使っています。可逆でない連鎖（系統走査の熱浴法など）では $G(k) \ge 0$ の保証がないので、次節の一般の窓かバッチ平均を使います。

4. 窓の偏りと分散

命題 4.1（窓つき推定量の偏り）✓

弱定常列 $(A_t)$ が $S_1 = \sum_{t\in\mathbb{Z}}|\Gamma(t)| < \infty$ を満たすとし、$0 \le W < N$ とする。

1. $0 \le t < N$ について $|\mathbb{E}[\hat\Gamma(t)] - \Gamma(t)| \le 3S_1/N$。
2. $\mathbb{E}[\hat\sigma_W^2] - \sigma_\infty^2 = -2\sum_{t > W}\Gamma(t) + R$、$|R| \le 3(2W+1)S_1/N$。したがって $\Gamma(0) > 0$ のとき $$\frac{\mathbb{E}[\hat\sigma_W^2]}{2\Gamma(0)} - \tau_{\mathrm{int}} = -\sum_{t > W}\rho(t) + R', \qquad |R'| \le \frac{3(2W+1)S_1}{2\Gamma(0)N}.$$
3. $(A_t)$ が 定理 3.1 の列で絶対スペクトルギャップが $\gamma_* > 0$ のとき、$|\rho(t)| \le (1 - \gamma_*)^t$ で、 $$\Bigl|\sum_{t > W}\rho(t)\Bigr| \le \frac{(1 - \gamma_*)^{W+1}}{\gamma_*} \le \frac{(1-\gamma_*)^W}{\gamma_*}.$$

証明（命題 4.1）

1. $1 \le s \le N - t$ とします。共分散の双線形性（命題 5.3［確率変数と期待値］）から

$$\mathbb{E}\bigl[(A_s - \bar A_N)(A_{s+t} - \bar A_N)\bigr] = \Gamma(t) - \mathrm{Cov}(A_s, \bar A_N) - \mathrm{Cov}(A_{s+t}, \bar A_N) + \mathrm{Var}(\bar A_N)$$

です（各因子の期待値は $0$ なので、左辺は共分散です）。$\mathrm{Cov}(A_s, \bar A_N) = N^{-1}\sum_{u=1}^N\Gamma(u - s)$ の絶対値は $S_1/N$ 以下です。$A_{s+t}$ についても同じです。$\mathrm{Var}(\bar A_N) = N^{-2}\sum_{u=1}^N\sum_{v=1}^N\Gamma(u - v)$ は、$u$ を固定した $v$ についての和が絶対値で $S_1$ 以下なので $S_1/N$ 以下です。したがって各 $s$ で期待値と $\Gamma(t)$ の差は $3S_1/N$ 以下で、$s$ についての平均 $\mathbb{E}[\hat\Gamma(t)]$ でも同じです。

2. $\sigma_\infty^2 = \Gamma(0) + 2\sum_{t\ge1}\Gamma(t)$ なので

$$\mathbb{E}[\hat\sigma_W^2] - \sigma_\infty^2 = \bigl(\mathbb{E}\hat\Gamma(0) - \Gamma(0)\bigr) + 2\sum_{t=1}^W\bigl(\mathbb{E}\hat\Gamma(t) - \Gamma(t)\bigr) - 2\sum_{t > W}\Gamma(t)$$

です。はじめの 2 つの和を $R$ とおくと、1 から $|R| \le (1 + 2W)\cdot3S_1/N$ です。両辺を $2\Gamma(0)$ で割れば後半です。

3. 定理 3.1 から $\rho(t) = \sum_{j\ge2}c_j\lambda_j^t/\sum_{j\ge2}c_j$ で、$|\lambda_j| \le 1 - \gamma_*$ なので $|\rho(t)| \le (1-\gamma_*)^t$ です。等比級数の和から $\sum_{t>W}(1-\gamma_*)^t = (1-\gamma_*)^{W+1}/\gamma_*$ です。

∎

偏りの主要部は捨てた裾 $-\sum_{t>W}\rho(t)$ です。正の相関だけなら、窓つき推定量は $\tau_{\mathrm{int}}$ を小さく見積もり、誤差棒は小さすぎる側にずれます。3 から、裾を $\varepsilon$ 以下にするには $W \ge \ln(1/(\varepsilon\gamma_*))/\gamma_*$ 程度、つまり緩和時間 $1/\gamma_*$ の数倍の窓が要ります。

注意 4.2（窓つき推定量の分散と、偏りとの兼ね合い）

命題 4.1 が評価しているのは分子 $\hat\sigma_W^2$ の期待値です。$\hat\tau_W$ は $\hat\Gamma(0)$ との比なので、比を取ることによる $O(1/N)$ の偏り（命題 6.6 と同じ仕組み）が加わります。

分散については、$\tau_{\mathrm{int}} \ll W \ll N$ の範囲で

$$\mathrm{Var}(\hat\tau_W) \approx \frac{2(2W+1)}{N}\,\tau_{\mathrm{int}}^2$$

という Madras–Sokal の公式が使われます（N. Madras, A. D. Sokal, J. Stat. Phys. 50 (1988) の付録。導出は 4 次のキュムラントを落とす近似によるもので、この記事では定理にしません）。相対誤差は $\sqrt{2(2W+1)/N}$ で、窓に比例して分散が増えます。偏りは $W$ とともに指数的に減り、分散は $W$ に比例して増えるので、$W$ は $\tau_{\mathrm{int}}$ の数倍が妥当です。Sokal の規則は「$W \ge c\,\hat\tau_W$ を満たす最小の $W$」（$c$ は 4〜10）を選びます。U. Wolff の Γ 法は、偏りと統計誤差の和の見積もりを最小にする $W$ を自動で選び、平均の非線形な関数にも適用できます（Comput. Phys. Commun. 156 (2004)）。§9 の実験で、$\rho = 0.9$ の AR(1) 過程についてこの公式と実測のばらつきを比べます。

窓 W の選び方。W までの自己相関を足し、それより先の裾を捨てる。捨てた裾の面積が偏り、足した項の個数に比例して分散が増える。

5. バッチ平均

自己共分散を 1 つずつ推定せずに $\sigma_\infty^2$ を求める方法がバッチ平均です。データを長さ $b$ の区間に切り、区間ごとの平均を新しいデータと見ます。$b$ が相関の長さより十分大きければ、区間の平均どうしはほぼ無相関で、独立な標本の式が使えます。

定義 5.1（バッチ平均）↗

$N = kb$（$k \ge 2$）とし、$Y_j = b^{-1}\sum_{s=(j-1)b+1}^{jb}A_s$（$j = 1, \ldots, k$）をバッチ平均という。$\bar Y = k^{-1}\sum_j Y_j = \bar A_N$ として

$$\hat\sigma_{\mathrm{BM}}^2 = \frac{b}{k-1}\sum_{j=1}^k(Y_j - \bar Y)^2$$

を、長期分散 $\sigma_\infty^2$のバッチ平均推定量という。

これはビニング推定量 $\hat\sigma_B^2$（定義 6.4［格子場のモンテカルロ法］。$\mathrm{Var}(\bar A_N)$ の推定量 $\frac{1}{k(k-1)}\sum_j(Y_j - \bar Y)^2$）の $N$ 倍です。期待値が $b \to \infty$ で $\sigma_\infty^2$ に収束することは 命題 6.5［格子場のモンテカルロ法］ にあります。ここでは収束の速さを示します。

次の定理は 2 つのことを言います。1 つは、バッチ長 $b$ を大きくすれば推定量の偏りが $1/b$ の速さで消えること。もう 1 つは、バッチの個数が少ないままでは、$b$ をいくら大きくしても推定値のばらつきが残ることです。後者は、独立な正規標本を 2 つのバッチに分けると、推定値が「標準正規分布に従う数の 2 乗」の定数倍になって、データを増やしても定数に近づかないことから分かります。

定理 5.2（バッチ平均の偏りと、バッチ数を固定したときの不一致）✓

弱定常列 $(A_t)$ が $C_1 = \sum_{t=1}^\infty t\,|\Gamma(t)| < \infty$ を満たすとする。

1. すべての $b \ge 1$ で $\bigl|b\,\mathrm{Var}(Y_1) - \sigma_\infty^2\bigr| \le 2C_1/b$。
2. $N = kb$、$k \ge 2$ のとき $\bigl|\mathbb{E}[\hat\sigma_{\mathrm{BM}}^2] - \sigma_\infty^2\bigr| \le 5C_1/b$。
3. $A_t$ が独立で同じ正規分布 $N(\mu, \sigma^2)$（$\sigma > 0$）に従い、$k = 2$ のとき、すべての $b$ で $\hat\sigma_{\mathrm{BM}}^2$ は $\sigma^2Z^2$（$Z \sim N(0,1)$）と同じ分布に従う。特に $b \to \infty$ としても $\hat\sigma_{\mathrm{BM}}^2$ は $\sigma_\infty^2 = \sigma^2$ に確率収束しない。

証明（定理 5.2）

1. $C_1 < \infty$ なら $\sum_t|\Gamma(t)| \le |\Gamma(0)| + 2C_1 < \infty$ です。§2.1 の分散の式（定理 6.2［格子場のモンテカルロ法］）を $N = b$ で使うと $b\,\mathrm{Var}(Y_1) = \sum_{|t|<b}(1 - |t|/b)\Gamma(t)$ なので

$$\sigma_\infty^2 - b\,\mathrm{Var}(Y_1) = \sum_{|t|\ge b}\Gamma(t) + \frac1b\sum_{|t|<b}|t|\,\Gamma(t)$$

です。第 1 項は、$|t| \ge b$ で $1 \le |t|/b$ なので絶対値が $b^{-1}\sum_{|t|\ge b}|t||\Gamma(t)|$ 以下です。第 2 項と合わせて

$$\bigl|\sigma_\infty^2 - b\,\mathrm{Var}(Y_1)\bigr| \le \frac1b\sum_{t\in\mathbb{Z}}|t|\,|\Gamma(t)| = \frac{2C_1}{b}.$$

2. $\mathbb{E}[Y_j] = \mathbb{E}[\bar Y] = \mu$ です。$\sum_j(Y_j - \bar Y)^2 = \sum_j(Y_j - \mu)^2 - k(\bar Y - \mu)^2$ の期待値を取ると、弱定常性から $\mathrm{Var}(Y_j) = \mathrm{Var}(Y_1)$ なので

$$\mathbb{E}[\hat\sigma_{\mathrm{BM}}^2] = \frac{b}{k-1}\bigl(k\,\mathrm{Var}(Y_1) - k\,\mathrm{Var}(\bar A_N)\bigr) = \frac{k\,v_b - v_N}{k-1}, \qquad v_m = m\,\mathrm{Var}(\bar A_m)$$

です（$bk\,\mathrm{Var}(\bar A_N) = v_N$、$\bar A_b$ は $Y_1$ のこと）。1 を $b$ と $N$ に使うと $|v_b - \sigma_\infty^2| \le 2C_1/b$、$|v_N - \sigma_\infty^2| \le 2C_1/N = 2C_1/(kb)$ です。したがって

$$\bigl|\mathbb{E}[\hat\sigma_{\mathrm{BM}}^2] - \sigma_\infty^2\bigr| = \frac{|k(v_b - \sigma_\infty^2) - (v_N - \sigma_\infty^2)|}{k-1} \le \frac{2C_1}{b}\cdot\frac{k + 1/k}{k-1}$$

です。$(k + 1/k)/(k-1) = 1 + (1 + 1/k)/(k - 1)$ は $k \ge 2$ で単調に減り、$k = 2$ での値は $5/2$ なので、右辺は $5C_1/b$ 以下です。

3. $k = 2$ では $\bar Y = (Y_1 + Y_2)/2$ から $(Y_1 - \bar Y)^2 + (Y_2 - \bar Y)^2 = (Y_1 - Y_2)^2/2$ で、$\hat\sigma_{\mathrm{BM}}^2 = b(Y_1 - Y_2)^2/2$ です。$(A_1, \ldots, A_{2b})$ は独立な正規確率変数を並べたものなので、同時密度は密度の積（補題 2.2［分布からの標本抽出］）で、これは $N(\mu\boldsymbol{1}, \sigma^2I)$ の密度です（定理 5.2［特性関数と Gauss ベクトル］ の 3）。同じ定理の 4 から、線形結合 $Y_1 - Y_2 = b^{-1}(\sum_{s\le b}A_s - \sum_{s>b}A_s)$ は平均 $0$、分散 $2b\sigma^2/b^2 = 2\sigma^2/b$ の正規分布に従います。$Z = (Y_1 - Y_2)\sqrt{b/(2\sigma^2)}$ とおくと $Z \sim N(0,1)$ で $\hat\sigma_{\mathrm{BM}}^2 = \sigma^2Z^2$ です。$\mathbb{P}(|\hat\sigma_{\mathrm{BM}}^2 - \sigma^2| > \sigma^2/2) = \mathbb{P}(|Z^2 - 1| > 1/2)$ は $b$ によらない正の数なので、確率収束しません。

∎

1 と 2 は「$b$ を大きくすれば偏りが消える」こと、3 は「$b$ だけ大きくしても推定量のばらつきは消えない」ことを、正規・独立でバッチが 2 個の場合に言っています。一般のバッチ数 $k$ でも、正規・独立な場合には $(k-1)\hat\sigma_{\mathrm{BM}}^2/\sigma^2$ が自由度 $k-1$ の $\chi^2$ 分布に従い（正規標本の標本分散の分布についての標準的な事実で、この記事では証明しません）、相対的なばらつきは $\sqrt{2/(k-1)}$ です。$k = 20$ でも誤差棒の誤差が 3 割あります。一致性には $b \to \infty$ と $k = N/b \to \infty$ の両方が要ります。

注意 5.3（バッチ長の選び方と一致性の条件）

偏りが $O(1/b)$、分散が $O(b/N)$ なので、平均二乗誤差を最小にするバッチ長は $b \propto N^{1/3}$ です。一般の Markov 連鎖でバッチ平均推定量が強一致性をもつための条件（幾何エルゴード性と積率の条件）と、スペクトル分散推定量との比較は J. M. Flegal, G. L. Jones, Ann. Statist. 38 (2010) にあります。証明は連鎖の強近似に基づくので、この記事では定理にしません。実務では、$b$ を 2 倍ずつ増やして推定値が頭打ちになる（定理 5.2 の偏りが消える）ことと、そのとき $k$ が数十以上残っていることの両方を確かめます。

例 5.4（AR(1) 過程では全部が閉じた形で分かる）

$|\varrho| < 1$ とし、$A_{t+1} = \varrho A_t + \sqrt{1 - \varrho^2}\,Z_t$（$Z_t$ は独立な標準正規、$A_1 \sim N(0,1)$ は $Z_t$ と独立）とします。$A_t \sim N(0,1)$ がすべての $t$ で成り立ち、$\mathrm{Cov}(A_s, A_{s+t}) = \varrho^t$ です（$A_{s+t} = \varrho^tA_s + (Z_s, \ldots, Z_{s+t-1}$ の線形結合$)$ で、後者は $A_s$ と独立だからです）。したがって

$$\rho(t) = \varrho^{|t|}, \qquad \tau_{\mathrm{int}} = \frac12 + \frac{\varrho}{1 - \varrho} = \frac{1 + \varrho}{2(1 - \varrho)}, \qquad C_1 = \sum_{t\ge1}t\varrho^t = \frac{\varrho}{(1-\varrho)^2}\ (\varrho > 0).$$

$\varrho = 0.9$ なら $\tau_{\mathrm{int}} = 9.5$ で、有効標本数は $N/19$ です。窓の偏りは $-\sum_{t>W}\varrho^t = -\varrho^{W+1}/(1-\varrho)$ で、$W = 20$ のとき $-1.094$ です。バッチ平均は、等比級数の和 $\sum_{|t|<b}(1 - |t|/b)\varrho^{|t|}$ を計算して（演習 10.1）

$$\frac{b\,\mathrm{Var}(Y_1)}{2\Gamma(0)} = \tau_{\mathrm{int}} - \frac{\varrho(1 - \varrho^b)}{b(1-\varrho)^2}$$

です。$b = 16$ で $4.917$、$b = 256$ で $9.148$ となり、偏りが $1/b$ で減ることが見えます。

6. ジャックナイフ

6.1. 分散推定量の一致性

報告したい量が平均の関数 $g(\mu)$ で、推定量が $g(\bar X_N)$ のとき、その誤差をデータの作り直しで求めるのがジャックナイフです。この節と次節では独立同分布な標本を考えます。相関のあるデータへの適用は §6.3 で述べます。

定義 6.1（ジャックナイフ分散推定量）↗

データ $X_1, \ldots, X_N$（$N \ge 2$）と統計量 $\hat\theta = T(X_1, \ldots, X_N)$ について、$i$ 番目を除いた $N - 1$ 個のデータで同じ統計量を計算した値を $\hat\theta_{(i)}$、その平均を $\bar\theta_{(\cdot)} = N^{-1}\sum_i\hat\theta_{(i)}$ とする。

$$v_{\mathrm J} = \frac{N-1}{N}\sum_{i=1}^N\bigl(\hat\theta_{(i)} - \bar\theta_{(\cdot)}\bigr)^2, \qquad \hat\theta_{\mathrm J} = N\hat\theta - (N-1)\bar\theta_{(\cdot)}$$

を、$\hat\theta$ の分散のジャックナイフ推定量、偏りを補正したジャックナイフ推定値という。

$\hat\theta = g(\bar X_N)$ のときは $\hat\theta_{(i)} = g(\bar X_{(i)})$、$\bar X_{(i)} = (N-1)^{-1}\sum_{j\ne i}X_j$ で、定義 6.6［格子場のモンテカルロ法］ と同じ定義です。係数 $(N-1)/N$ は、$\hat\theta_{(i)}$ どうしが $N - 2$ 個のデータを共有していて互いに近い（差が $1/N$ の大きさ）ことを補う因子です。$g(x) = x$ では $v_{\mathrm J}$ が $s_N^2/N$ に一致します（演習 10.2）。

補題 6.2（最大値は N より遅くしか増えない）✓

$X_1, X_2, \ldots$ が独立同分布で $\mathbb{E}|X_1| < \infty$ なら、確率 $1$ で $\max_{1\le i\le N}|X_i|/N \to 0$ である。

証明（補題 6.2）

$S_N = \sum_{i\le N}X_i$ とします。大数の強法則（定理 4.5［大数の法則と中心極限定理］）から確率 $1$ で $S_N/N \to \mu$ で、

$$\frac{X_N}{N} = \frac{S_N}{N} - \frac{N-1}{N}\cdot\frac{S_{N-1}}{N-1} \to \mu - \mu = 0$$

です。この事象の上で $\varepsilon > 0$ を取ると、ある $n_0$ があって $i \ge n_0$ なら $|X_i| \le \varepsilon i$ です。$N \ge n_0$ のとき

$$\frac1N\max_{i\le N}|X_i| \le \max\Bigl(\frac1N\max_{i<n_0}|X_i|,\ \max_{n_0\le i\le N}\frac{\varepsilon i}{N}\Bigr) \le \max\Bigl(\frac1N\max_{i<n_0}|X_i|,\ \varepsilon\Bigr)$$

で、$N \to \infty$ で第 1 の項は $0$ に向かうので、上極限は $\varepsilon$ 以下です。$\varepsilon$ は任意なので極限は $0$ です。

∎

定理 6.3（ジャックナイフ分散推定量の一致性）✓

$X_1, X_2, \ldots$ を独立同分布な実確率変数、$\mu = \mathbb{E}[X_1]$、$\sigma^2 = \mathrm{Var}(X_1) < \infty$ とする。$g$ を $\mu$ を含む開区間 $I$ 上の $C^1$ 級関数とし、$\hat\theta = g(\bar X_N)$ のジャックナイフ分散推定量を $v_{\mathrm J}$ とする。このとき確率 $1$ で、十分大きいすべての $N$ について $\bar X_N$ と $\bar X_{(1)}, \ldots, \bar X_{(N)}$ は $I$ に入り（したがって $v_{\mathrm J}$ が定まり）、

$$\lim_{N\to\infty}N\,v_{\mathrm J} = g'(\mu)^2\sigma^2.$$

証明の方針は次のとおりです。1 個除いた平均 $\bar X_{(i)}$ は全体の平均からほんのわずかしかずれないので、$g$ を平均値の定理で 1 次式に置き換えられます。すると $v_{\mathrm J}$ は「$g'(\bar X_N)^2$ と標本分散の積を $N$ で割ったもの」にほぼ等しくなり、残りの誤差は 0 に向かうことを示します。

証明（定理 6.3）

$\Delta_i = \bar X_{(i)} - \bar X_N$ とおきます。$\bar X_{(i)} = (N\bar X_N - X_i)/(N-1)$ なので

$$\Delta_i = \frac{\bar X_N - X_i}{N-1}, \qquad \sum_{i=1}^N\Delta_i = 0, \qquad (N-1)\sum_{i=1}^N\Delta_i^2 = \frac{1}{N-1}\sum_{i=1}^N(X_i - \bar X_N)^2 = s_N^2$$

です。以下、次の 3 つが成り立つ確率 $1$ の事象の上で考えます。(a) $\bar X_N \to \mu$（定理 4.5［大数の法則と中心極限定理］）、(b) $s_N^2 \to \sigma^2$（定理 2.2［Monte Carlo 積分と分散低減］ の 3）、(c) $\max_{i\le N}|X_i|/N \to 0$（補題 6.2）。

(a) と (c) から $\max_i|\Delta_i| \le (|\bar X_N| + \max_i|X_i|)/(N-1) \to 0$ です。$\eta_N = |\bar X_N - \mu| + \max_i|\Delta_i|$ とおくと $\eta_N \to 0$ で、$\bar X_N$ と各 $\bar X_{(i)}$ は $[\mu - \eta_N, \mu + \eta_N]$ に入ります。$I$ は $\mu$ を含む開区間なので、十分大きい $N$ でこの区間は $I$ に含まれます。以下そのような $N$ を考えます。

平均値の定理（定理 3.3［平均値の定理とテイラーの定理］）から、$\bar X_N$ と $\bar X_{(i)}$ の間の点 $\xi_i$ があって $g(\bar X_{(i)}) - g(\bar X_N) = g'(\xi_i)\Delta_i$ です（$\Delta_i = 0$ のときは $\xi_i = \bar X_N$ とします）。

$$a_i = g'(\bar X_N)\Delta_i, \qquad r_i = \bigl(g'(\xi_i) - g'(\bar X_N)\bigr)\Delta_i$$

とおくと $\hat\theta_{(i)} = g(\bar X_N) + a_i + r_i$ で、$\sum_ia_i = 0$ から $\bar\theta_{(\cdot)} = g(\bar X_N) + \bar r$（$\bar r = N^{-1}\sum_ir_i$）、したがって $\hat\theta_{(i)} - \bar\theta_{(\cdot)} = a_i + (r_i - \bar r)$ です。$\mathbb{R}^N$ のベクトル $a = (a_i)$、$r' = (r_i - \bar r)$ とユークリッドノルムについて

$$\sqrt{N\,v_{\mathrm J}} = \sqrt{N-1}\,\|a + r'\|, \qquad \bigl|\sqrt{N\,v_{\mathrm J}} - \sqrt{N-1}\,\|a\|\bigr| \le \sqrt{N-1}\,\|r'\|$$

です（三角不等式）。第 1 の量は

$$(N-1)\|a\|^2 = g'(\bar X_N)^2(N-1)\sum_i\Delta_i^2 = g'(\bar X_N)^2s_N^2 \to g'(\mu)^2\sigma^2$$

です（(a)(b) と $g'$ の連続性）。第 2 の量は、平均を引くと 2 乗和が減ること（$\sum_i(r_i - \bar r)^2 = \sum_ir_i^2 - N\bar r^2$）から

$$(N-1)\|r'\|^2 \le (N-1)\sum_ir_i^2 \le \varepsilon_N^2(N-1)\sum_i\Delta_i^2 = \varepsilon_N^2s_N^2, \qquad \varepsilon_N = \max_i|g'(\xi_i) - g'(\bar X_N)|$$

です。$\xi_i$ と $\bar X_N$ は $[\mu - \eta_N, \mu + \eta_N]$ にあるので $\varepsilon_N \le 2\sup_{|x - \mu|\le\eta_N}|g'(x) - g'(\mu)|$ で、$g'$ が $\mu$ で連続、$\eta_N \to 0$ なので $\varepsilon_N \to 0$ です。よって $(N-1)\|r'\|^2 \to 0$ で、$\sqrt{N\,v_{\mathrm J}} \to |g'(\mu)|\sigma$ です。

∎

中心極限定理（定理 5.6［大数の法則と中心極限定理］）と delta 法（定理 7.1［点推定と漸近理論］）から $\sqrt N(g(\bar X_N) - g(\mu))$ は $N(0, g'(\mu)^2\sigma^2)$ に分布収束します。定理 6.3 は、ジャックナイフが $g'$ を計算せずにこの分散を正しく推定することを言っています。$g'(\mu) \ne 0$ なら、Slutsky の定理（命題 6.1［大数の法則と中心極限定理］）から $(g(\bar X_N) - g(\mu))/\sqrt{v_{\mathrm J}}$ が標準正規分布に分布収束し、$g(\bar X_N) \pm z\sqrt{v_{\mathrm J}}$ が漸近的な信頼区間になります。ベクトル値の平均の関数（比、分散、対数の差）でも、同じ証明が勾配と共分散行列について通ります。

証明で使ったのは $g$ の微分可能性です。滑らかでない統計量では結論が壊れます。

例 6.4（中央値ではジャックナイフ分散推定量は一致しない）

$X_1, \ldots, X_N$ を $(0,1)$ 上の一様分布からの独立な標本、$N = 2m$ とし、順序統計量を $X_{(1)} < \cdots < X_{(N)}$ と書きます（確率 $1$ で値はすべて異なります）。標本中央値を、偶数個なら中央の 2 個の平均、奇数個なら中央の値とします。$\hat\theta = (X_{(m)} + X_{(m+1)})/2$ です。

1 個を除くと $2m - 1$ 個が残り、中央値は小さいほうから $m$ 番目です。$X_{(i)}$（$i \le m$）を除くと残りの $m$ 番目は $X_{(m+1)}$、$X_{(i)}$（$i \ge m+1$）を除くと $X_{(m)}$ です。したがって $\hat\theta_{(i)}$ は $X_{(m+1)}$ と $X_{(m)}$ を $m$ 回ずつ取り、$\bar\theta_{(\cdot)} = \hat\theta$ で

$$v_{\mathrm J} = \frac{N-1}{N}\cdot N\cdot\Bigl(\frac{X_{(m+1)} - X_{(m)}}{2}\Bigr)^2, \qquad N\,v_{\mathrm J} = \frac{N-1}{N}\cdot\frac{(N\,D_N)^2}{4}, \quad D_N = X_{(m+1)} - X_{(m)}$$

です。中央の間隔 $D_N$ の分布を求めます。$(X_{(1)}, \ldots, X_{(N)})$ は集合 $\{0 < x_1 < \cdots < x_N < 1\}$（体積 $1/N!$）の上の一様分布に従います（$N!$ 通りの並べ替えが同じ確率だからです）。$0 < x < 1$ について、事象 $D_N > x$ に対応する集合を、$m+1$ 番目以降の座標から $x$ を引く平行移動で移すと $\{0 < y_1 < \cdots < y_N < 1 - x\}$ になり、その体積は $(1-x)^N/N!$ です。よって

$$\mathbb{P}(D_N > x) = (1 - x)^N, \qquad \mathbb{P}(N D_N > x) = \Bigl(1 - \frac xN\Bigr)^N \to e^{-x}$$

で、$ND_N$ は平均 $1$ の指数分布に分布収束します。したがって $N\,v_{\mathrm J}$ は $E^2/4$（$E$ は平均 $1$ の指数分布）に分布収束します。極限は定数ではなく、$N\,v_{\mathrm J}$ はどんな定数にも確率収束しません。極限の平均は $\mathbb{E}[E^2]/4 = 1/2$ で、標本中央値の漸近分散 $1/(4N)$（密度が中央値で $1$ の場合）の $N$ 倍 $1/4$ とも合いません。

6.2. 比の推定量の偏り

2 つの平均の比 $\bar X/\bar Y$ は、再重みづけ・符号つきの重み（経路積分 Monte Carlo の章）・自己正規化した重点サンプリングで現れます。分母が揺らぐので $O(1/N)$ の偏りがあります。

次の補題は、平均の比の偏りを評価する命題の準備です。比を $1/(1+\delta)$ の級数に展開すると、中心化した標本平均の積の期待値が現れます。この補題は、そのうち $N$ のどの冪で効くかを、同じ添字が 2 回以上現れる組の数え上げだけで抑えます。

補題 6.5（中心化した標本平均の積の期待値）✓

$W_i = (W_i^{(1)}, \ldots, W_i^{(q)})$（$i = 1, \ldots, N$）を独立同分布な $\mathbb{R}^q$ 値確率変数とし、各成分が $\mathbb{E}[W_i^{(l)}] = 0$、$|W_i^{(l)}| \le B$ を満たすとする。$\bar W^{(l)} = N^{-1}\sum_iW_i^{(l)}$ について

$$\Bigl|\mathbb{E}\Bigl[\prod_{l=1}^q\bar W^{(l)}\Bigr]\Bigr| \le q^qB^qN^{-\lceil q/2\rceil}.$$

証明（補題 6.5）

積を展開すると $\mathbb{E}[\prod_l\bar W^{(l)}] = N^{-q}\sum_{(i_1, \ldots, i_q)}\mathbb{E}[W_{i_1}^{(1)}\cdots W_{i_q}^{(q)}]$ で、和は $\{1, \ldots, N\}^q$ の全体にわたります。ある値 $i$ が $(i_1, \ldots, i_q)$ の中にちょうど 1 回だけ現れる（$i_l = i$ となる $l$ が 1 つ）とき、独立性から期待値は $\mathbb{E}[W_i^{(l)}] = 0$ と残りの因子の期待値の積になり、$0$ です。したがって残るのは、現れる値がどれも 2 回以上現れる組だけで、そのような組に現れる相異なる値は $\lfloor q/2\rfloor$ 個以下です。相異なる値が $j$ 個以下の組の個数は、値の集合の選び方 $\binom Nj \le N^j$ と各位置への割り当て $j^q$ の積以下なので、残る組は $N^{\lfloor q/2\rfloor}q^q$ 個以下です。各項の絶対値は $B^q$ 以下なので、全体は $N^{-q}\cdot N^{\lfloor q/2\rfloor}q^qB^q = q^qB^qN^{-\lceil q/2\rceil}$ 以下です。

∎

命題 6.6（比の推定量の偏り）✓

$(X_i, Y_i)$（$i = 1, 2, \ldots$）を独立同分布な組とし、定数 $0 < c \le B$ があって $|X_i| \le B$、$c \le Y_i \le B$ とする。$\mu_x = \mathbb{E}X_1$、$\mu_y = \mathbb{E}Y_1$、$\sigma_y^2 = \mathrm{Var}(Y_1)$、$\sigma_{xy} = \mathrm{Cov}(X_1, Y_1)$、$r = \mu_x/\mu_y$、$\hat R_N = \bar X_N/\bar Y_N$ とおく。このとき

$$\mathbb{E}[\hat R_N] = r + \frac{c_1}{N} + \frac{c_2}{N^2} + O(N^{-3}), \qquad c_1 = \frac{\mu_x\sigma_y^2}{\mu_y^3} - \frac{\sigma_{xy}}{\mu_y^2}$$

が成り立つ。$c_2$ は $N$ によらない定数で、証明の中に式を与える。したがってジャックナイフ推定値 $\hat R_{\mathrm J} = N\hat R_N - (N-1)\bar R_{(\cdot)}$ は $\mathbb{E}[\hat R_{\mathrm J}] = r + O(N^{-2})$ を満たす。

証明（命題 6.6）

$D_i = X_i - rY_i$、$Z_i = (Y_i - \mu_y)/\mu_y$ とおくと、どちらも平均 $0$ で有界です（$|D_i| \le B(1 + |r|)$、$|Z_i| \le B/\mu_y$）。$\bar D = \bar X_N - r\bar Y_N$、$\delta = \bar Z = \bar Y_N/\mu_y - 1$ とすると

$$\hat R_N - r = \frac{\bar X_N - r\bar Y_N}{\bar Y_N} = \frac{\bar D}{\mu_y(1 + \delta)}, \qquad 1 + \delta = \frac{\bar Y_N}{\mu_y} \ge \frac{c}{\mu_y} > 0$$

です。有限等比級数の和 $\sum_{k=0}^4(-\delta)^k = (1 + \delta^5)/(1 + \delta)$ から

$$\frac{1}{1 + \delta} = 1 - \delta + \delta^2 - \delta^3 + \delta^4 - \frac{\delta^5}{1 + \delta}$$

なので、

$$\mu_y\,\mathbb{E}[\hat R_N - r] = \mathbb{E}[\bar D] - \mathbb{E}[\bar D\delta] + \mathbb{E}[\bar D\delta^2] - \mathbb{E}[\bar D\delta^3] + \mathbb{E}[\bar D\delta^4] - \mathbb{E}\Bigl[\frac{\bar D\delta^5}{1+\delta}\Bigr]$$

です。各項を計算します。$a_1 = \mathbb{E}[D_1Z_1]$、$a_2 = \mathbb{E}[D_1Z_1^2]$、$a_3 = \mathbb{E}[D_1Z_1^3]$、$s^2 = \mathbb{E}[Z_1^2]$ とおきます。補題 6.5 の証明と同じく、添字の組のうち期待値が残るのは、どの値も 2 回以上現れるものだけです。

• $\mathbb{E}[\bar D] = 0$。
• $\mathbb{E}[\bar D\delta] = N^{-2}\sum_{i,j}\mathbb{E}[D_iZ_j] = N^{-2}\cdot N a_1 = a_1/N$（$i = j$ の項だけが残る）。
• $\mathbb{E}[\bar D\delta^2] = N^{-3}\cdot Na_2 = a_2/N^2$（3 つの添字がすべて等しい項だけが残る）。
• $\mathbb{E}[\bar D\delta^3] = N^{-4}\bigl(Na_3 + 3N(N-1)a_1s^2\bigr)$。残るのは 4 つの添字がすべて等しい $N$ 項と、$D$ の添字が 3 つの $Z$ のうち 1 つ（3 通り）と等しく、残りの 2 つの $Z$ が別の同じ添字をもつ $3N(N-1)$ 項です。よって $\mathbb{E}[\bar D\delta^3] = 3a_1s^2/N^2 + (a_3 - 3a_1s^2)/N^3$。
• $|\mathbb{E}[\bar D\delta^4]| \le 5^5B'^5N^{-3}$（補題 6.5 を $q = 5$、$B' = \max(B(1 + |r|), B/\mu_y)$ で使う）。
• 最後の項は、$1 + \delta \ge c/\mu_y$ と、$L^2$ の内積 $\mathbb{E}[UV]$ についての Cauchy–Schwarz の不等式（命題 2.2［Hilbert 空間の基礎］）から $$\Bigl|\mathbb{E}\Bigl[\frac{\bar D\delta^5}{1 + \delta}\Bigr]\Bigr| \le \frac{\mu_y}{c}\sqrt{\mathbb{E}[\bar D^2]}\sqrt{\mathbb{E}[\delta^{10}]} \le \frac{\mu_y}{c}\sqrt{2^2B'^2N^{-1}}\sqrt{10^{10}B'^{10}N^{-5}} = O(N^{-3})$$ です（補題 6.5 を $q = 2$ と $q = 10$ で使う）。

まとめると

$$\mathbb{E}[\hat R_N] - r = -\frac{a_1}{\mu_yN} + \frac{a_2 - 3a_1s^2}{\mu_yN^2} + O(N^{-3})$$

で、$c_2 = (a_2 - 3a_1s^2)/\mu_y$ です。$c_1$ は、$a_1 = \mathrm{Cov}(X_1 - rY_1, Y_1)/\mu_y = (\sigma_{xy} - r\sigma_y^2)/\mu_y$ から

$$c_1 = -\frac{a_1}{\mu_y} = \frac{r\sigma_y^2 - \sigma_{xy}}{\mu_y^2} = \frac{\mu_x\sigma_y^2}{\mu_y^3} - \frac{\sigma_{xy}}{\mu_y^2}$$

です。

ジャックナイフについて。$\hat R_{(i)}$ は $N - 1$ 個の独立同分布な組から作った比なので $\mathbb{E}[\bar R_{(\cdot)}] = \mathbb{E}[\hat R_{N-1}]$ です。上の展開は 命題 6.7［格子場のモンテカルロ法］ の仮定（$\mathbb{E}f(\bar A_n) = f(\mu) + c/n + d/n^2 + O(n^{-3})$）を満たすので、$\mathbb{E}[\hat R_{\mathrm J}] = r + O(N^{-2})$ です。計算を書くと

$$N\Bigl(\frac{c_1}{N} + \frac{c_2}{N^2}\Bigr) - (N-1)\Bigl(\frac{c_1}{N-1} + \frac{c_2}{(N-1)^2}\Bigr) = \frac{c_2}{N} - \frac{c_2}{N-1} = -\frac{c_2}{N(N-1)}$$

で、$1/N$ の項が消え、$O(N^{-3})$ の剰余は $N$ 倍しても $O(N^{-2})$ です。

∎

偏り $c_1/N$ は標準誤差 $O(N^{-1/2})$ より小さいので、1 回の計算では見えません。問題になるのは、短い計算を多数回行って平均するときです。$N = 20$ の比を $10^4$ 個平均すると、統計誤差は $1/100$ になりますが偏りは $c_1/20$ のままで、偏りが誤差棒を超えます。並列計算で短い連鎖を多数走らせるときに起きる誤りです。§9 でこの命題の $c_1$・$c_2$ を数値と照合します。

6.3. 相関のあるデータでの使い方

定理 6.3 と 命題 6.6 は独立同分布な標本についての主張です。相関のある $N$ 個のデータから 1 個ずつ除くジャックナイフは、データを独立と見なした誤差（$\sqrt{2\tau_{\mathrm{int}}}$ 倍小さい）を返します。そこで、まず長さ $b$ のバッチ平均 $Y_1, \ldots, Y_k$ を作り、バッチを 1 個ずつ除くジャックナイフを使います。定理 5.2 から、$b$ が相関の長さより十分大きければバッチ平均の分散は $\sigma_\infty^2/b$ に近く、バッチ間の相関は端の $O(\tau_{\mathrm{int}}/b)$ の部分だけに由来します。バッチ長についての注意は §5 と同じで、$b$ を 2 倍ずつ増やして誤差が頭打ちになることを確かめます。

比熱 $\propto \langle E^2\rangle - \langle E\rangle^2$ は $(E, E^2)$ の平均の関数 $g(u, v) = v - u^2$、有効質量 $\ln(C(t)/C(t+1))$ は相関関数の平均の関数です。どちらも $g$ が滑らかなので、バッチ化したあとのジャックナイフが使えます。有効質量では $C(t)$ と $C(t+1)$ が強く正に相関していて、比の誤差は個々の誤差よりずっと小さくなります。ジャックナイフは同じバッチを同時に除くので、この相関が自動的に入ります。

7. ブートストラップ

7.1. 再標本化分布

定義 7.1（ブートストラップ（再標本化分布））↗

データ $x_1, \ldots, x_N$ を固定する。$I_1, \ldots, I_N$ を $\{1, \ldots, N\}$ 上の一様分布に従う独立な確率変数とし（その確率を $\mathbb{P}^*$ と書く）、$X_k^* = x_{I_k}$ を再標本という。統計量 $T$ について、$\mathbb{P}^*$ のもとでの $T(X_1^*, \ldots, X_N^*)$ の分布を $T$ の再標本化分布（ブートストラップ分布）という。

再標本は、経験分布（各 $x_i$ に確率 $1/N$）からの独立な $N$ 個の標本です。手元のデータを「母集団の縮図」と見なし、そこから何度も取り直したときの統計量のばらつきを、本物の標本を取り直したときのばらつきの代わりに使う、というのが誤差が分かる理由です。再標本化分布はデータの関数で、$N^N$ 通りの $(I_1, \ldots, I_N)$ を数え上げれば原理的には厳密に求まります。実際には $\mathbb{P}^*$ からの Monte Carlo（再標本を数百〜数千組作る）で近似します。「未知の分布 $F$ のもとでの $T$ の分布」を「経験分布 $F_N$ のもとでの $T$ の分布」で代用する、というのがブートストラップの考え方です。代用が正しいかどうかは統計量によります。

定義 7.2（ブロック・ブートストラップ）↗

相関のあるデータ $A_1, \ldots, A_N$ と $N = kb$ について、長さ $b$ の連続した区間 $(A_s, \ldots, A_{s+b-1})$（$s = 1, \ldots, N - b + 1$）から $k$ 個を一様に独立に選び、選んだ順に並べて長さ $N$ の再標本を作る方法を**（移動）ブロック・ブートストラップ**という。

データを 1 個ずつ選び直すと時間方向の相関が壊れます。ブロックごとに選べば、ブロックの中の相関が保たれます。$b$ の条件はバッチ平均と同じ（$b \to \infty$ かつ $N/b \to \infty$）です。移動ブロック法は H. R. Künsch, Ann. Statist. 17 (1989)、ブロック長を幾何分布で無作為にして再標本を定常にする定常ブートストラップは D. N. Politis, J. P. Romano, J. Amer. Statist. Assoc. 89 (1994) にあります。

7.2. 標本平均での一致性

2 つの分布関数 $F, G$ の Kolmogorov 距離を $\sup_{x\in\mathbb{R}}|F(x) - G(x)|$ と定めます。

定理 7.3（標本平均のブートストラップの一致性）✓

$X_1, X_2, \ldots$ を独立同分布な実確率変数とし、$\mu = \mathbb{E}X_1$、$\sigma^2 = \mathrm{Var}(X_1) > 0$、$\rho_3 = \mathbb{E}|X_1 - \mu|^3 < \infty$ とする。$\hat\sigma_N^2 = N^{-1}\sum_{i\le N}(X_i - \bar X_N)^2$ とおき、

$$F_N(x) = \mathbb{P}\Bigl(\frac{\sqrt N(\bar X_N - \mu)}{\sigma} \le x\Bigr), \qquad H_N(x) = \mathbb{P}^*\Bigl(\frac{\sqrt N(\bar X_N^* - \bar X_N)}{\hat\sigma_N} \le x\Bigr)$$

とする（$H_N$ はデータ $X_1, \ldots, X_N$ で決まる分布関数で、$\hat\sigma_N > 0$ のとき定義する。$\bar X_N^*$ は再標本の平均）。$C$ を Berry–Esseen の定理の絶対定数とすると、確率 $1$ で

$$\limsup_{N\to\infty}\sqrt N\sup_{x\in\mathbb{R}}|H_N(x) - F_N(x)| \le \frac{5C\rho_3}{\sigma^3}$$

が成り立つ。特に Kolmogorov 距離は確率 $1$ で $0$ に収束する。

証明（定理 7.3）

Berry–Esseen の定理（定理 5.9［大数の法則と中心極限定理］）を $(X_i)$ に使うと、すべての $N$ で

$$\sup_x|F_N(x) - \Phi(x)| \le \frac{C\rho_3}{\sigma^3\sqrt N}$$

です。次にデータ $X_1 = x_1, \ldots, X_N = x_N$ を固定し、$\mathbb{P}^*$ のもとで考えます。$X_k^* = x_{I_k}$（$k = 1, 2, \ldots$。添字の列を無限に延ばしておきます）は独立同分布で、

$$\mathbb{E}^*[X_1^*] = \frac1N\sum_ix_i = \bar x, \qquad \mathrm{Var}^*(X_1^*) = \frac1N\sum_i(x_i - \bar x)^2 = \hat\sigma_N^2, \qquad \mathbb{E}^*|X_1^* - \bar x|^3 = \frac1N\sum_i|x_i - \bar x|^3 = \hat\rho_N$$

です。$\hat\sigma_N > 0$ のとき、同じ定理を列 $(X_k^*)$ に、標本の個数 $N$ で使うと

$$\sup_x|H_N(x) - \Phi(x)| \le \frac{C\hat\rho_N}{\hat\sigma_N^3\sqrt N}$$

です。三角不等式から $\sqrt N\sup_x|H_N - F_N| \le C(\hat\rho_N/\hat\sigma_N^3 + \rho_3/\sigma^3)$ です。

$\hat\sigma_N$ と $\hat\rho_N$ の挙動を調べます。大数の強法則（定理 4.5［大数の法則と中心極限定理］）を $X_i$、$X_i^2$、$|X_i - \mu|^3$ に使うと、確率 $1$ で $\bar X_N \to \mu$、$\hat\sigma_N^2 = N^{-1}\sum X_i^2 - \bar X_N^2 \to \sigma^2 > 0$、$N^{-1}\sum_i|X_i - \mu|^3 \to \rho_3$ です。$u \mapsto |u|^3$ は凸なので $|a + b|^3 = 8|(a+b)/2|^3 \le 4(|a|^3 + |b|^3)$ で、$a = X_i - \mu$、$b = \mu - \bar X_N$ とすると

$$\hat\rho_N \le 4\Bigl(\frac1N\sum_i|X_i - \mu|^3 + |\bar X_N - \mu|^3\Bigr) \to 4\rho_3$$

です。したがって確率 $1$ で、十分大きい $N$ では $\hat\sigma_N > 0$ で $H_N$ が定義され、$\limsup_N\hat\rho_N/\hat\sigma_N^3 \le 4\rho_3/\sigma^3$ です。合わせて結論を得ます。

∎

この定理から、$\bar X_N^* - \bar X_N$ の再標本化分布の分位点を使って $\mu$ の信頼区間が作れます。$q_\alpha^*$ を $\bar X_N^*$ の再標本化分布の $\alpha$ 分位点とすると、区間 $[2\bar X_N - q_{1-\alpha/2}^*,\ 2\bar X_N - q_{\alpha/2}^*]$ の被覆確率は $1 - \alpha$ に収束します（$\hat\sigma_N \to \sigma$ と $\Phi$ の連続性から、標準化を外した分布どうしの Kolmogorov 距離も $0$ に向かうためです）。

注意 7.4（ブートストラップが正規近似より良くなる場合）

定理 7.3 の速さ $O(N^{-1/2})$ は正規近似と同じです。標本標準偏差で割った統計量 $\sqrt N(\bar X_N - \mu)/s_N$ については、分布が格子分布でなければ、再標本化分布との Kolmogorov 距離が $o(N^{-1/2})$ になります。歪度による $N^{-1/2}$ の補正項（Edgeworth 展開の第 1 項）までブートストラップが再現するためです。証明は Edgeworth 展開によるもので、K. Singh, Ann. Statist. 9 (1981) にあります。この記事では定理にしません。

7.3. 最大値での破綻

命題 7.5（一様分布の最大値ではブートストラップは一致しない）✓

$\theta > 0$ とし、$X_1, \ldots, X_N$ を $(0, \theta)$ 上の一様分布からの独立な標本、$M_N = \max_iX_i$、$M_N^* = \max_kX_k^*$（再標本の最大値）とする。

1. データの値がすべて異なるとき（確率 $1$）、$\mathbb{P}^*(M_N^* = M_N) = 1 - (1 - 1/N)^N \ge 1 - 1/e$。
2. $x \ge 0$ について $\mathbb{P}(N(\theta - M_N) > x) = (1 - x/(N\theta))^N$（$x \le N\theta$）で、$N \to \infty$ で $e^{-x/\theta}$ に収束する。特に $\mathbb{P}(N(\theta - M_N) \le 0) = 0$。
3. したがって、$N(\theta - M_N)$ の分布関数と、その再標本化による代用 $N(M_N - M_N^*)$ の $\mathbb{P}^*$ のもとでの分布関数の Kolmogorov 距離は、すべての $N$ で確率 $1$ で $1 - 1/e$ 以上である。

証明（命題 7.5）

1. 値がすべて異なるとき、最大値を与える添字 $i_0$ はただ 1 つです。$M_N^* = M_N$ は、$I_1, \ldots, I_N$ のうち少なくとも 1 つが $i_0$ であることと同値です。どれも $i_0$ でない確率は $(1 - 1/N)^N$ です。$\ln(1 - u) \le -u$（$0 \le u < 1$）から $(1 - 1/N)^N \le e^{-1}$ です。値がすべて異なる確率が $1$ であることは、$i \ne j$ について $(X_i, X_j)$ が正方形上の一様分布に従い、対角線の面積が $0$ であることから分かります。

2. $0 \le x \le N\theta$ で $\mathbb{P}(N(\theta - M_N) > x) = \mathbb{P}(X_i < \theta - x/N\ \forall i) = (1 - x/(N\theta))^N$ です。極限は $(1 - u/N)^N \to e^{-u}$ によります。$x = 0$ とすると $\mathbb{P}(N(\theta - M_N) > 0) = 1$ です。

3. $x = 0$ での 2 つの分布関数の値は、1 から $\mathbb{P}^*(N(M_N - M_N^*) \le 0) = \mathbb{P}^*(M_N^* = M_N) \ge 1 - 1/e$（$M_N^* \le M_N$ に注意）、2 から $\mathbb{P}(N(\theta - M_N) \le 0) = 0$ です。差は $1 - 1/e$ 以上です。

∎

再標本の最大値が、元データの上から j+1 番目の値になる確率の極限 e^{-j}(1-e^{-1})。最大値と一致する確率 0.632 の山があり、真の分布（連続）とは形が違う。

真の分布（指数分布）は連続なのに、再標本化分布は $0$ に確率 $0.632$ の原子をもちます。平均とのちがいは、最大値が経験分布の「滑らかな汎関数」でないことです。平均は全データの寄与の和ですが、最大値は 1 個のデータで決まり、再標本にそのデータが入るかどうかで値が跳びます。端点・極値・分位点の端・しきい値を超える回数のような統計量にブートストラップを使うときは、この種の破綻を疑ってください。再標本の大きさを $m \ll N$ に減らす「$N$ 個から $m$ 個」のブートストラップはこの例を直しますが、$m$ の選び方という新しい問題が生じます（P. J. Bickel, D. A. Freedman, Ann. Statist. 9 (1981) に反例の一般論があります）。

8. 熱化と収束の診断

8.1. 初期分布による偏り

Markov 連鎖 Monte Carlo では、連鎖を定常分布から始めることができません（それができるなら連鎖は要りません）。初期の標本は定常分布からずれています。

定義 8.1（熱化）↗

連鎖 $X_0, X_1, \ldots$ のはじめの $n_0$ 個を捨て、$X_{n_0}, \ldots, X_{n_0+N-1}$ だけで時間平均を作ることを熱化（burn-in）といい、$n_0$ を熱化の長さという。

命題 8.2（初期分布による時間平均の偏り）✓

§2.2 の設定で絶対スペクトルギャップが $\gamma_* > 0$ とし、$\pi_{\min} = \min_x\pi(x)$ とする。初期分布を任意の $\nu$ とした連鎖について、$f : S \to \mathbb{R}$ と $n_0 \ge 0$、$N \ge 1$ に対し

$$\Bigl|\mathbb{E}_\nu\Bigl[\frac1N\sum_{k=n_0}^{n_0+N-1}f(X_k)\Bigr] - \pi(f)\Bigr| \le \sqrt{\frac{1}{\pi_{\min}} - 1}\;\frac{(1-\gamma_*)^{n_0}}{\gamma_*N}\sqrt{\mathrm{Var}_\pi(f)}.$$

証明（命題 8.2）

$\tilde f = f - \pi(f)$、$h(x) = \nu(x)/\pi(x)$ とおきます。$\mathbb{P}_\nu(X_k = y) = (\nu P^k)(y)$（命題 2.3［Markov 連鎖と Poisson 過程］）から

$$\mathbb{E}_\nu[f(X_k)] - \pi(f) = \sum_x\nu(x)(P^k\tilde f)(x) = \langle h, P^k\tilde f\rangle_\pi$$

です。$\pi$ の定常性から $\pi(P^k\tilde f) = \pi(\tilde f) = 0$ なので $\langle 1, P^k\tilde f\rangle_\pi = 0$ で、$\langle h, P^k\tilde f\rangle_\pi = \langle h - 1, P^k\tilde f\rangle_\pi$ です。内積の Cauchy–Schwarz の不等式（定理 4.2［内積空間とグラム・シュミット直交化］）から

$$|\langle h - 1, P^k\tilde f\rangle_\pi| \le \|h - 1\|_\pi\|P^k\tilde f\|_\pi$$

です。補題 2.3 の展開 $\tilde f = \sum_{j\ge2}\langle\tilde f, \varphi_j\rangle_\pi\varphi_j$ から $\|P^k\tilde f\|_\pi^2 = \sum_{j\ge2}\lambda_j^{2k}\langle\tilde f, \varphi_j\rangle_\pi^2 \le (1-\gamma_*)^{2k}\|\tilde f\|_\pi^2$ です。また

$$\|h - 1\|_\pi^2 = \sum_x\frac{\nu(x)^2}{\pi(x)} - 2\sum_x\nu(x) + 1 = \sum_x\frac{\nu(x)^2}{\pi(x)} - 1 \le \frac{1}{\pi_{\min}}\sum_x\nu(x) - 1 = \frac{1}{\pi_{\min}} - 1$$

です（$\nu(x)^2 \le \nu(x)$）。$k = n_0, \ldots, n_0 + N - 1$ について平均すると、$\sum_{k\ge n_0}(1-\gamma_*)^k = (1-\gamma_*)^{n_0}/\gamma_*$ から結論を得ます。

∎

熱化なし（$n_0 = 0$）でも偏りは $O(1/(\gamma_*N))$ です。統計誤差は $\sqrt{2\tau_{\mathrm{int}}\mathrm{Var}_\pi(f)/N}$ で $N^{-1/2}$ の大きさなので、$N$ が緩和時間 $1/\gamma_*$ より十分長ければ、初期分布の偏りは統計誤差より高次です。それでも熱化を行う理由は 2 つあります。第 1 に、係数 $\sqrt{1/\pi_{\min} - 1}$ が巨大になりうること（初期配位が定常分布のもとで極端に珍しい場合）。$n_0$ を $\ln(1/\pi_{\min})/(2\gamma_*)$ の数倍に取れば因子 $(1-\gamma_*)^{n_0}$ がこれを打ち消します。第 2 に、初期の過渡的な部分が分散の推定を汚すこと。平衡値から大きく外れた初期の値は、自己共分散やバッチ平均の分散の推定値を膨らませ、誤差棒と $\tau_{\mathrm{int}}$ の推定を狂わせます。熱化は平均の偏りを消すためというより、誤差の推定を守るために行うものです。

捨てすぎにも代償があります。$n_0$ 個捨てると統計誤差は $\sqrt{(N + n_0)/N}$ 倍になります。緩和時間の 10〜20 倍を捨てれば上の因子は $e^{-10}$ 以下になるので、全長の半分を捨てるような習慣は、連鎖が緩和時間に比べて十分長い限り無駄です。

8.2. 複数の連鎖による診断

1 本の連鎖の内部だけを見ても、連鎖が状態空間の一部に閉じ込められていることは分かりません。初期値をばらつかせた複数の連鎖を比べるのが、Gelman と Rubin の診断です。

定義 8.3（Gelman–Rubin の統計量）↗

$m \ge 2$ 本の連鎖のそれぞれから長さ $n \ge 2$ のデータ $A_{j,1}, \ldots, A_{j,n}$（$j = 1, \ldots, m$）を得たとする。連鎖ごとの平均 $\bar A_j$、標本分散 $s_j^2 = (n-1)^{-1}\sum_t(A_{j,t} - \bar A_j)^2$、全体の平均 $\bar A = m^{-1}\sum_j\bar A_j$ について

$$W = \frac1m\sum_{j=1}^ms_j^2, \qquad \frac Bn = \frac{1}{m-1}\sum_{j=1}^m(\bar A_j - \bar A)^2, \qquad \hat V = \frac{n-1}{n}W + \frac Bn, \qquad \hat R = \sqrt{\frac{\hat V}{W}}$$

とおく。$W$ を連鎖内分散、$B/n$ を連鎖間分散、$\hat R$ を Gelman–Rubin の統計量（potential scale reduction factor）という。

命題 8.4（R̂ の収束）✓

定義 8.3 で、定数 $a$、$a_2$（$\sigma^2 = a_2 - a^2 > 0$）があって、どの連鎖 $j$ についても確率 $1$ で

$$\frac1n\sum_{t=1}^nA_{j,t} \to a, \qquad \frac1n\sum_{t=1}^nA_{j,t}^2 \to a_2 \qquad (n \to \infty)$$

が成り立つとする。このとき確率 $1$ で $W \to \sigma^2$、$B/n \to 0$、$\hat R \to 1$ である。

証明（命題 8.4）

$s_j^2 = \frac{n}{n-1}\bigl(n^{-1}\sum_tA_{j,t}^2 - \bar A_j^2\bigr)$ なので、仮定から確率 $1$ で $s_j^2 \to a_2 - a^2 = \sigma^2$ です。有限個の $j$ についての平均 $W$ も $\sigma^2$ に収束します。$\bar A_j \to a$ がすべての $j$ で成り立つので $\bar A \to a$、$B/n = (m-1)^{-1}\sum_j(\bar A_j - \bar A)^2 \to 0$ です。よって $\hat V \to \sigma^2$ で、$\sigma^2 > 0$ から $\hat R^2 = \hat V/W \to 1$ です。

∎

有限・既約な連鎖では、任意の初期分布について時間平均が $\pi(f)$ に概収束するので（Markov 連鎖 Monte Carlo の章のエルゴード定理 定理 3.2［Markov 連鎖 Monte Carlo］ を $f$ と $f^2$ に使います）、仮定は $a = \pi(f)$、$a_2 = \pi(f^2)$ で満たされます。$\hat V$ は、連鎖内分散と連鎖間分散を足して全体の分散 $\mathrm{Var}_\pi(f)$ を見積もったものです。連鎖が別々の領域に留まっていると、$W$ は領域内の分散しか見ないので小さく、$B/n$ は領域間の差で大きくなり、$\hat R > 1$ になります。実務では $\hat R < 1.01$ 程度を目安にします。

ただし逆は成り立ちません。

例 8.5（R̂ ≈ 1 でも収束していない例）

$S = \{1, 2, 3, 4\}$ とし、$\{1, 2\}$ の中と $\{3, 4\}$ の中でそれぞれ確率 $1/2$ ずつで動き、2 つの組の間は行き来しない遷移行列 $P$ を考えます。一様分布 $\pi$ は $P$ について可逆な定常分布で、$f(x) = x$ の平均は $\pi(f) = 2.5$ です。すべての連鎖を $\{1, 2\}$ の中から始めると、各連鎖の $f(X_t)$ は $1$ と $2$ を確率 $1/2$ で取る独立な列で、大数の強法則から時間平均は $1.5$、2 乗の時間平均は $2.5$ に概収束します。命題 8.4 の仮定が $a = 1.5$、$\sigma^2 = 0.25$ で満たされるので $\hat R \to 1$ ですが、時間平均は $\pi(f) = 2.5$ に収束していません。

この $P$ は既約でない極端な例です。2 つの組の間を確率 $\varepsilon$ で行き来するように直すと既約になりますが、$1/\varepsilon$ より短い計算では上と同じことが起きます。2 つの山をもつ分布で、すべての連鎖を同じ山から始めた場合がこれに当たります。$\hat R$ は「連鎖どうしが食い違っている」ことを検出する道具で、「どの連鎖も見ていない領域がある」ことは検出できません。初期値を目標分布より広くばらつかせることが、この診断の前提です。

9. 数値実験

コードは Python 3.9.6・NumPy 2.0.2・SciPy 1.13.1 で実行し、出力をそのまま載せています。乱数は numpy.random.default_rng(20261008) です。

9.1. AR(1) 過程での τ_int の推定

例 5.4 の過程（$\varrho = 0.9$、$\tau_{\mathrm{int}} = 9.5$）を $N = 10^5$ で作り、窓つき推定量（$W = 20, 50, 200, 1000$）、初期正値列、バッチ平均（$b = 16, \ldots, 4096$）を比べます。独立な 100 本について繰り返し、推定量の平均とばらつきを見ます。

import numpy as np
from scipy.signal import lfilter

def ar1(rng, n, rho):
    # 定常な AR(1): x_{t+1} = rho x_t + sqrt(1-rho^2) z_t, 分散 1
    z = rng.standard_normal(n)
    z[1:] *= np.sqrt(1.0 - rho * rho)
    return lfilter([1.0], [1.0, -rho], z)

def autocov(x, tmax):
    n = len(x)
    d = x - x.mean()
    f = np.fft.rfft(d, 2 * n)
    g = np.fft.irfft(f * np.conj(f))[: tmax + 1]
    return g / (n - np.arange(tmax + 1))

def tau_window(g, w):
    return 0.5 + g[1 : w + 1].sum() / g[0]

def tau_ips(g):
    # 初期正値列: G_k = g(2k) + g(2k+1) が正である間だけ足す
    s = 0.0
    for k in range(len(g) // 2):
        pair = g[2 * k] + g[2 * k + 1]
        if pair <= 0.0:
            break
        s += pair
    return s / g[0] - 0.5

def tau_batch(x, b):
    k = len(x) // b
    y = x[: k * b].reshape(k, b).mean(axis=1)
    return b * y.var(ddof=1) / (2.0 * x.var(ddof=1))

rho, n, reps = 0.9, 100_000, 100
rng = np.random.default_rng(20261008)
names = ["W=20", "W=50", "W=200", "W=1000", "IPS"] + [f"b={b}" for b in (16, 64, 256, 1024, 4096)]
out = np.empty((reps, len(names)))
naive = np.empty(reps); means = np.empty(reps)
for r in range(reps):
    x = ar1(rng, n, rho)
    g = autocov(x, 2000)
    row = [tau_window(g, w) for w in (20, 50, 200, 1000)] + [tau_ips(g)]
    row += [tau_batch(x, b) for b in (16, 64, 256, 1024, 4096)]
    out[r] = row
    means[r] = x.mean(); naive[r] = x.std(ddof=1) / np.sqrt(n)
print("true tau_int = 9.5")
for j, name in enumerate(names):
    print(f"{name:>7s}  first={out[0, j]:7.3f}  mean={out[:, j].mean():7.3f}  sd={out[:, j].std(ddof=1):6.3f}")
print(f"sd of sample means      = {means.std(ddof=1):.5f}")
print(f"mean naive std. error   = {naive.mean():.5f}")
print(f"ratio = {means.std(ddof=1) / naive.mean():.3f},  sqrt(2 tau) = {np.sqrt(19.0):.3f}")

true tau_int = 9.5
   W=20  first=  8.174  mean=  8.401  sd= 0.117
   W=50  first=  9.164  mean=  9.477  sd= 0.328
  W=200  first=  9.030  mean=  9.579  sd= 0.779
 W=1000  first=  7.470  mean=  9.508  sd= 1.931
    IPS  first=  9.183  mean=  9.661  sd= 0.471
   b=16  first=  4.874  mean=  4.918  sd= 0.038
   b=64  first=  7.777  mean=  8.097  sd= 0.215
  b=256  first=  8.714  mean=  9.254  sd= 0.592
 b=1024  first=  8.695  mean=  9.406  sd= 1.389
 b=4096  first=  4.558  mean=  9.299  sd= 3.248
sd of sample means      = 0.01437
mean naive std. error   = 0.00316
ratio = 4.541,  sqrt(2 tau) = 4.359

first は 1 本目の値、mean と sd は 100 本の平均と標準偏差です（mean の標準誤差は sd の $1/10$）。

• 窓の偏り。$W = 20$ の平均 $8.401$ は、命題 4.1 の予言 $9.5 - 1.094 = 8.406$ と標準誤差 $0.012$ の範囲で一致します。$W \ge 50$ では裾 $0.9^{W+1}/0.1$ が $0.05$ 以下で、偏りは見えません。
• 窓の分散。Madras–Sokal の公式（注意 4.2）の $\tau_{\mathrm{int}}\sqrt{2(2W+1)/N}$ は $W = 200$ で $0.85$、$W = 1000$ で $1.90$ で、実測の $0.78$、$1.93$ と合います。$W = 1000$ では 1 本の推定値が $7.47$ と、真値から $2$ 以上ずれています。窓を大きく取りすぎた代償です。
• 初期正値列は窓を自動で選び、ばらつき $0.47$ は $W = 50$ と $W = 200$ の間です。平均 $9.66$ は真値より標準誤差の $3.4$ 倍だけ大きく、小さな正の偏りが見えます（対の和が正である限り足し続ける規則なので、裾の正の雑音が負の雑音より多く和に入るためと考えられます）。
• バッチ平均。$b = 16, 64$ の平均は 例 5.4 の式の値 $4.917$、$8.095$ と一致し、偏りが $1/b$ で減っています。$b = 4096$ ではバッチが $24$ 個しかなく、ばらつきが $3.2$ で、1 本目の値 $4.56$ は真値の半分以下です（定理 5.2 の 3 と同じく、バッチ数が少ないことによるばらつきです）。
• 相関を無視した標準誤差は $0.00316$ で、標本平均の実際のばらつき $0.01437$（100 本からの推定なので相対誤差は約 $7\%$）の $1/4.5$ です。理論値は $1/\sqrt{2\tau_{\mathrm{int}}} = 1/4.36$ です。

9.2. 比の推定量の偏り

$Y \sim \mathrm U(0.2, 1.8)$、$X = Y + V$（$V \sim \mathrm U(0,1)$ は $Y$ と独立）とします。$\mu_x = 1.5$、$\mu_y = 1$、$r = 1.5$、$\sigma_y^2 = \sigma_{xy} = 1.6^2/12$ で、命題 6.6 の係数は $c_1 = 0.5\sigma_y^2 = 0.1067$ です。$D = (V - 1/2) - (Y - 1)/2$、$Z = Y - 1$ から $a_1 = -\sigma_y^2/2$、$a_2 = -\mathbb{E}[(Y-1)^3]/2 = 0$ で、$c_2 = -3a_1s^2 = 1.5\sigma_y^4 = 0.0683$ です。

import numpy as np

rng = np.random.default_rng(20261008)
s2 = 1.6**2 / 12                  # Var(Y), Y ~ U(0.2, 1.8), mu_y = 1
c1 = 1.5 * s2 - s2                # mu_x s_y^2/mu_y^3 - s_xy/mu_y^2
c2 = 3 * (0.5 * s2) * s2          # (a_2 - 3 a_1 s^2)/mu_y, a_1 = -s2/2, a_2 = 0
reps, chunk = 16_000_000, 200_000
print("   N  bias(ratio)   +-se     c1/N    c1/N+c2/N^2  bias(jack)    +-se    -c2/(N(N-1))")
for n in (20, 40, 80):
    plain, jack = [], []
    for _ in range(reps // chunk):
        y = rng.uniform(0.2, 1.8, size=(chunk, n))
        x = y + rng.uniform(0.0, 1.0, size=(chunk, n))
        sx, sy = x.sum(axis=1, keepdims=True), y.sum(axis=1, keepdims=True)
        r = (sx / sy)[:, 0]
        loo = ((sx - x) / (sy - y)).mean(axis=1)
        plain.append(r - 1.5)
        jack.append(n * r - (n - 1) * loo - 1.5)
    p, j = np.concatenate(plain), np.concatenate(jack)
    se = lambda v: v.std(ddof=1) / np.sqrt(len(v))
    print(f"{n:4d}  {p.mean():.6f}  {se(p):.6f}  {c1 / n:.6f}  {c1 / n + c2 / n**2:.6f}   "
          f"{j.mean():+.6f}  {se(j):.6f}  {-c2 / (n * (n - 1)):+.6f}")

   N  bias(ratio)   +-se     c1/N    c1/N+c2/N^2  bias(jack)    +-se    -c2/(N(N-1))
  20  0.005476  0.000021  0.005333  0.005504   -0.000226  0.000021  -0.000180
  40  0.002717  0.000015  0.002667  0.002709   -0.000038  0.000015  -0.000044
  80  0.001351  0.000010  0.001333  0.001344   -0.000004  0.000010  -0.000011

比の偏りは $c_1/N + c_2/N^2$ と標準誤差の $-1.3$、$+0.5$、$+0.7$ 倍の差で一致し、$N$ を 2 倍にすると半分になっています。ジャックナイフ補正後の偏りは $1/20$ 以下に減り、予言 $-c_2/(N(N-1))$ との差は標準誤差の $-2.2$、$+0.4$、$+0.7$ 倍です（$N = 20$ では、ここで落とした $O(N^{-3})$ の項が $10^{-5}$ の大きさで残ります）。

9.3. ブートストラップの区間の被覆率

指数分布の平均（真値 $1$）と、一様分布 $\mathrm U(0,1)$ の最大値（真値 $\theta = 1$）について、§7.2 の作り方の $95\%$ 区間 $[2\hat\theta - q_{0.975}^*, 2\hat\theta - q_{0.025}^*]$ が真値を含む割合を、2000 回の反復で求めます（再標本は 1000 組）。

import numpy as np

rng = np.random.default_rng(20261008)
nboot, reps = 1000, 2000
q = lambda v: np.quantile(v, [0.025, 0.975], method="inverted_cdf")
print("   N   cover(mean)  +-se    cover(max)  +-se    P*(M*=M_N)  1-(1-1/N)^N")
for n in (50, 200, 800):
    hit_mean = hit_max = 0
    atom = 0.0
    for _ in range(reps):
        idx = rng.integers(0, n, size=(nboot, n))
        # 平均: 指数分布 Exp(1)（真値 1）。区間 [2m - q_hi, 2m - q_lo]
        x = rng.exponential(1.0, n)
        m = x.mean()
        lo, hi = q(x[idx].mean(axis=1))
        hit_mean += (2 * m - hi <= 1.0 <= 2 * m - lo)
        # 最大値: U(0,1)（真値 theta = 1）。同じ作り方の区間
        u = rng.uniform(0.0, 1.0, n)
        mx = u.max()
        bm = u[idx].max(axis=1)
        lo, hi = q(bm)
        hit_max += (2 * mx - hi <= 1.0 <= 2 * mx - lo)
        atom += np.mean(bm == mx)
    p1, p2 = hit_mean / reps, hit_max / reps
    se = lambda p: np.sqrt(p * (1 - p) / reps)
    print(f"{n:4d}   {p1:.4f}     {se(p1):.4f}   {p2:.4f}     {se(p2):.4f}   {atom / reps:.4f}      {1 - (1 - 1 / n) ** n:.4f}")

   N   cover(mean)  +-se    cover(max)  +-se    P*(M*=M_N)  1-(1-1/N)^N
  50   0.9145     0.0063   0.8705     0.0075   0.6358      0.6358
 200   0.9395     0.0053   0.8810     0.0072   0.6326      0.6330
 800   0.9455     0.0051   0.8715     0.0075   0.6324      0.6324

平均の被覆率は $N$ とともに $0.95$ に近づきます（定理 7.3。$N = 50$ での不足は、指数分布の歪みによる $O(N^{-1/2})$ の誤差です）。最大値の被覆率は $N$ を増やしても $0.87$ 付近から動きません。再標本の最大値が標本の最大値に一致する割合は 命題 7.5 の $1 - (1 - 1/N)^N$ に一致しています。$0.87$ という値も説明できます。$\mathbb{P}^*(M_N^* \le X_{(N-j)}) = (1 - j/N)^N \approx e^{-j}$ で、$e^{-3} \approx 0.050 \ge 0.025 > e^{-4}$ なので下側の分位点は上から 4 番目の値 $X_{(N-3)}$ です。区間が $\theta$ を含む条件は $\theta - M_N \le M_N - X_{(N-3)}$ で、$N$ 倍すると左辺は平均 $1$ の指数分布、右辺はそれと独立な 3 個の指数分布の和に近づくので、確率は $1 - 2^{-3} = 0.875$ に近づきます。名目の $0.95$ にはなりません。

10. 演習

演習 10.1標準

例 5.4 の AR(1) 過程（$0 < \varrho < 1$）について、$\sigma_\infty^2 - b\,\mathrm{Var}(Y_1) = \dfrac{2\varrho(1 - \varrho^b)}{b(1-\varrho)^2}$ を示し、定理 5.2 の上界 $2C_1/b$ と比べよ。

解答

定理 5.2 の証明の等式から、$\Gamma(t) = \varrho^{|t|}$ として

$$\sigma_\infty^2 - b\,\mathrm{Var}(Y_1) = 2\sum_{t\ge b}\varrho^t + \frac2b\sum_{t=1}^{b-1}t\varrho^t$$

です。$\sum_{t\ge b}\varrho^t = \varrho^b/(1-\varrho)$。また $\sum_{t=1}^{b-1}t\varrho^t = \varrho\frac{d}{d\varrho}\sum_{t=0}^{b-1}\varrho^t = \varrho\,\dfrac{1 - b\varrho^{b-1} + (b-1)\varrho^b}{(1-\varrho)^2}$ です。通分すると

$$\frac{2b\varrho^b(1-\varrho) + 2\varrho - 2b\varrho^b + 2(b-1)\varrho^{b+1}}{b(1-\varrho)^2} = \frac{2\varrho - 2\varrho^{b+1}}{b(1-\varrho)^2}$$

です（$\varrho^b$ の項は $2b - 2b = 0$、$\varrho^{b+1}$ の項は $-2b + 2(b-1) = -2$）。上界は $2C_1/b = 2\varrho/(b(1-\varrho)^2)$ で、厳密な値はその $(1 - \varrho^b)$ 倍です。$b \to \infty$ で上界は等号に近づき、$O(1/b)$ という次数は改良できません。

演習 10.2易

$\hat\theta = \bar X_N$（$g(x) = x$）のとき、$v_{\mathrm J} = s_N^2/N$、$\hat\theta_{\mathrm J} = \bar X_N$ を示せ。

解答

定理 6.3 の証明の記号で $\hat\theta_{(i)} = \bar X_N + \Delta_i$、$\sum_i\Delta_i = 0$ なので $\bar\theta_{(\cdot)} = \bar X_N$ です。よって $\hat\theta_{\mathrm J} = N\bar X_N - (N-1)\bar X_N = \bar X_N$。分散は

$$v_{\mathrm J} = \frac{N-1}{N}\sum_i\Delta_i^2 = \frac{N-1}{N}\cdot\frac{\sum_i(X_i - \bar X_N)^2}{(N-1)^2} = \frac{s_N^2}{N}$$

で、独立な標本の平均の分散の不偏推定量に一致します。係数 $(N-1)/N$ はこの一致が成り立つように選ばれています。

演習 10.3標準

命題 7.5 の設定で、値がすべて異なるとき、再標本の最大値が上から $j+1$ 番目の値 $X_{(N-j)}$ に等しい確率を求め、$N \to \infty$ の極限を求めよ（$j \ge 0$ は固定）。

解答

$M_N^* \le X_{(N-j)}$ は、すべての $I_k$ が上位 $j$ 個以外の $N - j$ 個の添字のどれかであることと同値なので、確率は $(1 - j/N)^N$ です。したがって

$$\mathbb{P}^*(M_N^* = X_{(N-j)}) = \Bigl(1 - \frac jN\Bigr)^N - \Bigl(1 - \frac{j+1}{N}\Bigr)^N \to e^{-j} - e^{-(j+1)} = e^{-j}(1 - e^{-1})$$

です。再標本化分布は、上位の順序統計量の上に幾何分布の重み $(1 - e^{-1})e^{-j}$ を置いた離散分布に近づきます。$j = 0$ が命題の $1 - 1/e$ です。真の分布（連続）とは形が違い、どれだけ $N$ を増やしても近づきません。

演習 10.4標準

例 3.2 の 2 状態の連鎖（$0 < a \le 1/2$）を状態 $1$ から始め、$f = \varphi_2 = (1, -1)$ とする。熱化なしの時間平均 $N^{-1}\sum_{k=0}^{N-1}f(X_k)$ の偏りを厳密に求め、命題 8.2 の上界と比べよ。

解答

$\nu = \delta_1$ で $\mathbb{E}_\nu[f(X_k)] = (P^kf)(1) = \lambda_2^kf(1) = (1-2a)^k$、$\pi(f) = 0$ です。偏りは

$$\frac1N\sum_{k=0}^{N-1}(1-2a)^k = \frac{1 - (1-2a)^N}{2aN}$$

です。上界は、$\pi_{\min} = 1/2$ から $\sqrt{1/\pi_{\min} - 1} = 1$、$\mathrm{Var}_\pi(f) = 1$、$\gamma_* = 1 - |1 - 2a| = 2a$ なので $1/(2aN)$ です。厳密な値はその $(1 - (1-2a)^N)$ 倍で、$N \gg 1/a$ では上界がほぼ等号になります。統計誤差は $\sqrt{2\tau_{\mathrm{int}}/N} = \sqrt{(1-a)/(aN)}$ なので、偏りと統計誤差の比は $1/(2\sqrt{a(1-a)N})$ で、$N$ が緩和時間 $1/(2a)$ より十分大きければ偏りは無視できます。

参考文献

1. C. J. Geyer, “Practical Markov chain Monte Carlo”, Statistical Science 7 (1992), 473–483 — 可逆な連鎖の自己共分散の対の和の正値性・単調性・凸性と、初期正値列。
2. N. Madras, A. D. Sokal, “The pivot algorithm: a highly efficient Monte Carlo method for the self-avoiding walk”, Journal of Statistical Physics 50 (1988), 109–186 — 付録に窓つき推定量の分散の公式と窓の選び方。
3. U. Wolff, “Monte Carlo errors with less errors”, Computer Physics Communications 156 (2004), 143–153 — Γ 法と自動窓。arXiv:hep-lat/0306017
4. J. M. Flegal, G. L. Jones, “Batch means and spectral variance estimators in Markov chain Monte Carlo”, The Annals of Statistics 38 (2010), 1034–1070 — バッチ平均推定量の一致性の条件とバッチ長。
5. J. Shao, D. Tu, The Jackknife and Bootstrap, Springer, 1995 — ジャックナイフ分散推定量の一致性と中央値での不一致、ブートストラップの一致性と反例。
6. A. Gelman, D. B. Rubin, “Inference from iterative simulation using multiple sequences”, Statistical Science 7 (1992), 457–472 — 複数の連鎖による収束の診断。