# Markov 連鎖 Monte Carlo：提案の設計・Peskun の順序・混合時間 > 有限状態の可逆な Markov 連鎖について、時間平均の概収束、スペクトルギャップによる全変動距離と漸近分散の評価、採択関数の優劣を決める Peskun の順序、Cheeger の不等式、カップリングによる混合時間の上界を証明し、数値例で確かめる。 > https://rikai.mugen-giken.com/computer-science/stochastic-computation/markov-chain-monte-carlo ## 0. この記事の要点 複雑な確率分布から「くじ」を直接引くことができないとき、サイコロを振るように少しずつ状態を動かしていき、その軌跡の平均から答えを求める方法があります。この記事は、その方法がなぜ正しい値に近づくのか、どれくらい待てばよいのか、どんな動かし方が速いのかを、証明つきで調べます。読み終えると、動かし方の良し悪しを数や不等式で比べられるようになります。 - 目標の分布 $\pi$ から独立な標本が作れなくても、$\pi$ を[定常分布](/mathematics/probability/markov-chains-poisson#def-stationary)にもつ Markov 連鎖の**時間平均**は $\pi$ についての期待値に[概収束](/mathematics/probability/limit-theorems#def-convergence-ae)します（確率 $1$ で収束する、という意味です）。有限で既約な連鎖なら、周期的でも成り立ちます（）。 - 可逆な連鎖では、分布の収束の速さは[絶対スペクトルギャップ](/computer-science/stochastic-computation/statistical-error-analysis#def-absolute-gap) $\gamma^*$（遷移行列の固有値から決まる、収束の速さの目安）で決まります。$\|P^t(x,\cdot) - \pi\|_{\mathrm{TV}} \le \frac12\sqrt{(1-\pi(x))/\pi(x)}\,(1 - \gamma^*)^t$ で、混合時間は $1/\gamma^*$ の程度です（）。 - 定常な時間平均の漸近分散は $\sigma_f^2 = \sum_{j\ge2}\langle f, \varphi_j\rangle_\pi^2\frac{1+\lambda_j}{1-\lambda_j} \le (2/\gamma - 1)\mathrm{Var}_\pi(f)$ です。独立な標本に比べて、誤差の二乗は高々 $(2/\gamma - 1)$ 倍です（）。 - 同じ定常分布をもつ 2 つの可逆な連鎖で、**よく動くほう**（非対角成分が大きいほう）の漸近分散は、すべての $f$ について小さくなります（Peskun の順序。）。同じ提案に対する採択関数の中では Metropolis の $\min(1, r)$ が最良です（）。 - スペクトルギャップは、状態空間を 2 つに切ったときの流量の比（ボトルネック比 $\Phi_*$）で $\Phi_*^2/2 \le \gamma \le 2\Phi_*$ と挟まれます（）。2 つの山をもつ分布での局所的な更新は、山の間の確率に比例して遅くなります（）。 - 2 つの連鎖を同じ乱数で動かして合流させる**カップリング**は、固有値を求めずに混合時間の上界を与えます（）。 ## 1. 動機 [前章](/computer-science/stochastic-computation/monte-carlo-integration)の方法はどれも、目標の分布（または提案分布）から**独立な**標本を作れることを前提にしていました。ところが、統計力学の Boltzmann 分布 $\pi(x) \propto e^{-\beta H(x)}$ や Bayes 統計の事後分布では、状態 $x$ は数百から数百万の成分をもち、正規化定数は計算できず、[標本抽出の章](/computer-science/stochastic-computation/sampling-from-distributions)の逆関数法も棄却法も使えません（高次元では棄却法の採択率が指数的に小さくなります）。 Markov 連鎖 Monte Carlo 法（MCMC）は、独立性をあきらめます。$\pi$ を定常分布にもつ Markov 連鎖 $X_0, X_1, X_2, \ldots$ を作り、その時間平均 $\frac1N\sum_{k、、、）で一般の状態空間について扱われています。この章で使う有限状態の場合は、§2 で定義と証明を述べ直します。 ## 2. 準備：有限状態の可逆な連鎖 ### 2.1. 記号と仮定 $S$ を有限集合とします。遷移行列・Markov 連鎖・定常分布・可逆性（詳細つり合い）の定義は[Markov 連鎖と Poisson 過程](/mathematics/probability/markov-chains-poisson)（、）に従います。詳細つり合い $\pi(x)P(x,y) = \pi(y)P(y,x)$ を満たせば $\pi$ は定常分布です（）。 初期分布 $\mu$ の連鎖についての確率と期待値を $\mathbb{P}_\mu$、$\mathbb{E}_\mu$、1 点 $x$ から出発する場合を $\mathbb{P}_x$、$\mathbb{E}_x$ と書きます。Markov 連鎖の定義の条件付き確率を順に掛け合わせると、道の確率の公式 $$ \mathbb{P}_\mu(X_0 = x_0, X_1 = x_1, \ldots, X_k = x_k) = \mu(x_0)P(x_0, x_1)\cdots P(x_{k-1}, x_k) $$ が得られます（$k$ についての帰納法。途中の確率が $0$ なら両辺とも $0$ です）。$k$ 歩の遷移確率は行列の冪 $P^k(x, y)$ です（）。 遷移行列 $P$ が**既約**であるとは、任意の $x, y \in S$ に対して、ある $k \ge 1$ で $P^k(x, y) > 0$ となることをいう。既約な $P$ が**非周期的**であるとは、ある $m \ge 1$ で $P^m$ のすべての成分が正であることをいう。 既約で、ある状態に自己ループ $P(x_*, x_*) > 0$ があれば非周期的です（）。非周期性は、通常は「戻り時刻の最大公約数が $1$」と定義されます。有限で既約な連鎖ではこの定義と同値であることが知られていますが、この章では上の形を定義として使います。 有限集合 $S$ 上の既約な遷移行列 $P$ は定常分布 $\pi$ をもち、すべての $x$ で $\pi(x) > 0$ である。 $P_L = \frac12(I + P)$ と置きます。$P_L$ は遷移行列で、$P_L^k(x,y) \ge 2^{-k}P^k(x,y)$ から既約、$P_L(x,x) \ge \frac12 > 0$ です。 から $P_L$ は定常分布 $\pi$ をもちます。$\pi P_L = \pi$ は $\pi P = \pi$ と同値なので、$\pi$ は $P$ の定常分布です。$\pi(x_0) > 0$ となる $x_0$ を取ります。任意の $y$ に対して $P^k(x_0, y) > 0$ となる $k$ があり、$\pi = \pi P^k$ から $\pi(y) \ge \pi(x_0)P^k(x_0, y) > 0$ です。 定常分布の一意性は から従います。 ### 2.2. 採択関数の族と熱浴法 MCMC で使う連鎖の作り方を、有限状態の場合に述べます。 $\pi$ を $S$ 上の正の確率分布、$Q$ を $S$ 上の遷移行列（**提案**）とする。$x \ne y$、$Q(x,y) > 0$ となる組について $$ r(x, y) = \frac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)} \in [0, \infty) $$ と置く。関数 $a(x,y) \in [0,1]$（**採択関数**）が、$x \ne y$ なるすべての組で $$ \pi(x)Q(x,y)\,a(x,y) = \pi(y)Q(y,x)\,a(y,x) $$ を満たすとき、遷移行列 $$ P_a(x,y) = Q(x,y)\,a(x,y) \quad (x \ne y), \qquad P_a(x,x) = 1 - \sum_{y \ne x}Q(x,y)\,a(x,y) $$ を、提案 $Q$ と採択関数 $a$ による **Metropolis–Hastings 型の連鎖**という。「$x$ から $Q(x,\cdot)$ に従って $y$ を提案し、確率 $a(x,y)$ で $y$ に移り、そうでなければ $x$ にとどまる」連鎖である。次の 2 つが代表的である。 - **Metropolis–Hastings の採択関数**: $a_{\mathrm{MH}}(x,y) = \min(1, r(x,y))$。 - **Barker の採択関数**: $a_{\mathrm{B}}(x,y) = \dfrac{r(x,y)}{1 + r(x,y)} = \dfrac{\pi(y)Q(y,x)}{\pi(x)Q(x,y) + \pi(y)Q(y,x)}$。 （$Q(x,y) = 0$ の組では $a(x,y)$ の値は $P_a$ に影響しない。） の設定で次が成り立つ。 1. $P_a$ は $\pi$ について可逆である。したがって $\pi$ は $P_a$ の定常分布である。 2. $a_{\mathrm{MH}}$ と $a_{\mathrm{B}}$ は採択関数の条件を満たす。 3. 任意の採択関数 $a$ と、$x \ne y$、$Q(x,y) > 0$ なる組について $a(x,y) \le a_{\mathrm{MH}}(x,y)$。したがって $x \ne y$ で $P_a(x,y) \le P_{a_{\mathrm{MH}}}(x,y)$。 **1.** $x \ne y$ について $\pi(x)P_a(x,y) = \pi(x)Q(x,y)a(x,y) = \pi(y)Q(y,x)a(y,x) = \pi(y)P_a(y,x)$ です。$x = y$ では自明です。定常性は です。 **2.** $u = \pi(x)Q(x,y)$、$v = \pi(y)Q(y,x)$ と置きます。$u, v > 0$ のとき $r(x,y) = v/u$、$r(y,x) = u/v$ で、$u\min(1, v/u) = \min(u, v) = v\min(1, u/v)$、また $u\cdot\frac{v}{u+v} = v\cdot\frac{u}{u+v}$ です。$u > 0 = v$ のときは $r(x,y) = 0$ から $a(x,y) = 0$ で、両辺とも $0$ です（$u = 0$ の側の値は何でも構いません）。$u = v = 0$ のときも両辺 $0$ です。どちらの採択関数も $[0,1]$ に値を取ります。 **3.** 条件から $u\,a(x,y) = v\,a(y,x) \le v$ なので $a(x,y) \le v/u = r(x,y)$、また $a(x,y) \le 1$ です。よって $a(x,y) \le \min(1, r(x,y))$ です。 $Q$ が対称（$Q(x,y) = Q(y,x)$）なら $r(x,y) = \pi(y)/\pi(x)$ で、$\pi(x) \propto e^{-\beta H(x)}$ のとき $a_{\mathrm{MH}} = \min(1, e^{-\beta(H(y) - H(x))})$ です。これが本来の Metropolis 法です。どの採択関数も $\pi$ の比しか使わないので、正規化定数は不要です。3 は、採択関数の条件を満たす範囲で Metropolis–Hastings の採択確率が最大であることを述べています。それが推定の精度の意味でも最良であることが、§5 の主題です。 $S = S_1 \times \cdots \times S_n$（各 $S_i$ は有限集合）、$\pi$ を $S$ 上の正の確率分布とする。$x \in S$ の第 $i$ 成分だけを $s \in S_i$ に替えた状態を $x^{i \to s}$ と書く。遷移行列 $$ P_i(x, y) = \frac{\pi(y)}{\sum_{s \in S_i}\pi(x^{i\to s})} \quad (y \text{ は } x \text{ と第 } i \text{ 成分以外で一致}), \qquad P_i(x,y) = 0 \quad (\text{それ以外}) $$ を成分 $i$ の**熱浴更新**（Gibbs 更新）という。他の成分を固定したときの $\pi$ の条件付き分布から、第 $i$ 成分を引き直す操作である。$P_{\mathrm{RS}} = \frac1n\sum_{i=1}^nP_i$ を**無作為走査**、$P_{\mathrm{SS}} = P_1P_2\cdots P_n$ を**系統走査**の熱浴法という。 熱浴更新は、提案を「第 $i$ 成分を条件付き分布から引く」とした Metropolis–Hastings 法で、採択確率が常に $1$ になる場合にあたります。走査の順序による違いは §8 で調べます。 ### 2.3. スペクトル分解 $S$ 上の実数値関数 $f, g$ に対して、$\pi$ で重みを付けた内積と、平均・分散を $$ \langle f, g\rangle_\pi = \sum_x\pi(x)f(x)g(x), \qquad \pi(f) = \langle f, \mathbf{1}\rangle_\pi, \qquad \mathrm{Var}_\pi(f) = \langle f, f\rangle_\pi - \pi(f)^2 $$ と書きます。この内積を入れた関数の空間を $\ell^2(\pi)$ と書きます。ここでの狙いは、重み $\pi$ を内積に取り込むと、可逆な遷移行列が「対称行列」と同じ扱いになり、スペクトル定理（実対称行列は正規直交な固有ベクトルで対角化できる）が使える、という点です。遷移行列は関数に $(Pf)(x) = \sum_yP(x,y)f(y)$ で作用します。$P$ が $\pi$ について可逆であることは、$\langle f, Pg\rangle_\pi = \langle Pf, g\rangle_\pi$、すなわち $P$ が $\ell^2(\pi)$ の自己共役作用素であることと同値です（両辺を書き下すと $\sum_{x,y}\pi(x)P(x,y)f(x)g(y)$ と $\sum_{x,y}\pi(y)P(y,x)f(x)g(y)$ です）。 $P$ を、正の確率分布 $\pi$ について可逆な $S$ 上の遷移行列とし、$d = |S|$ とする。 1. 実数 $1 = \lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_d \ge -1$ と $S$ 上の実数値関数 $\varphi_1 = \mathbf{1}, \varphi_2, \ldots, \varphi_d$ で、$P\varphi_j = \lambda_j\varphi_j$、$\langle\varphi_i, \varphi_j\rangle_\pi = \delta_{ij}$ を満たすものが存在する。 2. 任意の $f$ について $f = \sum_j\langle f, \varphi_j\rangle_\pi\varphi_j$、$\langle f, f\rangle_\pi = \sum_j\langle f, \varphi_j\rangle_\pi^2$、$\mathrm{Var}_\pi(f) = \sum_{j \ge 2}\langle f, \varphi_j\rangle_\pi^2$ であり、すべての $t \ge 0$ と $x, y \in S$ について $$ P^t(x, y) = \pi(y)\sum_{j=1}^d\lambda_j^t\,\varphi_j(x)\varphi_j(y). $$ 特に $\sum_{j=1}^d\varphi_j(x)^2 = 1/\pi(x)$ である。 3. $P$ が既約ならば $\lambda_2 < 1$。さらに非周期的ならば $\lambda_d > -1$。 **1.** $D = \mathrm{diag}(\pi)$ とし、$A = D^{1/2}PD^{-1/2}$、すなわち $A(x,y) = \pi(x)P(x,y)/\sqrt{\pi(x)\pi(y)}$ と置きます。可逆性から $A$ は実対称行列なので、実の固有値 $\lambda_1 \ge \cdots \ge \lambda_d$ と、$\mathbb{R}^d$ の正規直交基底をなす固有ベクトル $v_1, \ldots, v_d$ があります（）。$\varphi_j = D^{-1/2}v_j$ と置くと、$P\varphi_j = D^{-1/2}Av_j = \lambda_j\varphi_j$、$\langle\varphi_i, \varphi_j\rangle_\pi = \sum_x v_i(x)v_j(x) = \delta_{ij}$ です。 $|\lambda_j| \le 1$ を示します。$|\varphi_j(x)|$ が最大になる $x^*$ を取ると、$|\lambda_j||\varphi_j(x^*)| = |\sum_yP(x^*,y)\varphi_j(y)| \le \sum_yP(x^*,y)|\varphi_j(x^*)| = |\varphi_j(x^*)|$ で、$\varphi_j(x^*) \ne 0$ なので $|\lambda_j| \le 1$ です。$P\mathbf{1} = \mathbf{1}$ なので $1$ は固有値で、最大の固有値は $\lambda_1 = 1$ です。$\langle\mathbf{1}, \mathbf{1}\rangle_\pi = 1$ なので、固有値 $1$ の固有空間の正規直交基底を $\mathbf{1}$ を含むように取り直して $\varphi_1 = \mathbf{1}$ とできます。 **2.** $(\varphi_j)$ は $\ell^2(\pi)$ の正規直交基底なので、展開と $\langle f, f\rangle_\pi = \sum_j\langle f, \varphi_j\rangle_\pi^2$ が成り立ちます。$\langle f, \varphi_1\rangle_\pi = \pi(f)$ を引いて分散の式を得ます。$P^tf = \sum_j\lambda_j^t\langle f, \varphi_j\rangle_\pi\varphi_j$ を、$y$ の指示関数 $f = \mathbf{1}_{\{y\}}$ に使うと、$\langle\mathbf{1}_{\{y\}}, \varphi_j\rangle_\pi = \pi(y)\varphi_j(y)$、$(P^t\mathbf{1}_{\{y\}})(x) = P^t(x,y)$ から主張の式になります。$t = 0$、$y = x$ とすると $1 = \pi(x)\sum_j\varphi_j(x)^2$ です。 **3.** $Pf = f$ とします。Dirichlet 形式の表示（）から $$ 0 = \langle f, (I - P)f\rangle_\pi = \frac12\sum_{x,y}\pi(x)P(x,y)\bigl(f(x) - f(y)\bigr)^2 $$ なので、$P(x,y) > 0$ なら $f(x) = f(y)$ です。既約性から任意の 2 点は正の確率の道でつながるので、$f$ は定数です。固有値 $1$ の固有空間は 1 次元で、$\lambda_2 < 1$ です。 非周期的とし、$Pf = -f$ とします。$P^2f = f$ です。$P^2$ は $\pi$ について可逆な遷移行列で、ある $m$ で $(P^2)^m = P^mP^m$ の全成分が正なので既約です。前半を $P^2$ に使うと $f$ は定数で、$Pf = f$ です。$f = -f$ から $f = 0$ となり、$-1$ は固有値ではありません。 の記号で、 $$ \gamma = 1 - \lambda_2, \qquad \lambda_* = \max\{|\lambda_2|, |\lambda_d|\}, \qquad \gamma^* = 1 - \lambda_* $$ を、それぞれ**スペクトルギャップ**、第 2 固有値の絶対値、**絶対スペクトルギャップ**という（$d = 1$ のときは $\lambda_* = 0$ とする）。$\gamma^* > 0$ のとき $t_{\mathrm{rel}} = 1/\gamma^*$ を**緩和時間**という。遷移行列 $P_L = \frac12(I + P)$ を $P$ の**怠惰な連鎖**という。 $\gamma$ は Dirichlet 形式による定義 $\lambda(P) = \min\mathcal{E}_P(f,f)/\mathrm{Var}_\pi(f)$ と一致します（ の 1）。怠惰な連鎖の固有値は $\frac12(1 + \lambda_j) \ge 0$ なので、負の固有値がなくなり、$\gamma^*(P_L) = \gamma(P_L) = \gamma(P)/2$ です。周期性や負の固有値を避けるために、理論ではしばしば怠惰な連鎖を使います。 ## 3. 時間平均の収束 MCMC の推定量は時間平均です。その収束は、分布の収束（§4）とは別の、より弱い仮定で成り立ちます。 次の補題と定理の方針は、連鎖が状態 $o$ に戻るたびに連鎖を区切る、というものです。区切られた各区間は互いに独立で同じ分布に従うので、独立な標本に対する大数の法則が使えます。補題は、$o$ に必ず戻ってくること（しかも戻るまでの時間の平均が有限であること）を保証します。 $P$ を有限集合 $S$ 上の既約な遷移行列、$o \in S$ とし、$T_o = \min\{k \ge 1 : X_k = o\}$（そのような $k$ がなければ $\infty$）と置く。このとき整数 $m \ge 1$ と $\varepsilon > 0$ があって、すべての $y \in S$ と $k \ge 0$ について $$ \mathbb{P}_y(T_o > km) \le (1 - \varepsilon)^k $$ である。特に、任意の初期分布について確率 $1$ で $T_o < \infty$ であり、$\mathbb{E}_y[T_o] \le m/\varepsilon$ である。 各 $y$ について、$P^{k_y}(y, o) > 0$ となる $k_y \ge 1$ を選び、$m = \max_yk_y$、$\varepsilon = \min_yP^{k_y}(y, o) > 0$ と置きます。$\mathbb{P}_y(T_o \le m) \ge \mathbb{P}_y(X_{k_y} = o) = P^{k_y}(y, o) \ge \varepsilon$ なので、すべての $y$ で $\mathbb{P}_y(T_o > m) \le 1 - \varepsilon$ です。 道の確率の公式から、$\mathbb{P}_y(T_o > (k+1)m)$ は、$y$ から出発して $(k+1)m$ 歩の間 $o$ を踏まない道の確率の和です。最初の $km$ 歩（$o$ を踏まず、$z$ で終わる）と残りの $m$ 歩（$z$ から出発して $o$ を踏まない）に分けると、 $$ \mathbb{P}_y(T_o > (k+1)m) = \sum_z\mathbb{P}_y(T_o > km,\ X_{km} = z)\,\mathbb{P}_z(T_o > m) \le (1 - \varepsilon)\,\mathbb{P}_y(T_o > km) $$ です。帰納法で主張を得ます。$k \to \infty$ として $\mathbb{P}_y(T_o = \infty) = 0$、また $\mathbb{E}_y[T_o] = \sum_{j \ge 0}\mathbb{P}_y(T_o > j) \le m\sum_{k \ge 0}(1 - \varepsilon)^k = m/\varepsilon$ です（$\mathbb{P}_y(T_o > j)$ は $j$ について減少するので、$km \le j < (k+1)m$ の $m$ 個の項をそれぞれ $(1-\varepsilon)^k$ で抑えました）。 $P$ を有限集合 $S$ 上の既約な遷移行列、$\pi$ をその定常分布とする。任意の初期分布 $\mu$ と任意の $f : S \to \mathbb{R}$ について、 $$ \mathbb{P}_\mu\Bigl(\lim_{N\to\infty}\frac1N\sum_{k=0}^{N-1}f(X_k) = \pi(f)\Bigr) = 1 $$ である。特に、定常分布は一意である。非周期性は仮定しない。 状態 $o \in S$ を 1 つ固定します。 **区間への分割。** $o$ への到達時刻を $\tau_0 = \min\{k \ge 0 : X_k = o\}$、$\tau_{j+1} = \min\{k > \tau_j : X_k = o\}$ と定めます。$o$ から始まり、途中で $o$ を踏まない有限の語 $w = (o, z_1, \ldots, z_{l-1})$（$l \ge 1$、$z_i \ne o$）の全体を $\mathcal{W}$ とし、 $$ p(w) = P(o, z_1)P(z_1, z_2)\cdots P(z_{l-1}, o) $$ と置きます（$l = 1$ のときは $p(w) = P(o,o)$）。$p(w)$ は、$o$ から出発した連鎖が $w$ をたどって時刻 $l$ にはじめて $o$ に戻る確率なので、 から $\sum_{w \in \mathcal{W}}p(w) = \mathbb{P}_o(T_o < \infty) = 1$ です。連鎖の第 $j$ 区間を $W_j = (X_{\tau_j}, \ldots, X_{\tau_{j+1}-1})$ とします。 $o$ を含まない語 $a = (a_0, \ldots, a_{k-1})$（$k \ge 0$）と $w_0, \ldots, w_j \in \mathcal{W}$ に対して、事象「$(X_0, \ldots, X_{\tau_0 - 1}) = a$ かつ $W_0 = w_0, \ldots, W_j = w_j$」は、連鎖の最初の部分が $a, w_0, \ldots, w_j$ をつなげた語で、その次が $o$ であることと同じです。道の確率の公式から、その確率は $$ q(a)\,p(w_0)\cdots p(w_j), \qquad q(a) = \mu(a_0)P(a_0, a_1)\cdots P(a_{k-1}, o) $$ です（$k = 0$ のときは $q(a) = \mu(o)$）。$\sum_aq(a) = \mathbb{P}_\mu(\tau_0 < \infty) = 1$ なので（。$X_0 = o$ なら $\tau_0 = 0$、そうでなければ $\tau_0 = T_o$）、$a$ について和を取ると $$ \mathbb{P}_\mu(W_0 = w_0, \ldots, W_j = w_j) = p(w_0)\cdots p(w_j) $$ です。すなわち、区間 $W_0, W_1, \ldots$ は独立で、どれも分布 $p$ に従います。また確率 $1$ ですべての $\tau_j$ は有限です。 **大数の法則。** 語 $w$ の長さを $\ell(w)$、$w$ に沿った $f$ の和を $B(w)$ と書き、$\ell_j = \ell(W_j) = \tau_{j+1} - \tau_j$、$B_j = B(W_j)$ と置きます。$(\ell_j)$ と $(B_j)$ はそれぞれ独立同分布の列で、$\mathbb{E}[\ell_0] = \mathbb{E}_o[T_o] < \infty$（）、$|B_j| \le \|f\|_\infty\ell_j$（$\|f\|_\infty = \max|f|$）から $B_j$ も可積分です。大数の強法則（）から、確率 $1$ で $$ \frac1J\sum_{j。優関数は定数 $\|f\|_\infty$）から期待値も $c(f)$ に収束し、$c(f) = \pi(f)$ です。 定常分布が $\pi, \pi'$ の 2 つあれば、どちらについても $\pi(f) = c(f) = \pi'(f)$ がすべての $f$ で成り立つので、$\pi = \pi'$ です。 この定理は、MCMC の推定量が**一致性をもつ**ことを保証します。初期状態は何でもよく、連鎖が周期的でも構いません（周期的な連鎖では $X_N$ の分布は収束しませんが、時間平均は収束します）。必要なのは、既約性と、$\pi$ が定常分布であることだけです。ただし、定理は収束の**速さ**について何も言いません。有限の $N$ での誤差を評価するには、次の 2 つの節の定量的な結果が要ります。 ## 4. 全変動距離と混合時間 $S$ 上の確率分布 $\mu, \nu$ の**全変動距離**を $\|\mu - \nu\|_{\mathrm{TV}} = \frac12\sum_x|\mu(x) - \nu(x)|$ と定める。定常分布 $\pi$ をもつ遷移行列 $P$ について $$ d(t) = \max_{x \in S}\|P^t(x, \cdot) - \pi\|_{\mathrm{TV}}, \qquad t_{\mathrm{mix}}(\varepsilon) = \min\{t \ge 0 : d(t) \le \varepsilon\} $$ と置き、$t_{\mathrm{mix}}(\varepsilon)$ を**混合時間**という（そのような $t$ がなければ $\infty$）。 1. $\|\mu - \nu\|_{\mathrm{TV}} = \max_{A \subseteq S}\bigl(\mu(A) - \nu(A)\bigr)$。 2. 遷移行列 $P$ について $\|\mu P - \nu P\|_{\mathrm{TV}} \le \|\mu - \nu\|_{\mathrm{TV}}$。したがって $d(t)$ は $t$ について単調非増加である。 3. $d(t) \le \max_{x,y}\|P^t(x,\cdot) - P^t(y,\cdot)\|_{\mathrm{TV}}$。 **1.** $B = \{x : \mu(x) > \nu(x)\}$ と置きます。$\sum_x(\mu(x) - \nu(x)) = 0$ から $\sum_{x \in B}(\mu - \nu)(x) = \sum_{x \notin B}(\nu - \mu)(x)$ で、両者の和が $\sum_x|\mu(x) - \nu(x)|$ なので、$\mu(B) - \nu(B) = \|\mu - \nu\|_{\mathrm{TV}}$ です。任意の $A$ について $\mu(A) - \nu(A) = \sum_{x\in A}(\mu - \nu)(x) \le \sum_{x \in A \cap B}(\mu-\nu)(x) \le \mu(B) - \nu(B)$ です。 **2.** $\sum_y|\sum_x(\mu - \nu)(x)P(x,y)| \le \sum_x|\mu - \nu|(x)\sum_yP(x,y) = \sum_x|\mu(x) - \nu(x)|$ です。$\mu = P^t(x,\cdot)$、$\nu = \pi = \pi P$ とすると $\|P^{t+1}(x,\cdot) - \pi\|_{\mathrm{TV}} \le \|P^t(x,\cdot) - \pi\|_{\mathrm{TV}}$ です。 **3.** $\pi = \pi P^t$ から $P^t(x, \cdot) - \pi = \sum_y\pi(y)\bigl(P^t(x,\cdot) - P^t(y,\cdot)\bigr)$ です。三角不等式と $\sum_y\pi(y) = 1$ から従います。 $P$ を、正の確率分布 $\pi$ について可逆で既約な $S$ 上の遷移行列とし、$\pi_{\min} = \min_x\pi(x)$ と置く。 の記号を使う。 1. すべての $x \in S$ と $t \ge 0$ について $$ \|P^t(x, \cdot) - \pi\|_{\mathrm{TV}} \le \frac12\sqrt{\frac{1 - \pi(x)}{\pi(x)}}\,(1 - \gamma^*)^t. $$ 2. $\gamma^* > 0$ ならば、$0 < \varepsilon < 1$ について $$ t_{\mathrm{mix}}(\varepsilon) \le \max\Bigl\{0,\ \Bigl\lceil\frac{1}{\gamma^*}\ln\frac{1}{2\varepsilon\sqrt{\pi_{\min}}}\Bigr\rceil\Bigr\}. $$ 3. すべての $t \ge 0$ で $\lambda_*^t \le 2d(t)$ である。したがって $0 < \gamma^* < 1$、$0 < \varepsilon < \frac12$ ならば $$ t_{\mathrm{mix}}(\varepsilon) \ge \Bigl(\frac{1}{\gamma^*} - 1\Bigr)\ln\frac{1}{2\varepsilon}. $$ の $\lambda_j, \varphi_j$ を使います。 **1.** $x$ を固定し、$h(y) = P^t(x,y)/\pi(y) - 1$ と置きます。 の 2 と $\lambda_1 = 1$、$\varphi_1 = \mathbf{1}$ から $h = \sum_{j \ge 2}\lambda_j^t\varphi_j(x)\,\varphi_j$ で、正規直交性から $\langle h, h\rangle_\pi = \sum_{j\ge2}\lambda_j^{2t}\varphi_j(x)^2$ です。Cauchy–Schwarz の不等式（$\ell^2(\pi)$ での $|h|$ と $\mathbf{1}$ の内積）から $$ 4\|P^t(x,\cdot) - \pi\|_{\mathrm{TV}}^2 = \Bigl(\sum_y\pi(y)|h(y)|\Bigr)^2 \le \langle h, h\rangle_\pi = \sum_{j \ge 2}\lambda_j^{2t}\varphi_j(x)^2 \le \lambda_*^{2t}\sum_{j \ge 2}\varphi_j(x)^2 $$ です。$\sum_{j\ge1}\varphi_j(x)^2 = 1/\pi(x)$ と $\varphi_1(x)^2 = 1$ から $\sum_{j\ge2}\varphi_j(x)^2 = (1 - \pi(x))/\pi(x)$ です。平方根を取り、$\lambda_* = 1 - \gamma^*$ を代入します。 **2.** $1 - \gamma^* \le e^{-\gamma^*}$ と 1 から $d(t) \le \frac{1}{2\sqrt{\pi_{\min}}}e^{-\gamma^*t}$ です。右辺が $\varepsilon$ 以下になるのは $t \ge \frac{1}{\gamma^*}\ln\frac{1}{2\varepsilon\sqrt{\pi_{\min}}}$ のときです。$t$ は $0$ 以上の整数なので、この値が負のとき（$2\varepsilon\sqrt{\pi_{\min}} > 1$。状態数が 3 以下の場合にしか起きません）は $t = 0$ で既に成り立ち、$t_{\mathrm{mix}}(\varepsilon) = 0$ です。 **3.** $d = 1$ なら $\lambda_* = 0$ で自明なので $d \ge 2$ とし、$|\lambda_j| = \lambda_*$ となる $j \ge 2$ を取ります。$\varphi_j$ は $\varphi_1 = \mathbf{1}$ と直交するので $\sum_y\pi(y)\varphi_j(y) = 0$ です。よって任意の $x$ で $$ \lambda_j^t\varphi_j(x) = (P^t\varphi_j)(x) = \sum_y\bigl(P^t(x,y) - \pi(y)\bigr)\varphi_j(y) $$ で、絶対値を取ると $\lambda_*^t|\varphi_j(x)| \le \|\varphi_j\|_\infty\sum_y|P^t(x,y) - \pi(y)| \le 2d(t)\|\varphi_j\|_\infty$ です。$|\varphi_j(x)| = \|\varphi_j\|_\infty > 0$ となる $x$ を取って $\lambda_*^t \le 2d(t)$ です。 $t = t_{\mathrm{mix}}(\varepsilon)$（有限としてよい）では $\lambda_*^t \le 2\varepsilon$ なので、$t\ln(1/\lambda_*) \ge \ln(1/(2\varepsilon))$ です。$\ln u \le u - 1$ を $u = 1/\lambda_*$ に使うと $\ln(1/\lambda_*) \le \frac{1}{\lambda_*} - 1 = \frac{\gamma^*}{1 - \gamma^*}$ なので、$t \ge \frac{1-\gamma^*}{\gamma^*}\ln\frac{1}{2\varepsilon}$ です。 1 と 3 を合わせると、混合時間は緩和時間 $t_{\mathrm{rel}} = 1/\gamma^*$ で上下から挟まれます。上界には $\ln(1/\pi_{\min})$ の因子が付きます。これは、最悪の初期状態（確率の小さい 1 点）から出発することの代価で、状態数が指数的に大きい系では無視できない大きさです（$\{0,1\}^n$ なら $\ln(1/\pi_{\min}) = n\ln2$）。 $S = \{1, 2\}$、$P = \begin{pmatrix}1-a & a\\ b & 1-b\end{pmatrix}$（$0 < a, b \le 1$）は $\pi = (\frac{b}{a+b}, \frac{a}{a+b})$ について可逆で、固有値は $1$ と $\lambda_2 = 1 - a - b$ です。$\varphi_2 = (\sqrt{a/b}, -\sqrt{b/a})$ が $\ell^2(\pi)$ で正規化された固有関数で（$\pi(\varphi_2) = 0$、$\pi(\varphi_2^2) = \frac{b}{a+b}\cdot\frac ab + \frac{a}{a+b}\cdot\frac ba = 1$）、 の 2 から $$ P^t(1, 2) - \pi(2) = \pi(2)\lambda_2^t\varphi_2(1)\varphi_2(2) = -\frac{a}{a+b}\lambda_2^t, \qquad \|P^t(1,\cdot) - \pi\|_{\mathrm{TV}} = \frac{a}{a+b}|\lambda_2|^t $$ です。 の 1 の右辺は $\frac12\sqrt{a/b}\,|\lambda_2|^t$ です。相加・相乗平均の不等式から $\frac{a}{a+b} \le \frac{a}{2\sqrt{ab}} = \frac12\sqrt{a/b}$ で、$a = b$ のとき等号が成り立ちます。対称な 2 状態の連鎖では、定理の評価はこれ以上改善できません。$a + b = 1$ のときは $\lambda_2 = 0$ で、1 歩で定常分布に達します（独立な標本です）。$a = b = 1$ のときは $\lambda_2 = -1$、$\gamma^* = 0$ で、連鎖は周期 2 で分布は収束しません。 - **完全グラフ。** $P(x,y) = 1/d$（すべての $x, y$）では $P^2 = P$ で、固有値は $1$ と $0$ だけです。$\gamma = \gamma^* = 1$ で、1 歩で混合します。 - **円周。** $\mathbb{Z}/n$ 上の単純ランダムウォーク（$n$ は奇数）では $\gamma = 1 - \cos(2\pi/n) \approx 2\pi^2/n^2$ です（）。拡散的な動きでは、端から端までの距離の 2 乗に比例する時間がかかります。 - **超立方体。** $S = \{0,1\}^n$ の上で、「座標 $i$ を一様に選び、その成分を一様な乱数ビットで置き換える」連鎖（$\frac12$ の確率でとどまる怠惰なランダムウォーク）を考えます。$\pi$ は一様分布です。$u \in \{0,1\}^n$ に対して $\chi_u(x) = (-1)^{\langle u, x\rangle}$ と置くと、選ばれた座標 $i$ が $u_i = 0$ なら $\chi_u$ の値は変わらず、$u_i = 1$ なら新しいビットについての平均が $0$ なので、 $$ (P\chi_u)(x) = \Bigl(1 - \frac{|u|}{n}\Bigr)\chi_u(x), \qquad |u| = u_1 + \cdots + u_n $$ です。$2^n$ 個の $\chi_u$ は $\ell^2(\pi)$ の正規直交基底なので、固有値は $1 - k/n$（$k = 0, \ldots, n$、重複度 $\binom nk$）で尽きます。$\gamma = \gamma^* = 1/n$ です。 の 2 は $t_{\mathrm{mix}}(\varepsilon) \le \lceil n(\frac n2\ln2 + \ln\frac{1}{2\varepsilon})\rceil$ を与えますが、これは $n^2$ の程度で、真の大きさ $n\ln n$ の程度（§7）より粗い評価です。 ## 5. 漸近分散と Peskun の順序 ### 5.1. 定常な時間平均の分散 $\pi$ を定常分布にもつ連鎖を $X_0 \sim \pi$ から出発させる。$f : S \to \mathbb{R}$ に対して、極限 $$ \sigma_f^2 = \lim_{N\to\infty}N\,\mathrm{Var}_\pi\Bigl(\frac1N\sum_{k=0}^{N-1}f(X_k)\Bigr) $$ が存在するとき、これを $f$ の**漸近分散**という（$\mathrm{Var}_\pi$ の添字は、ここでは初期分布が $\pi$ であることを表す）。 独立な標本なら $\sigma_f^2 = \mathrm{Var}_\pi(f)$ です（ の 2）。比 $\sigma_f^2/\mathrm{Var}_\pi(f)$ は、「連鎖の何歩が独立な標本 1 個に相当するか」を表します。 $P$ を、正の確率分布 $\pi$ について可逆で既約な $S$ 上の遷移行列とする。任意の $f : S \to \mathbb{R}$ について漸近分散が存在し、 の記号で $$ \sigma_f^2 = \sum_{j=2}^{d}\langle f, \varphi_j\rangle_\pi^2\,\frac{1 + \lambda_j}{1 - \lambda_j} \le \Bigl(\frac{2}{\gamma} - 1\Bigr)\mathrm{Var}_\pi(f) $$ である。$f$ が $\lambda_2$ の固有関数のとき等号が成り立つ。 $a_j = \langle f, \varphi_j\rangle_\pi$ と置きます。$X_0 \sim \pi$ のとき、$\mathbb{E}_\pi[f(X_0)f(X_k)] = \sum_x\pi(x)f(x)(P^kf)(x) = \langle f, P^kf\rangle_\pi = \sum_j\lambda_j^ka_j^2$（ の 2）、$\mathbb{E}_\pi[f(X_k)] = \pi(f) = a_1$ なので、自己共分散は $$ C(k) := \mathrm{Cov}_\pi\bigl(f(X_0), f(X_k)\bigr) = \sum_{j \ge 2}a_j^2\lambda_j^k $$ です。定常性から $\mathrm{Cov}_\pi(f(X_i), f(X_{i+k})) = C(k)$ で、分散の和の公式（ の 4）から $$ N\,\mathrm{Var}_\pi\Bigl(\frac1N\sum_{k の 3）。$\lambda \in [-1, 1)$ について、部分和 $s_m = \sum_{k=1}^m\lambda^k = \lambda(1 - \lambda^m)/(1 - \lambda)$ を使うと $$ \sum_{k=1}^{N-1}\Bigl(1 - \frac kN\Bigr)\lambda^k = \frac1N\sum_{m=1}^{N-1}s_m = \frac{\lambda}{1-\lambda}\Bigl(\frac{N-1}{N} - \frac1N\sum_{m=1}^{N-1}\lambda^m\Bigr) $$ です（最初の等号は、$\lambda^k$ が $s_k, s_{k+1}, \ldots, s_{N-1}$ の $N - k$ 個に現れることから）。$|\lambda| < 1$ なら $\sum_m\lambda^m$ は有界、$\lambda = -1$ なら $\sum_{m=1}^{N-1}(-1)^m \in \{-1, 0\}$ なので、どちらの場合も $\frac1N\sum_{m 式の形から 3 つのことが読み取れます。第一に、$\lambda_j$ が $1$ に近い成分（ゆっくり緩和するモード）が分散を支配します。第二に、**負の固有値は分散を減らします**。$\lambda_j < 0$ の成分の係数は $1$ より小さく、独立な標本より精度がよくなります（隣り合う標本が負に相関する、対称変量と同じ効果です）。分布の収束（）では $|\lambda_j|$ が効くので負の固有値は害ですが、時間平均の精度では益です。第三に、上界に現れるのは絶対ギャップ $\gamma^*$ ではなく $\gamma$ です。2 状態の連鎖（）では $j = 2$ しかないので、すべての $f$ で $\sigma_f^2 = (\frac2\gamma - 1)\mathrm{Var}_\pi(f)$ と等号が成り立ちます。 ### 5.2. Peskun の順序 同じ $\pi$ を定常分布にもつ連鎖は無数にあります。どれを使うべきかを決める一般的な原理が、次の定理です。準備として、[正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite)な作用素の逆についての補題を示します。漸近分散は、§5.1 の式の書き換えで $(I - P)$ の逆を含む形になります。そのため、「$I - P$ が大きいほど逆は小さい」という大小関係が、連鎖の比較に必要になります。 $V$ を有限次元の実内積空間、$A, B$ を $V$ 上の自己共役な線形写像とし、すべての $g \ne 0$ で $\langle g, Ag\rangle \ge \langle g, Bg\rangle > 0$ とする。このとき $A, B$ は可逆で、すべての $f \in V$ について $$ \langle f, A^{-1}f\rangle \le \langle f, B^{-1}f\rangle. $$ $Bg = 0$ なら $\langle g, Bg\rangle = 0$ から $g = 0$ なので、$B$ は単射で、有限次元なので可逆です。$A$ も同様です。$B$ が自己共役で $B^{-1}$ もそうであることを使うと、任意の $f, g$ について $$ \bigl\langle g - B^{-1}f,\ B(g - B^{-1}f)\bigr\rangle = \langle g, Bg\rangle - 2\langle f, g\rangle + \langle f, B^{-1}f\rangle \ge 0 $$ で、等号は $g = B^{-1}f$ のときです。すなわち $$ \langle f, B^{-1}f\rangle = \max_{g \in V}\bigl(2\langle f, g\rangle - \langle g, Bg\rangle\bigr) $$ です。$A$ についても同じ式が成り立ちます。仮定から、各 $g$ で $2\langle f, g\rangle - \langle g, Ag\rangle \le 2\langle f, g\rangle - \langle g, Bg\rangle$ なので、最大値を比べて結論を得ます。 $P_1, P_2$ を、同じ正の確率分布 $\pi$ について可逆で既約な $S$ 上の遷移行列とし、すべての $x \ne y$ で $$ P_1(x, y) \ge P_2(x, y) $$ とする。このとき次が成り立つ。 1. すべての $g$ で $\langle g, (I - P_1)g\rangle_\pi \ge \langle g, (I - P_2)g\rangle_\pi$。したがって $\gamma(P_1) \ge \gamma(P_2)$。 2. すべての $f : S \to \mathbb{R}$ で $\sigma_f^2(P_1) \le \sigma_f^2(P_2)$。 **1.** Dirichlet 形式の表示（） $$ \langle g, (I - P_i)g\rangle_\pi = \frac12\sum_{x \ne y}\pi(x)P_i(x,y)\bigl(g(x) - g(y)\bigr)^2 $$ には非対角成分しか現れず（$x = y$ の項は $0$）、各項の係数は非負です。$P_1(x,y) \ge P_2(x,y)$ から不等式が従います。スペクトルギャップはこの量と $\mathrm{Var}_\pi(g)$ の比の最小値なので（ の 1）、$\gamma(P_1) \ge \gamma(P_2)$ です。 **2.** $H_0 = \{g : \pi(g) = 0\}$ を、内積 $\langle\cdot,\cdot\rangle_\pi$ をもつ実内積空間とみます。$P = P_1$ または $P_2$ について、定常性から $\pi(Pg) = \pi(g)$ なので $I - P$ は $H_0$ を $H_0$ に写し、可逆性から自己共役です。この制限を $A_P$ と書きます。 の基底では、$H_0$ は $\varphi_2, \ldots, \varphi_d$ で張られ、$A_P\varphi_j = (1 - \lambda_j)\varphi_j$ です。既約性から $1 - \lambda_j \ge \gamma > 0$ なので、$g \in H_0$、$g \ne 0$ について $\langle g, A_Pg\rangle_\pi = \sum_{j\ge2}(1 - \lambda_j)\langle g, \varphi_j\rangle_\pi^2 > 0$ で、$A_P^{-1}\varphi_j = \varphi_j/(1 - \lambda_j)$ です。 $f_0 = f - \pi(f) \in H_0$ と置くと、$j \ge 2$ で $\langle f_0, \varphi_j\rangle_\pi = \langle f, \varphi_j\rangle_\pi =: a_j$ です。 の式で $\frac{1+\lambda_j}{1-\lambda_j} = \frac{2}{1-\lambda_j} - 1$ と書き直すと $$ \sigma_f^2(P) = 2\sum_{j \ge 2}\frac{a_j^2}{1 - \lambda_j} - \sum_{j \ge 2}a_j^2 = 2\,\langle f_0, A_P^{-1}f_0\rangle_\pi - \langle f_0, f_0\rangle_\pi $$ です。1 から、すべての $g \in H_0$、$g \ne 0$ で $\langle g, A_{P_1}g\rangle_\pi \ge \langle g, A_{P_2}g\rangle_\pi > 0$ なので、 から $\langle f_0, A_{P_1}^{-1}f_0\rangle_\pi \le \langle f_0, A_{P_2}^{-1}f_0\rangle_\pi$ です。第 2 項 $\langle f_0, f_0\rangle_\pi = \mathrm{Var}_\pi(f)$ は連鎖によらないので、$\sigma_f^2(P_1) \le \sigma_f^2(P_2)$ です。 「現在の状態にとどまる確率が小さい連鎖ほど、**どんな量を測っても**精度がよい」という、強い主張です。仮定は成分ごとの不等式なので、検証は簡単です。ただし、比較できるのは同じ移動先への確率が一斉に大きい場合だけで、移動先の集合が違う連鎖（局所的な提案と大域的な提案）は、この順序では比べられません。 $\pi$ を正の確率分布、$Q$ を提案とし、$a$ を採択関数（）で $P_a$ が既約になるものとする。 1. $P_{a_{\mathrm{MH}}}$ は既約で、すべての $f$ について $\sigma_f^2(P_{a_{\mathrm{MH}}}) \le \sigma_f^2(P_a)$、また $\gamma(P_{a_{\mathrm{MH}}}) \ge \gamma(P_a)$。 2. Barker の採択関数による連鎖 $P_{a_{\mathrm{B}}}$ が既約ならば、すべての $f$ で $$ \sigma_f^2(P_{a_{\mathrm{MH}}}) \le \sigma_f^2(P_{a_{\mathrm{B}}}) \le 2\,\sigma_f^2(P_{a_{\mathrm{MH}}}) + \mathrm{Var}_\pi(f). $$ **1.** の 3 から、$x \ne y$ で $P_{a_{\mathrm{MH}}}(x,y) \ge P_a(x,y)$ です。$P_a$ で正の確率をもつ道は $P_{a_{\mathrm{MH}}}$ でも正の確率をもつので、$P_{a_{\mathrm{MH}}}$ は既約です。どちらも $\pi$ について可逆なので（同じ命題の 1）、 が使えます。 **2.** 左の不等式は 1 を $a = a_{\mathrm{B}}$ に使ったものです（$a_{\mathrm{B}}$ が採択関数であることは の 2）。右の不等式を示します。$r \ge 0$ について、$r \le 1$ なら $\frac{r}{1+r} \ge \frac r2$、$r \ge 1$ なら $\frac{r}{1+r} \ge \frac12$ なので、$a_{\mathrm{B}} \ge \frac12a_{\mathrm{MH}}$、したがって $x \ne y$ で $P_{a_{\mathrm{B}}}(x,y) \ge \frac12P_{a_{\mathrm{MH}}}(x,y)$ です。 の証明の記号で、Dirichlet 形式の表示から $\langle g, A_{\mathrm{B}}g\rangle_\pi \ge \langle g, \frac12A_{\mathrm{MH}}g\rangle_\pi > 0$（$g \in H_0$、$g \ne 0$）です。 を $A = A_{\mathrm{B}}$、$B = \frac12A_{\mathrm{MH}}$ に使うと $\langle f_0, A_{\mathrm{B}}^{-1}f_0\rangle_\pi \le 2\langle f_0, A_{\mathrm{MH}}^{-1}f_0\rangle_\pi$ です。よって $$ \sigma_f^2(P_{a_{\mathrm{B}}}) = 2\langle f_0, A_{\mathrm{B}}^{-1}f_0\rangle_\pi - \mathrm{Var}_\pi(f) \le 4\langle f_0, A_{\mathrm{MH}}^{-1}f_0\rangle_\pi - \mathrm{Var}_\pi(f) = 2\,\sigma_f^2(P_{a_{\mathrm{MH}}}) + \mathrm{Var}_\pi(f) $$ です。 次のコードは、20 状態の 2 つ山の分布について、Metropolis と Barker の遷移行列を作り、固有値から漸近分散を厳密に計算して順序を確かめます。あわせて、定常分布から出発する独立な 400 本の連鎖（長さ $N = 2\times10^4$）の時間平均の分散を、 の証明の有限 $N$ の式と比べます。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np n = 20 x = np.arange(n) V = 4.0 * ((x - 9.5) ** 2 / 36.0 - 1.0) ** 2 # 2 つ山（谷は x = 3.5, 15.5、障壁の高さは約 4） beta = 0.5 pi = np.exp(-beta * V); pi /= pi.sum() def chain(accept): # 提案は左右の隣に 1/2 ずつ（端の外への提案は棄却） P = np.zeros((n, n)) for i in range(n): for j in (i - 1, i + 1): if 0 <= j < n: P[i, j] = 0.5 * accept(pi[j] / pi[i]) P[i, i] = 1.0 - P[i].sum() return P def spectrum(P): # pi で対称化して固有分解 d = np.sqrt(pi) lam, U = np.linalg.eigh(d[:, None] * P / d[None, :]) return lam[::-1], (U / d[:, None])[:, ::-1] # 固有値（降順）と l^2(pi) で正規直交な固有関数 def variances(P, f, N): # 漸近分散と、有限の N での N * Var(時間平均) lam, phi = spectrum(P) a2, l = (phi.T @ (pi * f))[1:] ** 2, lam[1:] s_inf = np.sum(a2 * (1 + l) / (1 - l)) s_N = np.sum(a2 * ((1 + l) / (1 - l) - 2 * l * (1 - l**N) / (N * (1 - l) ** 2))) return float(s_inf), float(s_N), 1 - lam[1] f = (x >= 10).astype(float) # 右の山にいる割合 var_pi = float(pi @ f**2 - (pi @ f) ** 2) chains = {"Metropolis": chain(lambda r: min(1.0, r)), "Barker": chain(lambda r: r / (1.0 + r))} print(f"Var_pi(f) = {var_pi:.4f}") rng = np.random.default_rng(20261008) N, R = 20000, 400 for name, P in chains.items(): s_inf, s_N, gap = variances(P, f, N) cum = np.cumsum(P, axis=1); cum[:, -1] = 1.0 X = rng.choice(n, size=R, p=pi) # 定常分布から出発する独立な R 本 tot = np.zeros(R) for _ in range(N): tot += f[X] X = (rng.random(R)[:, None] > cum[X]).sum(axis=1) v = N * (tot / N).var(ddof=1) se = v * np.sqrt(2 / (R - 1)) print(f"{name:10s} gamma = {gap:.5f} sigma_f^2 = {s_inf:7.2f} (2/gamma-1)Var = {(2/gap-1)*var_pi:7.2f}" f" N*Var: 厳密 = {s_N:7.2f}, 標本 = {v:7.2f} ± {se:.2f}") D = chains["Metropolis"] - chains["Barker"] print("非対角成分がすべて Metropolis >= Barker:", bool(np.all(D[~np.eye(n, dtype=bool)] >= 0))) ``` ```text Var_pi(f) = 0.2500 Metropolis gamma = 0.00556 sigma_f^2 = 86.05 (2/gamma-1)Var = 89.66 N*Var: 厳密 = 85.28, 標本 = 82.38 ± 5.83 Barker gamma = 0.00316 sigma_f^2 = 152.09 (2/gamma-1)Var = 157.93 N*Var: 厳密 = 149.69, 標本 = 135.73 ± 9.61 非対角成分がすべて Metropolis >= Barker: True ``` 漸近分散は Metropolis が $86.05$、Barker が $152.09$ で、 の $86.05 \le 152.09 \le 2\times86.05 + 0.25$ が成り立っています。どちらも の上界 $(2/\gamma - 1)\mathrm{Var}_\pi(f)$ の $96\%$ で、「右の山にいるか」という量が最も遅いモードにほぼ一致することを示しています。$\sigma_f^2/\mathrm{Var}_\pi(f) \approx 344$ なので、この連鎖の 344 歩が独立な標本 1 個に相当します。シミュレーションの値は、有限の $N$ での厳密な値から標準誤差の $0.5$ 倍と $1.5$ 倍の位置にあります（標準誤差は、400 個の標本分散のばらつきの正規近似 $\sqrt{2/(R-1)}$ によります）。 採択率が高いほどよい連鎖とは限りません。$S = \{-30, \ldots, 30\}$ の上の、標準偏差 $5$ の正規分布を離散化した $\pi$ に対して、提案を「現在地から $\pm w$ 以内の一様分布」とした Metropolis 法を考え、$f(x) = x$ の漸近分散を固有値から厳密に計算します。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np n = 61 x = np.arange(n) - 30 pi = np.exp(-0.5 * (x / 5.0) ** 2); pi /= pi.sum() # 標準偏差 5 の離散化した正規分布 f = x.astype(float) var_pi = float(pi @ f**2 - (pi @ f) ** 2) print(" 幅 w 採択率 gamma sigma_f^2 / Var_pi(f)") for w in (1, 2, 4, 8, 12, 16, 24, 40): P = np.zeros((n, n)) for i in range(n): for k in range(-w, w + 1): # 提案は i-w, ..., i+w の一様分布（範囲の外は棄却） j = i + k if k != 0 and 0 <= j < n: P[i, j] = min(1.0, pi[j] / pi[i]) / (2 * w + 1) P[i, i] = 1.0 - P[i].sum() acc = float(pi @ (1 - np.diag(P))) d = np.sqrt(pi) lam, U = np.linalg.eigh(d[:, None] * P / d[None, :]) a = (U / d[:, None]).T.dot(pi * f) s2 = float(np.sum(a[:-1] ** 2 * (1 + lam[:-1]) / (1 - lam[:-1]))) print(f" {w:3d} {acc:.3f} {1 - lam[-2]:.5f} {s2 / var_pi:8.2f}") ``` ```text 幅 w 採択率 gamma sigma_f^2 / Var_pi(f) 1 0.613 0.01221 162.55 2 0.705 0.03377 57.88 4 0.715 0.09413 19.83 8 0.622 0.22515 7.31 12 0.518 0.31904 4.58 16 0.429 0.37455 3.71 24 0.304 0.25578 4.08 40 0.185 0.15473 7.05 ``` 幅が小さいと採択率は高い（$w = 4$ で $0.72$）のに、1 歩が小さすぎて漸近分散は大きくなります。幅が大きすぎると、提案が裾に飛んで棄却され、連鎖が動かなくなります。漸近分散が最小になるのは $w = 16$（分布の標準偏差の約 3 倍）で、そのときの採択率は $0.43$ です。幅の違う提案は移動先の集合が違うので では比べられず、このように固有値を計算する（あるいは実際に走らせて次章の方法で漸近分散を推定する）必要があります。 ## 6. ボトルネック比と Cheeger の不等式 実際の MCMC では、状態数が大きすぎて固有値は計算できません。スペクトルギャップを、状態空間の「くびれ」という幾何学的な量で評価するのがこの節です。 たとえば山が 2 つある分布で、山の間の谷を通る確率が小さければ、連鎖は片方の山に長くとどまり、混合は遅くなります。ボトルネック比は、この「出にくさ」を、集合 $A$ から出る確率の最小値として数にしたものです。 $P$ を $\pi$ について可逆な遷移行列とする。$Q(x,y) = \pi(x)P(x,y)$ を辺 $(x,y)$ の**流量**といい、$A, B \subseteq S$ に対して $Q(A, B) = \sum_{x\in A, y\in B}Q(x,y)$ と書く。可逆性から $Q(x,y) = Q(y,x)$ である。$\pi(A) > 0$ のとき $$ \Phi(A) = \frac{Q(A, A^c)}{\pi(A)}, \qquad \Phi_* = \min\bigl\{\Phi(A) : A \subseteq S,\ 0 < \pi(A) \le \tfrac12\bigr\} $$ と置き、$\Phi_*$ を**ボトルネック比**（conductance）という。 $\Phi(A)$ は、定常な連鎖が $A$ の中にいるという条件のもとで、次の 1 歩で $A$ から出る確率です。$\Phi_*$ が小さいことは、出にくい集合があることを意味します。 $|S| \ge 2$ とし、$P$ を正の確率分布 $\pi$ について可逆で既約な $S$ 上の遷移行列とする。このとき $$ \frac{\Phi_*^2}{2} \le \gamma \le 2\Phi_*. $$ Dirichlet 形式を $\mathcal{E}(f,f) = \langle f, (I-P)f\rangle_\pi = \frac12\sum_{x,y}Q(x,y)(f(x) - f(y))^2$ と書きます（）。$\gamma = \min\mathcal{E}(f,f)/\mathrm{Var}_\pi(f)$ です（ の 1）。 **上界 $\gamma \le 2\Phi_*$。** $0 < \pi(A) \le \frac12$ なる $A$ を取り、$f = \mathbf{1}_A$ とします。$(f(x) - f(y))^2$ は、$x, y$ の一方だけが $A$ に属するとき $1$、それ以外は $0$ なので、$\mathcal{E}(f,f) = \frac12(Q(A, A^c) + Q(A^c, A)) = Q(A, A^c)$ です。$\mathrm{Var}_\pi(f) = \pi(A)(1 - \pi(A)) \ge \frac12\pi(A)$ なので、$\gamma \le Q(A, A^c)/(\frac12\pi(A)) = 2\Phi(A)$ です。$A$ について最小を取ります。 **下界 $\gamma \ge \Phi_*^2/2$。** 2 段階に分けます。 *第 1 段：台の小さい非負の関数への帰着。* の $\varphi_2$ を取ります。$\varphi_2$ は $\mathbf{1}$ と直交し零でないので、正の値と負の値の両方を取ります。$\{\varphi_2 > 0\}$ と $\{\varphi_2 < 0\}$ は交わらないので、どちらかは $\pi$ の測度が $\frac12$ 以下です。必要なら $\varphi_2$ を $-\varphi_2$ に替えて、$\pi(\{\varphi_2 > 0\}) \le \frac12$ とします。$f = \max(\varphi_2, 0)$ と置きます。$f \ge 0$、$f \ne 0$ です。$f \ge \varphi_2$ から $Pf \ge P\varphi_2 = \lambda_2\varphi_2$ なので、$f(x) > 0$ である点では $(Pf)(x) \ge \lambda_2f(x)$、すなわち $$ \bigl((I - P)f\bigr)(x) \le (1 - \lambda_2)f(x) = \gamma f(x) \qquad (f(x) > 0) $$ です。$f(x) = 0$ の点は内積に寄与しないので、$\mathcal{E}(f,f) = \sum_{x : f(x) > 0}\pi(x)f(x)\bigl((I-P)f\bigr)(x) \le \gamma\langle f, f\rangle_\pi$ です。したがって $\gamma \ge \mathcal{E}(f,f)/\langle f, f\rangle_\pi$ で、次を示せば十分です。 *第 2 段：$f \ge 0$、$f \ne 0$、$\pi(\{f > 0\}) \le \frac12$ ならば $\mathcal{E}(f,f) \ge \frac12\Phi_*^2\langle f, f\rangle_\pi$。* $g = f^2$ と置き、$t \ge 0$ に対して $A_t = \{x : g(x) > t\}$ とします。$A_t \subseteq \{f > 0\}$ なので $\pi(A_t) \le \frac12$ で、$A_t$ が空でなければ $Q(A_t, A_t^c) \ge \Phi_*\pi(A_t)$ です（空なら両辺 $0$）。$g(x) > g(y)$ のとき $g(x) - g(y) = \int_0^\infty\mathbf{1}\{g(y) \le t < g(x)\}\,dt$ なので、 $$ \sum_{(x,y)\,:\,g(x) > g(y)}Q(x,y)\bigl(g(x) - g(y)\bigr) = \int_0^\infty Q(A_t, A_t^c)\,dt \ge \Phi_*\int_0^\infty\pi(A_t)\,dt = \Phi_*\sum_x\pi(x)g(x) = \Phi_*\langle f, f\rangle_\pi $$ です（$x \in A_t$ かつ $y \notin A_t$ は $g(y) \le t < g(x)$ と同値。最後から 2 番目の等号は $\int_0^\infty\mathbf{1}\{t < g(x)\}dt = g(x)$）。一方、$g(x) - g(y) = (f(x) - f(y))(f(x) + f(y))$ と Cauchy–Schwarz の不等式から、左辺は $$ \Bigl(\sum_{f(x) > f(y)}Q(x,y)(f(x) - f(y))^2\Bigr)^{1/2}\Bigl(\sum_{f(x) > f(y)}Q(x,y)(f(x) + f(y))^2\Bigr)^{1/2} $$ 以下です。$Q$ が対称なので、第 1 の和は $\frac12\sum_{x,y}Q(x,y)(f(x)-f(y))^2 = \mathcal{E}(f,f)$ に等しく、第 2 の和は $\frac12\sum_{x,y}Q(x,y)(f(x)+f(y))^2$ 以下です。$(f(x) + f(y))^2 \le 2f(x)^2 + 2f(y)^2$ と $\sum_yQ(x,y) = \pi(x)$、$\sum_xQ(x,y) = \pi(y)$ から、 $$ \frac12\sum_{x,y}Q(x,y)(f(x) + f(y))^2 \le \sum_{x,y}Q(x,y)\bigl(f(x)^2 + f(y)^2\bigr) = 2\langle f, f\rangle_\pi $$ です。まとめると $\Phi_*\langle f, f\rangle_\pi \le \sqrt{\mathcal{E}(f,f)}\sqrt{2\langle f, f\rangle_\pi}$ で、両辺を二乗して $\langle f,f\rangle_\pi > 0$ で割れば主張を得ます。 （下界の証明の方針は、固有関数の正の部分を取り出し、$f^2$ の水準集合 $A_t$ ごとに流量の評価 $Q(A_t, A_t^c) \ge \Phi_*\pi(A_t)$ を使って足し合わせる、というものです。）上界は、くびれが 1 つでもあればギャップが小さいことを示す**否定的な**結果に使います。下界は、すべての集合について流量を下から評価できれば混合が速いことを示す**肯定的な**結果に使います（下界は 2 乗で効くので、一般には粗い評価です）。 $S = \{0, 1, \ldots, L\}$ とし、$P$ を正の確率分布 $\pi$ について可逆で既約な遷移行列で、$|x - y| \ge 2$ ならば $P(x,y) = 0$（**局所的**）であるものとする。$1 \le m \le L$ が $A = \{0, \ldots, m-1\}$ について $\pi(A) \le \frac12$ を満たすならば $$ \gamma \le 2\Phi_* \le \frac{2\pi(m)}{\pi(A)} $$ である。特に $\pi(x) \propto e^{-\beta V(x)}$ で、$\Delta = V(m) - \min_{x < m}V(x)$ と置くと、$\gamma \le 2e^{-\beta\Delta}$ であり、$\gamma^* > 0$ かつ $0 < \varepsilon < \frac12$ のとき $t_{\mathrm{mix}}(\varepsilon) \ge \bigl(\frac12e^{\beta\Delta} - 1\bigr)\ln\frac{1}{2\varepsilon}$ である。 局所性から、$A$ と $A^c = \{m, \ldots, L\}$ をつなぐ辺は $(m-1, m)$ だけなので、$Q(A, A^c) = Q(m-1, m) = Q(m, m-1) = \pi(m)P(m, m-1) \le \pi(m)$ です。よって $\Phi_* \le \Phi(A) \le \pi(m)/\pi(A)$ で、 の上界から最初の主張を得ます。 $V$ の $A$ の上での最小点を $x_A$ とすると、$\pi(A) \ge \pi(x_A)$ から $\pi(m)/\pi(A) \le \pi(m)/\pi(x_A) = e^{-\beta(V(m) - V(x_A))} = e^{-\beta\Delta}$ です。$\gamma^* \le \gamma \le 2e^{-\beta\Delta}$ なので、$\gamma^* < 1$ なら の 3 から $t_{\mathrm{mix}}(\varepsilon) \ge (\frac{1}{\gamma^*} - 1)\ln\frac{1}{2\varepsilon} \ge (\frac12e^{\beta\Delta} - 1)\ln\frac{1}{2\varepsilon}$ です（$\gamma^* = 1$ は $2e^{-\beta\Delta} \ge 1$ のときにしか起きず、そのとき右辺は $0$ 以下です）。 $\Delta$ は、左の山の底から見た障壁の高さです。隣にしか動けない連鎖は、採択関数をどう工夫しても、障壁の上の状態 $m$ を通らなければ山を移れません。その状態の確率が $e^{-\beta\Delta}$ の程度しかないので、混合に $e^{\beta\Delta}$ の程度の時間がかかります。低温（大きい $\beta$）の統計力学の模型や、多峰の事後分布で MCMC が「片方の山に留まったまま収束したように見える」のはこのためです。次のコードは、14 状態の 2 つ山の分布で、$\beta$ を増やしたときのギャップとボトルネック比（全部分集合についての最小値）を計算します。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np n = 14 x = np.arange(n) V = ((x - 6.5) ** 2 / 16.0 - 1.0) ** 2 # 谷は x = 2.5, 10.5 付近、障壁の高さ 1（x = 6, 7 で V = 0.969） def metropolis(pi): P = np.zeros((n, n)) for i in range(n): for j in (i - 1, i + 1): if 0 <= j < n: P[i, j] = 0.5 * min(1.0, pi[j] / pi[i]) P[i, i] = 1.0 - P[i].sum() return P masks = ((np.arange(1, 2**n - 1)[:, None] >> np.arange(n)) & 1).astype(float) # 空でも全体でもない部分集合 print(" beta Phi*^2/2 gamma 2 Phi* 2 pi(m)/pi(A)") for beta in (1, 2, 4, 6, 8): pi = np.exp(-beta * V); pi /= pi.sum() P = metropolis(pi) d = np.sqrt(pi) lam = np.linalg.eigvalsh(d[:, None] * P / d[None, :]) gamma = 1 - lam[-2] F = pi[:, None] * P # 辺の流量 Q(x, y) = pi(x) P(x, y) piA = masks.dot(pi) flow = np.einsum("ai,ij,aj->a", masks, F, 1 - masks) ok = piA <= 0.5 + 1e-12 phi = (flow[ok] / piA[ok]).min() m = 7 # 切り口 A = {0,...,6} print(f"{beta:5.1f} {phi**2/2:.3e} {gamma:.3e} {2*phi:.3e} {2*pi[m]/pi[:m].sum():.3e}") ``` ```text beta Phi*^2/2 gamma 2 Phi* 2 pi(m)/pi(A) 1.0 1.156e-03 2.110e-02 9.617e-02 1.923e-01 2.0 3.255e-04 1.191e-02 5.103e-02 1.021e-01 4.0 1.492e-05 2.829e-03 1.092e-02 2.185e-02 6.0 4.917e-07 5.608e-04 1.983e-03 3.967e-03 8.0 1.434e-08 1.017e-04 3.387e-04 6.774e-04 ``` すべての $\beta$ で $\Phi_*^2/2 \le \gamma \le 2\Phi_* \le 2\pi(m)/\pi(A)$ が成り立っています。$\gamma$ は上界 $2\Phi_*$ の $0.2$〜$0.3$ 倍で、下界よりずっと上界に近く、$\beta$ を $2$ 増やすごとにおよそ $1/5$ になります（この例では $\Delta = 0.969 - 0.055 = 0.914$ で、$e^{-2\Delta} \approx 0.16$）。最小を与える集合は左半分 $\{0, \ldots, 6\}$ で、最後の列が $2\Phi_*$ のちょうど 2 倍なのは、$Q(A, A^c) = \pi(7)P(7,6) = \frac12\pi(7)$ だからです。 ## 7. カップリング 混合時間の上界を、固有値を使わずに得る方法です。同じ遷移行列に従う 2 つの連鎖を、**同じ乱数を使って**同時に動かし、出会うまでの時間を評価します。 遷移行列 $P$ の**カップリング**とは、$S \times S$ 上の遷移行列 $K$ で、すべての $x, y, x', y' \in S$ について $$ \sum_{y'}K\bigl((x,y), (x',y')\bigr) = P(x, x'), \qquad \sum_{x'}K\bigl((x,y), (x',y')\bigr) = P(y, y') $$ を満たし、かつ $x' \ne y'$ ならば $K((x,x), (x',y')) = 0$（一度出会ったら離れない）となるものをいう。$K$ に従う連鎖 $(X_t, Y_t)$ について、$\tau = \min\{t \ge 0 : X_t = Y_t\}$ を**合流時刻**という。 $K$ を $P$ のカップリング、$(X_t, Y_t)$ を $(x, y)$ から出発する $K$ の連鎖とする。 1. $(X_t)$ は $x$ から出発する遷移行列 $P$ の Markov 連鎖、$(Y_t)$ は $y$ から出発する遷移行列 $P$ の Markov 連鎖である。 2. すべての $t \ge 0$ で $\|P^t(x, \cdot) - P^t(y, \cdot)\|_{\mathrm{TV}} \le \mathbb{P}_{(x,y)}(\tau > t)$。 3. $P$ が定常分布 $\pi$ をもつとき、$d(t) \le \max_{x,y}\mathbb{P}_{(x,y)}(\tau > t)$。 **1.** $K$ の連鎖の道の確率の公式から $$ \mathbb{P}(X_0 = x_0, \ldots, X_t = x_t) = \sum_{y_1, \ldots, y_t}\prod_{s=1}^tK\bigl((x_{s-1}, y_{s-1}), (x_s, y_s)\bigr) \qquad (x_0 = x,\ y_0 = y) $$ です。$y_t$ は最後の因子にしか現れないので、先に $y_t$ について和を取ると、周辺の条件から最後の因子は $P(x_{t-1}, x_t)$ になります。これは $y_{t-1}$ によらないので、次に $y_{t-1}$ について和を取ることができ、順に繰り返して $\prod_{s=1}^tP(x_{s-1}, x_s)$ を得ます。これは $x$ から出発する $P$ の連鎖の道の確率で、条件付き確率を計算すると Markov 連鎖の定義（）の式が成り立ちます。$(Y_t)$ も同様です。 **2.** カップリングの第 3 の条件から、$X_s = Y_s$ なら確率 $1$ で $X_{s+1} = Y_{s+1}$ なので、帰納的に、$t \ge \tau$ では $X_t = Y_t$ です。したがって $\{X_t \ne Y_t\} \subseteq \{\tau > t\}$ です（確率 $0$ の事象を除いて）。1 から $\mathbb{P}(X_t \in A) = P^t(x, A)$、$\mathbb{P}(Y_t \in A) = P^t(y, A)$ なので、任意の $A \subseteq S$ について $$ P^t(x, A) - P^t(y, A) = \mathbb{P}(X_t \in A) - \mathbb{P}(Y_t \in A) \le \mathbb{P}(X_t \in A,\ Y_t \notin A) \le \mathbb{P}(X_t \ne Y_t) \le \mathbb{P}(\tau > t) $$ です。$A$ について最大を取ると、 の 1 から左辺は全変動距離になります。 **3.** の 3 と 2 から従います。 の超立方体 $\{0,1\}^n$ 上の連鎖について、すべての $t \ge 0$ で $d(t) \le n(1 - \frac1n)^t \le ne^{-t/n}$ である。したがって $0 < \varepsilon < 1$ について $$ t_{\mathrm{mix}}(\varepsilon) \le \bigl\lceil n\ln n + n\ln(1/\varepsilon)\bigr\rceil. $$ $(I_t, B_t)$（$t = 1, 2, \ldots$）を、$\{1, \ldots, n\} \times \{0, 1\}$ 上の一様分布に従う独立な確率変数の列とします。$(X_0, Y_0) = (x, y)$ から、時刻 $t$ に $X_{t-1}$ と $Y_{t-1}$ の**両方の**第 $I_t$ 成分を $B_t$ で置き換えて $(X_t, Y_t)$ を作ります。$(X_t, Y_t)$ は、現在の状態と独立な乱数で次の状態が決まるので $S \times S$ 上の Markov 連鎖で、その遷移行列 $K$ は、片方だけを見れば元の連鎖の 1 歩なので周辺の条件を満たし、$X_{t-1} = Y_{t-1}$ なら $X_t = Y_t$ です。よって $K$ はカップリングです。 座標 $i$ が時刻 $s$ に選ばれる（$I_s = i$）と、$X_s$ と $Y_s$ の第 $i$ 成分はともに $B_s$ になり、以後は同じ値で置き換えられ続けるので、ずっと一致します。したがって、すべての座標が少なくとも 1 回選ばれた時点で $X = Y$ となり、$\{\tau > t\} \subseteq \bigcup_{i=1}^n\{I_1 \ne i, \ldots, I_t \ne i\}$ です。各事象の確率は $(1 - \frac1n)^t$ なので、 $$ \mathbb{P}_{(x,y)}(\tau > t) \le n\Bigl(1 - \frac1n\Bigr)^t \le ne^{-t/n} $$ です。 の 3 から $d(t) \le ne^{-t/n}$ で、$t \ge n\ln n + n\ln(1/\varepsilon)$ なら右辺は $\varepsilon$ 以下です。 証明に出てくる「すべての座標が選ばれるまでの時間」は、クーポン収集問題の待ち時間です。スペクトルによる上界（）が $n^2$ の程度だったのに対して、カップリングは $n\ln n$ の程度を与えます。差は、 の $\ln(1/\pi_{\min}) = n\ln2$ の因子から来ています。次のコードは、$n = 8, 10$ で全変動距離を行列の冪から厳密に計算し、2 つの上界と比べます。 ```python # Python 3.9.6, NumPy 2.0.2 import numpy as np def hypercube(n): # 座標を 1 つ選び、一様なビットで置き換える N = 2**n P = np.zeros((N, N)) for s in range(N): P[s, s] += 0.5 for i in range(n): P[s, s ^ (1 << i)] += 0.5 / n return P def first(bound, eps=0.25): t = 0 while bound(t) > eps: t += 1 return t for n in (8, 10): N, P = 2**n, hypercube(n) Pt, d = np.eye(N), [] for t in range(81): d.append(0.5 * np.abs(Pt - 1.0 / N).sum(axis=1).max()) Pt = Pt.dot(P) spec = lambda t: 0.5 * np.sqrt(N - 1) * (1 - 1 / n) ** t coup = lambda t: n * (1 - 1 / n) ** t print(f"n = {n}") print(" t d(t) 厳密 スペクトルの上界 カップリングの上界") for t in (10, 20, 30, 40, 60, 80): print(f" {t:3d} {d[t]:.5f} {min(spec(t), 1):.5f} {min(coup(t), 1):.5f}") print(f" t_mix(1/4): 厳密 = {next(t for t in range(81) if d[t] <= 0.25)}, スペクトル <= {first(spec)}," f" カップリング <= {first(coup)}, 下界 (1/gamma* - 1) ln 2 = {(n - 1) * np.log(2):.2f}") ``` ```text n = 8 t d(t) 厳密 スペクトルの上界 カップリングの上界 10 0.30782 1.00000 1.00000 20 0.07734 0.55259 0.55367 30 0.02001 0.14537 0.14566 40 0.00524 0.03824 0.03832 60 0.00036 0.00265 0.00265 80 0.00003 0.00018 0.00018 t_mix(1/4): 厳密 = 12, スペクトル <= 26, カップリング <= 26, 下界 (1/gamma* - 1) ln 2 = 4.85 n = 10 t d(t) 厳密 スペクトルの上界 カップリングの上界 10 0.44288 1.00000 1.00000 20 0.15534 1.00000 1.00000 30 0.05289 0.67793 0.42391 40 0.01827 0.23638 0.14781 60 0.00221 0.02874 0.01797 80 0.00027 0.00349 0.00218 t_mix(1/4): 厳密 = 16, スペクトル <= 40, カップリング <= 36, 下界 (1/gamma* - 1) ln 2 = 6.24 ``` どの $t$ でも厳密な $d(t)$ は 2 つの上界の下にあり、混合時間は の 3 の下界より上にあります。2 つの上界はどちらも $(1 - \frac1n)^t$ に比例し、係数がスペクトルでは $\frac12\sqrt{2^n - 1}$、カップリングでは $n$ です。$n = 8$ では係数が $7.98$ と $8$ でたまたまほぼ一致し、$n = 10$ では $16.0$ と $10$ で、$n$ が大きいほどカップリングの優位が広がります。 ## 8. 走査の順序 多成分の系では、どの成分をどの順に更新するかを決める必要があります。熱浴法（）で 2 つの流儀を比べます。 の設定で次が成り立つ。 1. 各 $P_i$ は $\pi$ について可逆で、$P_i^2 = P_i$ である。したがってすべての $f$ で $\langle f, P_if\rangle_\pi = \langle P_if, P_if\rangle_\pi \ge 0$。 2. 無作為走査 $P_{\mathrm{RS}}$ は $\pi$ について可逆で、すべての $f$ で $\langle f, P_{\mathrm{RS}}f\rangle_\pi \ge 0$ である。したがって固有値はすべて $[0, 1]$ に属し、$\gamma^* = \gamma$ である。 3. 系統走査 $P_{\mathrm{SS}}$ は $\pi$ を定常分布にもつ。$n = 2$ のとき、$P_{\mathrm{SS}} = P_1P_2$ が $\pi$ について可逆であることは $P_1P_2 = P_2P_1$ と同値であり、一般には可逆でない。 $\pi_i(x) = \sum_{s \in S_i}\pi(x^{i\to s})$ と置きます。$\pi_i(x)$ は $x$ の第 $i$ 成分によりません。$x$ と $y$ が第 $i$ 成分以外で一致することを $x \sim_i y$ と書くと、$P_i(x,y) = \pi(y)/\pi_i(x)$（$x \sim_i y$）で、このとき $\pi_i(x) = \pi_i(y)$ です。 **1.** $x \sim_i y$ なら $\pi(x)P_i(x,y) = \pi(x)\pi(y)/\pi_i(x)$ で、$x, y$ について対称です。$x \not\sim_i y$ なら両辺 $0$ です。よって可逆です。$x \sim_i z$ について $$ (P_iP_i)(x,z) = \sum_{y \sim_i x}\frac{\pi(y)}{\pi_i(x)}\cdot\frac{\pi(z)}{\pi_i(y)} = \frac{\pi(z)}{\pi_i(x)}\cdot\frac{\sum_{y\sim_ix}\pi(y)}{\pi_i(x)} = \frac{\pi(z)}{\pi_i(x)} = P_i(x,z) $$ です（$x \not\sim_i z$ なら両辺 $0$）。$P_i$ は自己共役なので、$\langle f, P_if\rangle_\pi = \langle f, P_iP_if\rangle_\pi = \langle P_if, P_if\rangle_\pi$ です。 **2.** 可逆な遷移行列の平均は可逆で、$\langle f, P_{\mathrm{RS}}f\rangle_\pi = \frac1n\sum_i\langle P_if, P_if\rangle_\pi \ge 0$ です。 の固有関数について $\lambda_j = \langle\varphi_j, P_{\mathrm{RS}}\varphi_j\rangle_\pi \ge 0$ なので、$\lambda_* = \lambda_2$ です。 **3.** 1 と から各 $i$ で $\pi P_i = \pi$ なので、$\pi P_1\cdots P_n = \pi$ です。$n = 2$ とします。可逆性は $\ell^2(\pi)$ での自己共役性と同値で、$P_1, P_2$ は自己共役なので、$\langle f, P_1P_2g\rangle_\pi = \langle P_2P_1f, g\rangle_\pi$、すなわち $P_1P_2$ の共役は $P_2P_1$ です。よって $P_1P_2$ が可逆であることは $P_1P_2 = P_2P_1$ と同値です。 可逆でない例を挙げます。$S = \{0,1\}^2$、$\pi(00) = \pi(11) = 0.4$、$\pi(01) = \pi(10) = 0.1$ とします（$x_1x_2$ と並べて書きます）。$P_1$ は $x_1$ を、$P_2$ は $x_2$ を引き直します。 $$ P_1(00, 00) = \frac{0.4}{0.4 + 0.1} = 0.8, \qquad P_2(00, 01) = \frac{0.1}{0.4 + 0.1} = 0.2, \qquad P_1(01, 01) = \frac{0.1}{0.1 + 0.4} = 0.2, \qquad P_2(01, 00) = 0.8 $$ です。$00$ から $01$ へは「$P_1$ で動かず、$P_2$ で $x_2$ を変える」道しかなく（$P_1$ で $10$ に移ると、$P_2$ では $x_1 = 1$ のままです）、$(P_1P_2)(00, 01) = 0.8 \times 0.2 = 0.16$ です。同様に $(P_1P_2)(01, 00) = P_1(01,01)P_2(01,00) = 0.2 \times 0.8 = 0.16$ です。$\pi(00) \times 0.16 = 0.064$ と $\pi(01) \times 0.16 = 0.016$ は等しくないので、詳細つり合いは成り立ちません。 系統走査は可逆でないので、この章の §4〜§6 の定理（固有値が実数であることに依存しています）はそのままでは使えません。それでも $\pi$ が定常分布であることは変わらず、既約であれば から時間平均は正しい値に収束します。実用の格子計算では、乱数の消費が少なくメモリの参照が規則的な系統走査がよく使われます。可逆な連鎖がほしいときは、$P_1\cdots P_nP_n\cdots P_1$ のように往復で走査します（共役が自分自身になります）。 周期境界の 1 次元 Ising 模型は、$S = \{-1, +1\}^n$、$\pi(s) \propto \exp\bigl(\beta\sum_{i=1}^ns_is_{i+1}\bigr)$（$s_{n+1} = s_1$）です。スピン $i$ を含む項は $\beta s_ih_i$（$h_i = s_{i-1} + s_{i+1} \in \{-2, 0, 2\}$）だけなので、熱浴更新は $$ \mathbb{P}(s_i = +1 \mid \text{他のスピン}) = \frac{e^{\beta h_i}}{e^{\beta h_i} + e^{-\beta h_i}} = \frac{1}{1 + e^{-2\beta h_i}} $$ で $s_i$ を引き直します。現在の値と逆の値になる確率は、$r = \pi(s^{(i)})/\pi(s) = e^{-2\beta s_ih_i}$（$s^{(i)}$ はスピン $i$ を反転した配位）として $r/(1 + r)$ です。これは「スピン $i$ の反転」を提案とした **Barker の採択関数**にほかなりません。したがって、無作為走査の熱浴法は、提案 $Q(s, s^{(i)}) = 1/n$ と採択関数 $a_{\mathrm{B}}$ による連鎖で、同じ提案に Metropolis の採択関数 $\min(1, r)$ を使った連鎖と で比べられます。2 値のスピンでは、どの量についても Metropolis 法の漸近分散は熱浴法以下で、熱浴法は高々 2 倍（と $\mathrm{Var}_\pi(f)$）しか悪くなりません。成分が 3 つ以上の値を取る模型では、熱浴法は反転より多くの移動先をもつので、この比較は当てはまりません。 連続な状態空間では、遷移行列が遷移核に、固有値の有限和がスペクトル分解に替わります。定理の形は似ていますが、スペクトルギャップが $0$ になりうる、1 点の確率が $0$ なので の $\pi_{\min}$ が使えない、といった違いが生じます。 - **収束の理論。** 既約性にあたる条件（Harris 再帰性）と、ドリフト条件・小集合の条件による幾何エルゴード性の理論は、Meyn–Tweedie の本にまとめられています。 の一般の状態空間への拡張は Tierney によります。 - **中心極限定理。** 可逆な連鎖で $\sigma_f^2 < \infty$ なら、$\sqrt N(\frac1N\sum f(X_k) - \pi(f))$ は $N(0, \sigma_f^2)$ に分布収束します（Kipnis–Varadhan）。証明はマルチンゲール近似によるもので、ここでは定理として扱いません。時間平均に誤差を付ける実際の方法は次章で扱います。 - **最適な採択率。** $d$ 次元の積の形の目標に対する無作為歩行 Metropolis 法では、$d \to \infty$ で提案の幅を $d^{-1/2}$ に比例させるのが最適で、そのときの採択率は約 $0.234$ になります（Roberts–Gelman–Gilks）。 の 1 次元の例で最適な採択率が $0.43$ だったことと矛盾しません（低次元では最適値はもっと高くなります）。 - **ボトルネックへの対策。** の遅さは局所的な更新に共通です。温度の違う複数の連鎖を走らせて配位を交換する方法（レプリカ交換法）や、模型の構造を使った大域的な更新（クラスター更新）は、移動先の集合そのものを変えて障壁を迂回します。 ## 9. 演習 $P$ を可逆で既約な遷移行列、$P_L = \frac12(I + P)$ を怠惰な連鎖とする。すべての $f$ について $\sigma_f^2(P_L) = 2\sigma_f^2(P) + \mathrm{Var}_\pi(f)$ を示せ。これは と整合するか。 $P_L$ は $P$ と同じ固有関数 $\varphi_j$ をもち、固有値は $\mu_j = \frac12(1 + \lambda_j)$ です。 の係数は $$ \frac{1 + \mu_j}{1 - \mu_j} = \frac{3 + \lambda_j}{1 - \lambda_j} = \frac{2(1 + \lambda_j)}{1 - \lambda_j} + 1 $$ なので、$a_j^2$ を掛けて $j \ge 2$ について足すと $\sigma_f^2(P_L) = 2\sigma_f^2(P) + \mathrm{Var}_\pi(f)$ です。$x \ne y$ で $P(x,y) \ge \frac12P(x,y) = P_L(x,y)$ なので、 は $\sigma_f^2(P) \le \sigma_f^2(P_L)$ を予言し、上の等式はそれを定量的にしたものです。怠惰な連鎖は半分の時間を止まって過ごすので、同じ精度に約 2 倍の歩数が要ります。理論の都合（負の固有値を消す）のための道具で、実際の計算でわざわざ怠惰にする理由はありません。 $S = \{0, 1\}$、$\pi = (1 - p, p)$（$0 < p \le \frac12$）とする。(1) 提案を「反転」（$Q(0,1) = Q(1,0) = 1$）とした Metropolis 法の遷移行列 $P_{\mathrm{M}}$ と、熱浴法（毎回 $\pi$ から引き直す）の遷移行列 $P_{\mathrm{H}}$ を書け。(2) それぞれの $\lambda_2$ と、$\sigma_f^2/\mathrm{Var}_\pi(f)$ を求めよ。(3) 独立な標本より精度がよくなる理由を説明せよ。 (1) $r(0,1) = p/(1-p) \le 1$、$r(1,0) = (1-p)/p \ge 1$ なので $$ P_{\mathrm{M}} = \begin{pmatrix} 1 - \frac{p}{1-p} & \frac{p}{1-p} \\ 1 & 0\end{pmatrix}, \qquad P_{\mathrm{H}} = \begin{pmatrix} 1 - p & p \\ 1 - p & p\end{pmatrix} $$ です。$P_{\mathrm{H}}$ は Barker の採択関数（$r/(1+r)$ は $0 \to 1$ で $p$、$1 \to 0$ で $1 - p$）による連鎖でもあります。 (2) から $\lambda_2 = 1 - a - b$ です。$P_{\mathrm{H}}$ では $a + b = 1$ で $\lambda_2 = 0$、$\sigma_f^2/\mathrm{Var}_\pi(f) = 1$（独立な標本そのものです）。$P_{\mathrm{M}}$ では $\lambda_2 = 1 - \frac{p}{1-p} - 1 = -\frac{p}{1-p}$ で、 $$ \frac{\sigma_f^2}{\mathrm{Var}_\pi(f)} = \frac{1 + \lambda_2}{1 - \lambda_2} = \frac{(1 - 2p)/(1-p)}{1/(1-p)} = 1 - 2p $$ です。 (3) Metropolis 法では、状態 $1$ にいると必ず $0$ に戻るので、隣り合う標本が負に相関します。 の式で $\lambda_2 < 0$ の係数が $1$ より小さくなる場合です。$p = \frac12$ では $\lambda_2 = -1$、連鎖は $0, 1, 0, 1, \ldots$ と決定的に交代し、時間平均は（偶数歩で）正確に $\frac12$ になって分散は $0$ です。このとき連鎖は周期的で分布は収束しませんが（$\gamma^* = 0$）、時間平均の精度は最良です。分布の収束と時間平均の精度が別の問題であることを示す例です。 $S = \{0, 1, 2\}$、 $$ P = \begin{pmatrix} \frac12 & \frac12 & 0 \\ \frac14 & \frac12 & \frac14 \\ 0 & \frac12 & \frac12\end{pmatrix} $$ について、定常分布、固有値、$\Phi_*$ を求め、 の両側の不等式を確かめよ。 $\pi = (\frac14, \frac12, \frac14)$ が詳細つり合い $\pi(0)P(0,1) = \frac18 = \pi(1)P(1,0)$、$\pi(1)P(1,2) = \frac18 = \pi(2)P(2,1)$ を満たします。固有関数は $\mathbf{1}$（固有値 $1$）、$(1, 0, -1)$（$P(1,0,-1)^{\mathsf{T}} = (\frac12, 0, -\frac12)^{\mathsf{T}}$ で固有値 $\frac12$）、$(1, -1, 1)$（$P(1,-1,1)^{\mathsf{T}} = (0, 0, 0)^{\mathsf{T}}$ で固有値 $0$）です。よって $\gamma = \frac12$ です。 $\pi(A) \le \frac12$ となる空でない $A$ は $\{0\}, \{1\}, \{2\}, \{0, 2\}$ です。$\Phi(\{0\}) = \frac{1/8}{1/4} = \frac12$、$\Phi(\{2\}) = \frac12$、$\Phi(\{1\}) = \frac{\pi(1)(P(1,0) + P(1,2))}{\pi(1)} = \frac12$、$\Phi(\{0,2\}) = \frac{1/8 + 1/8}{1/2} = \frac12$ なので $\Phi_* = \frac12$ です。不等式は $\Phi_*^2/2 = \frac18 \le \gamma = \frac12 \le 2\Phi_* = 1$ で、成り立っています。 $S = \{1, \ldots, q\}^n$（$q \ge 2$）、$\pi = \nu^{\otimes n}$（$\nu$ は $\{1, \ldots, q\}$ 上の正の確率分布）に対する無作為走査の熱浴法の混合時間が、$q$ と $\nu$ によらず $t_{\mathrm{mix}}(\varepsilon) \le \lceil n\ln n + n\ln(1/\varepsilon)\rceil$ を満たすことを示せ。 積測度では、成分 $i$ の条件付き分布は他の成分によらず $\nu$ なので、連鎖の 1 歩は「座標 $I$ を一様に選び、その成分を $\nu$ に従う $B$ で置き換える」です。 の証明と同じく、2 つの連鎖に同じ $(I_t, B_t)$ を使うカップリングを作ると、座標 $i$ が一度選ばれた後は両者の第 $i$ 成分が一致し続けます。合流時刻はすべての座標が選ばれる時刻以下で、$\mathbb{P}(\tau > t) \le n(1 - \frac1n)^t \le ne^{-t/n}$ です。 の 3 から主張が従います。一方 の 2 では $\pi_{\min} = (\min\nu)^n$ が現れ、$\nu$ が偏っているほど上界が悪くなります。カップリングの評価にはこの依存がありません。成分の間に相互作用がある（$\pi$ が積でない）と、同じ乱数を使っても条件付き分布が 2 つの連鎖で異なるので、この単純な議論は使えなくなります。 ## 参考文献 1. D. A. Levin, Y. Peres, E. L. Wilmer, *Markov Chains and Mixing Times*, American Mathematical Society, 2009（第 2 版 2017）— 第 4〜5 章（全変動距離、カップリング）、第 7 章（下界とボトルネック比）、第 12〜13 章（固有値、Dirichlet 形式、Cheeger の不等式）。 2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, "Equation of state calculations by fast computing machines", *Journal of Chemical Physics* 21 (1953), 1087–1092. 3. W. K. Hastings, "Monte Carlo sampling methods using Markov chains and their applications", *Biometrika* 57 (1970), 97–109. 4. P. H. Peskun, "Optimum Monte-Carlo sampling using Markov chains", *Biometrika* 60 (1973), 607–612. 5. L. Tierney, "A note on Metropolis–Hastings kernels for general state spaces", *Annals of Applied Probability* 8 (1998), 1–9. 6. A. Sinclair, M. Jerrum, "Approximate counting, uniform generation and rapidly mixing Markov chains", *Information and Computation* 82 (1989), 93–133. 7. G. O. Roberts, A. Gelman, W. R. Gilks, "Weak convergence and optimal scaling of random walk Metropolis algorithms", *Annals of Applied Probability* 7 (1997), 110–120. 8. C. Kipnis, S. R. S. Varadhan, "Central limit theorem for additive functionals of reversible Markov processes and applications to simple exclusions", *Communications in Mathematical Physics* 104 (1986), 1–19. 9. S. P. Meyn, R. L. Tweedie, *Markov Chains and Stochastic Stability*, 2nd ed., Cambridge University Press, 2009.