CS236 補充課 D · 前置數學

Markov 鏈與 MCMC：算不出機率，也能取樣

Energy-based model 只能算出機率的相對大小，無法直接取樣。MCMC 的想法是：設計一個隨機漫步，讓它走久了之後，停留在每個地方的時間比例剛好等於目標機率。這一頁從條件機率和一個天氣的例子開始，一步一步走到 Metropolis–Hastings 和 Gibbs sampling，每條公式都會推導出來。

先修：補充課 A（特徵值、變異數、取條件）服務：第 11–13 講預計閱讀 2–3 小時，可以分次讀

讀完這頁你應該能夠

說明為什麼只知道相對機率也能取樣。
手算 Markov 鏈走一步後的分佈，並解出穩態分佈。
說明第二大特徵值為什麼決定收斂速度。
寫出 detailed balance，證明它保證目標分佈是穩態。
從 detailed balance 推導 Metropolis–Hastings 的接受率，並說明為什麼不需要歸一化常數。
解釋 burn-in、自相關和有效樣本數，並推導有效樣本數的公式。
解釋 Gibbs sampling 的步驟，以及為什麼它在變數高度相關時很慢。
解釋多峰分佈為什麼會困住鏈，以及提高溫度為什麼有幫助。

00為什麼需要 MCMC

第 11 講的 energy-based model 定義

$$\begin{aligned}p_\theta(\mathbf x) &= \frac1Ze^{-E_\theta(\mathbf x)}\\Z &= \int e^{-E_\theta(\mathbf x)}\,d\mathbf x\end{aligned}$$

$E_\theta$ 是一個神經網路，隨時可以算。但 $Z$ 是對所有可能的 $\mathbf x$ 積分。如果把每個像素切成 10 個格子，一張 100 個像素的小圖就有 $10^{100}$ 種組合，根本加不完。沒有 $Z$，就不知道任何一個點的確切機率，也沒辦法像高斯那樣直接抽樣。

但有一件事永遠算得出來：兩個點的機率比，因為 $Z$ 會消掉：

$$\frac{p(\mathbf x')}{p(\mathbf x)} = \frac{e^{-E(\mathbf x')}/Z}{e^{-E(\mathbf x)}/Z} = e^{-(E(\mathbf x') - E(\mathbf x))}$$

MCMC（Markov chain Monte Carlo）只用這個比值就能取樣。想法是：讓一個漫步者在空間裡走，每一步只比較「下一個位置」和「現在的位置」誰的機率高，並按一個精心設計的規則決定要不要走過去。走得夠久之後，漫步者停留在每個地方的時間比例，就等於目標機率。

要做到這件事，需要回答三個問題，也就是這一頁的主線：

隨機漫步走久了會發生什麼？（第 2–3 節）
怎麼設計規則，讓它走久了之後的分佈剛好是我們要的？（第 4–5 節）
實際上要走多久？什麼時候會失敗？（第 6–8 節）

01需要的機率知識

只需要兩條規則。

條件機率 $P(B\mid A)$ 是「已知 $A$ 發生，$B$ 發生的機率」。兩件事都發生的機率是：

$$P(A\text{ 且 }B) = P(A)\,P(B\mid A)$$

全機率：要算 $B$ 的機率，就把所有可能的情況加起來：

$$P(B) = \sum_iP(A_i)\,P(B\mid A_i)$$

手算

今天晴天的機率是 0.8。晴天之後明天晴的機率是 0.9，雨天之後明天晴的機率是 0.5。明天晴天的機率是：

$$0.8\times0.9 + 0.2\times0.5 = 0.82$$

這個計算就是 Markov 鏈走一步。

02Markov 鏈：下一步只看現在在哪裡

定義

一個 Markov 鏈是一串隨機狀態 $x^{(1)}, x^{(2)}, \dots$，下一步的機率只依賴現在的狀態，不管過去怎麼走到這裡。就像下棋的下一步只看目前的棋盤，不必知道之前怎麼下。

有限個狀態時，把規則排成轉移矩陣 $T$：$T_{ij}$ 是「現在在 $i$，下一步到 $j$」的機率。每一列是「從某個狀態出發的所有去向」，所以每一列加起來是 1。

例子：天氣

晴天之後 90% 還是晴天；雨天之後一半一半。

$$T = \begin{bmatrix}0.9&0.1\\0.5&0.5\end{bmatrix}\quad\begin{matrix}\leftarrow\text{今天晴}\\\leftarrow\text{今天雨}\end{matrix}$$

同一條鏈的圖。每個狀態射出去的箭頭，機率加起來是 1。

走一步：矩陣乘法就是全機率

把現在在每個狀態的機率排成一個橫的向量（列向量）$\boldsymbol\mu_t$。第 1 節的全機率寫成：

$$\mu_{t+1}(j) = \sum_i\mu_t(i)\,T_{ij}$$

這正好是矩陣乘法的定義，所以可以寫成：

$$\boldsymbol\mu_{t+1} = \boldsymbol\mu_t\,T$$

走 $n$ 步就是乘 $n$ 次：$\boldsymbol\mu_{t+n} = \boldsymbol\mu_tT^n$。

手算

第一天確定是晴天，$\boldsymbol\mu_0 = (1, 0)$。

第 1 天：晴的機率 $1\cdot0.9 + 0\cdot0.5 = 0.9$，雨的機率 $1\cdot0.1 + 0\cdot0.5 = 0.1$。
第 2 天：晴 $0.9\cdot0.9 + 0.1\cdot0.5 = 0.86$，雨 $0.9\cdot0.1 + 0.1\cdot0.5 = 0.14$。
第 3 天：$(0.844, 0.156)$
第 4 天：$(0.8376, 0.1624)$

數字越來越接近某個固定的值，而且從雨天出發也一樣：晴的機率依序是 $0, 0.5, 0.7, 0.78, \dots$。

穩態分佈

那個固定的值叫穩態分佈 $\boldsymbol\pi$，它滿足「再走一步也不變」：

穩態 $$\boldsymbol\pi = \boldsymbol\pi\,T$$

手算：解出穩態

設 $\boldsymbol\pi = (a, b)$，$a + b = 1$。第一個座標的方程是 $a = 0.9a + 0.5b$，移項得 $0.1a = 0.5b$，所以 $a = 5b$。代入 $a + b = 1$：$b = \tfrac16,\ a = \tfrac56$。

所以長期來看，83.3% 的日子是晴天。驗算：$\tfrac56\cdot0.9 + \tfrac16\cdot0.5 = 0.8333$ ✓。

注意 $0.1a = 0.5b$ 的意思：每天從晴變雨的機率流量，等於從雨變晴的機率流量。第 4 節會用到這個看法。

時間平均等於穩態

穩態分佈有兩種意思，而它們相同：

很多個平行世界，在同一時刻，有多少比例是晴天。
單一世界，過很多天，有多少比例的日子是晴天。

第二種才是 MCMC 能用的原因：只要讓一條鏈一直走，把走過的位置都記下來，它們就是目標分佈的樣本，可以用來估計期望值。

Demo D1 · 三狀態鏈：分佈的演化與一個漫步者的足跡

怎麼看：這是 Andrieu et al. (2003) 圖 4 的例子，$T = \begin{bmatrix}0&1&0\\0&0.1&0.9\\0.6&0.4&0\end{bmatrix}$。左圖的圓圈大小是第 $t$ 步的機率，箭頭上的數字是轉移機率。右圖是機率隨步數的變化，大約 10 步之後就貼上穩態。切換起始分佈，答案都一樣。把 $n$ 拉大，看讀數最後一行：漫步者停留的時間比例，會越來越接近 $\boldsymbol\pi$。

起始分佈

走了 t 步：t =

一個漫步者走 n 步：n =

x₁x₂x₃虛線：穩態 π

論文寫的穩態是 (0.2, 0.4, 0.4)、第一步是 (0.2, 0.6, 0.2)，都是四捨五入後的數字；精確值見讀數。

什麼時候不會收斂

收斂需要兩個條件（Andrieu et al. §3）：

不可約（irreducible）：每個狀態都能走到每個狀態。反例：$T = \begin{bmatrix}1&0\\0&1\end{bmatrix}$，兩個狀態完全不相通，從哪裡出發就永遠待在哪裡。
非週期（aperiodic）：不會被困在固定的循環裡。反例：$T = \begin{bmatrix}0&1\\1&0\end{bmatrix}$，每天必定換天氣。從晴天出發，機率永遠在 $(1,0)$ 和 $(0,1)$ 之間跳，不會收斂。

記號提醒：這裡沿用 Andrieu et al. 的「列向量 × 矩陣，每列加總為 1」。第 18 講和補充課 G 的速率矩陣用的是「矩陣 × 行向量」，兩者是轉置關係。

一句話帶走：Markov 鏈走一步就是乘上轉移矩陣；在溫和的條件下，不論從哪裡出發都會收斂到同一個穩態，而且時間平均等於穩態平均。

小練習：$T = \begin{bmatrix}0.8&0.2\\0.4&0.6\end{bmatrix}$ 的穩態是什麼？

$0.2a = 0.4b$，所以 $a = 2b$，$\boldsymbol\pi = (\tfrac23, \tfrac13)$。

03收斂有多快：看第二大特徵值

穩態方程 $\boldsymbol\pi T = \boldsymbol\pi$ 就是特徵向量的定義，只是向量在左邊（補充課 A 第 4 節）：$\boldsymbol\pi$ 是 $T$ 的特徵值 1 對應的特徵向量。其他特徵值的絕對值都小於 1。

現在把分佈拆成「穩態」加上「偏離穩態的部分」。每走一步，穩態不變，偏離的部分就乘上特徵值：

$$\boldsymbol\mu_t = \boldsymbol\pi + (\text{偏離})\cdot\lambda_2^t$$

所以第二大的特徵值 $|\lambda_2|$ 決定收斂速度。

手算：天氣的例子

$T$ 的 trace 是 $1.4$，也就是特徵值的和。其中一個是 1，所以 $\lambda_2 = 0.4$。

檢查第 2 節的數字。下雨機率和穩態 $\tfrac16$ 的差距：

差距 $|\mu_t(\text{雨}) - \tfrac16|$ 每天乘上 0.4：$0.1667\to0.0667\to0.0267\to\cdots$，就是 $\lambda_2$。

最簡單的一般例子：兩個狀態，每一步有機率 $s$ 留在原地。$T = \begin{bmatrix}s&1-s\\1-s&s\end{bmatrix}$，trace 是 $2s$，所以 $\lambda_2 = 2s - 1$。$s$ 越接近 1，鏈越「黏」，收斂越慢，而且相鄰的樣本幾乎一樣。

自相關

自相關 $\rho_k$ 是「相隔 $k$ 步的兩個樣本有多像」，定義是相隔 $k$ 步的相關係數。$\rho_k = 1$ 是完全一樣，$\rho_k = 0$ 是毫無關係。獨立抽樣時 $\rho_k = 0$。這條兩狀態鏈的理論值是 $\lambda_2^k$，自相關消失得越慢，同樣多的樣本裡能提供的資訊越少，第 6 節會把它變成一個數字。

Demo D2 · 黏著程度與收斂、自相關

怎麼看：上面的色帶是一條 300 步的路徑。$s$ 接近 0.5 時顏色快速交替，$s$ 接近 1 時一大段都是同一個顏色。左下是和穩態的距離（log 刻度，所以是直線），右下是自相關。把 $s$ 從 0.5 拉到 0.99，看兩張圖怎麼變慢。

留在原地的機率 s =

一條 300 步的路徑（藍 = 狀態 0，橘 = 狀態 1）

左：和穩態的距離（從狀態 0 出發），log 刻度右：自相關（實測）理論 λ₂ᵏ

一句話帶走：偏離穩態的部分每步乘上 $\lambda_2$，所以 $|\lambda_2|$ 越接近 1，收斂越慢，樣本也越相關。

小練習：$s = 0.95$ 時 $\lambda_2$ 是多少？差距大約要幾步才縮成原來的一半？

$\lambda_2 = 0.9$。要 $0.9^t = 0.5$，$t = \ln0.5/\ln0.9\approx6.6$，大約 7 步。

04Detailed balance：每一對狀態之間的流量互相抵銷

為什麼需要它

前兩節是「給定規則，算穩態」。MCMC 要反過來：給定想要的穩態 $\boldsymbol\pi$，設計規則 $T$。直接解 $\boldsymbol\pi T = \boldsymbol\pi$ 很難，因為它把所有狀態綁在一起。我們需要一個容易檢查的條件。

機率流量

把機率想成水。穩態時，從 $i$ 流到 $j$ 的流量是 $\pi_iT_{ij}$（在 $i$ 的機率 × 走到 $j$ 的機率）。$\boldsymbol\pi T = \boldsymbol\pi$ 的意思是每個狀態的總流入等於總流出，水位才不會變，這叫 global balance。

一個更強、也更容易檢查的條件是：每一對狀態之間，來回的流量都相等：

Detailed balance $$\pi_i\,T_{ij} = \pi_j\,T_{ji}\qquad\text{對所有 }i, j$$

左：只要求總量平衡。右：每一對都平衡，總量當然也平衡。

手算

天氣的例子：$\pi_{\text{晴}}T_{\text{晴雨}} = \tfrac56\times0.1 = 0.0833$，$\pi_{\text{雨}}T_{\text{雨晴}} = \tfrac16\times0.5 = 0.0833$。相等 ✓。

為什麼它保證 $\boldsymbol\pi$ 是穩態

把 detailed balance 兩邊對 $i$ 加總：

左邊 $\sum_i\pi_iT_{ij}$，就是「走一步之後在 $j$ 的機率」，也就是 $(\boldsymbol\pi T)_j$。
右邊 $\sum_i\pi_jT_{ji} = \pi_j\sum_iT_{ji}$。$\pi_j$ 和 $i$ 無關，可以提出來。
$\sum_iT_{ji}$ 是「從 $j$ 出發的所有去向」，加起來是 1。所以右邊就是 $\pi_j$。

於是 $(\boldsymbol\pi T)_j = \pi_j$，也就是 $\boldsymbol\pi T = \boldsymbol\pi$。設計 MCMC 時，只要讓規則對目標分佈滿足 detailed balance，目標分佈就自動是穩態。

它只牽涉兩個狀態，不必解整個方程，這就是它好用的地方。滿足 detailed balance 的鏈也叫可逆的（reversible）：倒著播放和正著播放，統計上看不出差別。

Demo D3 · 流量圖：箭頭粗細 = 機率流量 πᵢTᵢⱼ

怎麼看：比較每一對狀態之間的兩支箭頭。Metropolis 鏈（下一節的方法）每一對都一樣粗。繞圈的鏈順時針的箭頭明顯比較粗，有淨的循環流量，但每個狀態的總流入仍然等於總流出，所以 (1/3, 1/3, 1/3) 還是穩態：detailed balance 是充分但不必要的條件。

哪一條鏈

出處：Andrieu et al. (2003) §3 的 detailed balance 條件與加總證明。

一句話帶走：每一對狀態的來回流量相等，總流入就等於總流出，目標分佈自動成為穩態。

小練習：想要 $\boldsymbol\pi = (0.25, 0.75)$，已經決定 $T_{12} = 0.6$。$T_{21}$ 該是多少？

$0.25\times0.6 = 0.75\times T_{21}$，所以 $T_{21} = 0.2$。

05Metropolis–Hastings：先提議，再決定要不要接受

想法

目標：從只知道 $\tilde p(x)\propto p(x)$ 的分佈取樣。直接設計 $T$ 滿足 detailed balance 不容易，所以把一步拆成兩個動作：

提議：隨便用一個容易抽樣的分佈 $q(x'\mid x)$ 提出一個候選點，例如 $x' = x + \sigma\epsilon$。
接受或拒絕：以機率 $A(x, x')$ 移過去，否則留在原地。

提議分佈是隨便選的，所以不會滿足 detailed balance。接受率 $A$ 的任務，就是把差距補回來。

推導接受率

從 $x$ 走到另一個點 $x'$ 的機率是「提議 × 接受」：$T(x\to x') = q(x'\mid x)\,A(x, x')$。
detailed balance 要求 $$\begin{aligned}&p(x)\,q(x'\mid x)\,A(x, x')\\={}&p(x')\,q(x\mid x')\,A(x', x)\end{aligned}$$
把 $A$ 集中到一邊，定義比值 $r$： $$\frac{A(x, x')}{A(x', x)} = \frac{p(x')\,q(x\mid x')}{p(x)\,q(x'\mid x)} = r$$
只要 $A$ 的比值等於 $r$ 就好，選法有很多種。我們希望盡量多接受（否則鏈一直不動），所以讓其中較大的那個等於 1：
- $r\ge1$：取 $A(x, x') = 1$，$A(x', x) = 1/r$。
- $r \lt 1$：取 $A(x, x') = r$，$A(x', x) = 1$。
兩種情況合起來：

Metropolis–Hastings 接受率 $$A(x, x') = \min\Big\{1,\ \frac{\tilde p(x')\,q(x\mid x')}{\tilde p(x)\,q(x'\mid x)}\Big\}$$

比值裡 $Z$ 消掉了，所以用 $\tilde p$ 就夠。「留在原地」的機率不影響 detailed balance，因為 $x = x'$ 時兩邊本來就一樣。

提議分佈對稱時（$q(x'\mid x) = q(x\mid x')$，例如上面的高斯隨機漫步），兩個 $q$ 也消掉：

$$A = \min\Big\{1,\ \frac{\tilde p(x')}{\tilde p(x)}\Big\}$$

這就是第 11 講投影片 24 的 Metropolis 演算法：往機率更高的地方一定接受，往機率低的地方以比值的機率接受。

往低處不是一定拒絕，所以鏈會探索整個分佈，而不是爬到峰頂就停住。低多少，接受機率就小多少。

手算一：連續的一步

目標 $\tilde p(x) = e^{-x^2/2}$，現在 $x = 1$。

提議 $x' = 0.4$：$\tilde p$ 變大，一定接受。
提議 $x' = 2$：比值 $e^{-2}/e^{-0.5} = e^{-1.5} = 0.223$。抽一個 $[0,1]$ 的亂數 $u$，$u \lt 0.223$ 就走過去，否則留在 $x = 1$，並把 1 再記一次。

手算二：三個狀態，算出整個 T

$\boldsymbol\pi = (0.2, 0.3, 0.5)$，提議：等機率選另外兩個狀態之一（每個 $\tfrac12$）。

$T_{12} = \tfrac12\times\min(1, \tfrac{0.3}{0.2}) = 0.5$；$T_{21} = \tfrac12\times\min(1, \tfrac{0.2}{0.3}) = \tfrac13$。
流量：$0.2\times0.5 = 0.1$，$0.3\times\tfrac13 = 0.1$。相等 ✓。
其他對也一樣，留在原地的機率就是剩下的部分。完整的矩陣：

$$T = \begin{bmatrix}0&\tfrac12&\tfrac12\\\tfrac13&\tfrac16&\tfrac12\\\tfrac15&\tfrac3{10}&\tfrac12\end{bmatrix}$$

驗算：$\boldsymbol\pi T$ 的第一個座標是 $0.3\times\tfrac13 + 0.5\times\tfrac15 = 0.2$ ✓。整個過程只用到機率的比值。

Demo D4 · 一步一步跑 Metropolis

怎麼看：這就是上面手算的三狀態鏈。按「走一步」，讀數會寫出這一步的提議、比值和決定。多走幾步，再按「走 1000 步」，看長條（實際停留的比例）怎麼逼近灰色的目標 $(0.2, 0.3, 0.5)$。注意演算法從頭到尾都只用了比值。

停留的比例目標 π目前所在的狀態

連續的例子：步長怎麼選

目標 $\tilde p(x) = 0.3e^{-0.2x^2} + 0.7e^{-0.2(x-10)^2}$：兩個峰在 0 和 10，各自的標準差約 1.58，右邊的峰佔 70%。提議是 $x' = x + \sigma\epsilon$。

Demo D5 · MH 取樣一個雙峰分佈（Andrieu et al. 圖 6–7 的目標）

怎麼看：拖動 $\sigma$（對數刻度），看軌跡和直方圖的三種情況：

$\sigma$ 太小（例如 0.3）：幾乎每次都接受，但每步只挪一點點，鏈困在一個峰裡，根本看不到另一個峰。估出來的 $P(x > 5)$ 完全錯。
$\sigma$ 適中（3 到 10 之間）：兩個峰之間來回跳，直方圖貼著目標。Andrieu et al. 用的是 $\sigma = 10$。
$\sigma$ 太大（例如 80）：提議點大多落在機率極低的地方，接受率很低，鏈常常原地不動，軌跡變成一階一階的水平線段。

提議步長 σ =

起點 x₀ =

丟掉前面幾步（burn-in）=

軌跡與直方圖目標 p(x)灰色區域：被丟掉的 burn-in自相關

出處：Andrieu et al. (2003) §3.1 的 MH 演算法與收斂條件。收斂還需要「提議能到達所有地方」（不可約）和「有機會拒絕」（非週期）。

一句話帶走：先隨便提議，再用 $\min\{1, r\}$ 接受，這個比值剛好補足 detailed balance，而且 $Z$ 消掉了。

小練習：$\tilde p(x) = e^{-x^4}$，現在 $x = 0$，對稱提議 $x' = 1$。接受機率是多少？反過來從 $1$ 提議 $0$ 呢？

$\tilde p(1)/\tilde p(0) = e^{-1}\approx0.368$。反方向比值大於 1，一定接受。

06怎麼知道跑夠了：burn-in、自相關、有效樣本數

Burn-in

起點通常不在高機率區（Demo D5 預設從 $-20$ 出發），前面一段樣本反映的是起點而不是 $p$，應該丟掉。在 Demo D5 把 burn-in 拉到 0，看直方圖左邊多出來的雜訊。

有效樣本數

MCMC 的樣本彼此相關，$N$ 個相關樣本提供的資訊比 $N$ 個獨立樣本少。少多少？我們看 $N$ 個樣本的平均的變異數（補充課 A 第 8 節）。

獨立時：$\operatorname{Var}(\bar x) = \sigma^2/N$。
相關時，交叉項不會消失： $$\operatorname{Var}\Big(\sum_ix_i\Big) = \sum_i\sum_j\operatorname{Cov}(x_i, x_j)$$
相隔 $k$ 步的共變異數是 $\sigma^2\rho_k$，$N$ 很大時，每個 $k$ 都約有 $N$ 對（前後各一次）： $$\begin{aligned}&\sum_i\sum_j\operatorname{Cov}(x_i, x_j)\\\approx{}&N\sigma^2\Big(1 + 2\sum_{k\ge1}\rho_k\Big)\end{aligned}$$
除以 $N^2$，就是平均的變異數。和獨立的 $\sigma^2/N_{\text{eff}}$ 對照：

有效樣本數 $$N_{\text{eff}} = \frac{N}{1 + 2\sum_{k\ge1}\rho_k}$$

手算

Demo D2 的兩狀態鏈：$\rho_k = \lambda^k$。等比級數 $\sum_{k\ge1}\lambda^k = \frac{\lambda}{1-\lambda}$，所以

$$1 + \frac{2\lambda}{1-\lambda} = \frac{1+\lambda}{1-\lambda}$$$$N_{\text{eff}} = N\,\frac{1-\lambda}{1+\lambda}$$

$s = 0.99$ 時 $\lambda = 0.98$，10,000 個樣本只相當於 $10000\times\frac{0.02}{1.98}\approx101$ 個獨立樣本。

沒有保證

實務上會同時跑好幾條從不同起點出發的鏈，看它們是否給出一樣的答案。但沒有一個診斷方法能保證「已經收斂」，只能找出「明顯沒收斂」。

出處：Andrieu et al. (2003) §4（burn-in、mixing time、用特徵值界定收斂速度）。$N_{\text{eff}}$ 的公式是 MCMC 文獻的標準定義；本頁的 demo 用 Geyer 的 initial positive sequence 截斷 $\sum\rho_k$。

一句話帶走：丟掉前面的 burn-in，用有效樣本數衡量真正有多少資訊；診斷只能抓出沒收斂，不能保證已收斂。

07Gibbs sampling：一次只更新一個座標

方法

多維時，提議分佈很難選：步長太大幾乎都被拒絕，太小又走不動。但有時候，每個變數「給定其他變數」的條件分佈很容易取樣。那就輪流更新：

$$\begin{aligned}x_1 &\sim p(x_1\mid x_2)\\x_2 &\sim p(x_2\mid x_1)\end{aligned}$$

更新 $x_1$ 時用的是最新的 $x_2$，反之亦然。重複這兩步。

為什麼這樣是對的：它是 MH 的特例

把「更新 $x_1$」當成提議：$\mathbf x = (x_1, x_2)\to\mathbf x' = (x_1', x_2)$，提議機率 $q = p(x_1'\mid x_2)$。代入接受率，用 $p(x_1, x_2) = p(x_2)\,p(x_1\mid x_2)$ 展開：

$$r = \frac{p(x_2)\,p(x_1'\mid x_2)\cdot p(x_1\mid x_2)}{p(x_2)\,p(x_1\mid x_2)\cdot p(x_1'\mid x_2)} = 1$$

分子分母完全相同，所以永遠接受（Andrieu et al. §3.4）。第 11 講的 Ising model、RBM 都用它（RBM 的 block Gibbs：給定隱藏層、整層可見層一起更新）。

弱點：變數高度相關時很慢

每次只能沿座標軸走，在細長的橢圓裡要走很多小碎步。

灰色是分佈的等高線。沿對角線一步就能走很遠，但 Gibbs 只能橫著和直著走，每一步都被橢圓的窄邊限制住。

推導：慢多少

二維標準高斯，相關係數 $\rho$。由補充課 A 第 11 節，條件分佈是 $x_1\mid x_2\sim\mathcal N(\rho x_2,\ 1-\rho^2)$，反之亦然。所以一輪更新是（$\epsilon, \epsilon'$ 是獨立的雜訊）：

$$x_2^{\text{新}} = \rho x_1 + \sqrt{1-\rho^2}\,\epsilon$$ $$x_1^{\text{新}} = \rho x_2^{\text{新}} + \sqrt{1-\rho^2}\,\epsilon'$$

把第一式代入第二式，$x_1^{\text{新}} = \rho^2x_1 + \text{雜訊}$。所以 $x_1$ 的一步自相關就是 $\rho^2$。$\rho = 0.99$ 時是 $0.98$，由第 6 節，有效樣本比例只有 $\frac{1-\rho^2}{1+\rho^2}\approx1\%$。

Demo D6 · 二維高斯上的 Gibbs：沿座標軸走的鋸齒

怎麼看：藍色是前 30 步的路徑。$\rho = 0$ 時一步就能跳到任何地方。把 $\rho$ 拉到 0.99，路徑變成一小段鋸齒，比對讀數裡的自相關和 $\rho^2$。

相關係數 ρ =

前 30 步的鋸齒路徑3,000 個樣本真實分佈的 1σ、2σ 橢圓

一句話帶走：Gibbs 輪流從條件分佈抽樣，永遠接受；但變數高度相關時，只能走小碎步。

08多峰與溫度：為什麼鏈會被困住，怎麼解

為什麼會被困住

如果兩個峰之間隔著一大片機率極低的區域，局部的提議（小步長的 MH、Gibbs，以及補充課 F 的 Langevin）要跨過去，就必須連續接受很多次「往低機率走」的提議。

手算

能量 $E = -\log\tilde p$。如果低谷比峰頂高 $\Delta E = 20$，一口氣跳到低谷的接受機率是 $e^{-20}\approx2\times10^{-9}$，大約要試五億次。鏈會在一個峰裡待很久，估出來的峰的比例完全錯。

這正是第 13 講「三個陷阱」裡的第三個：score（或局部的 MCMC）看不出遠方另一個峰有多大。

加熱

一個辦法是改從 $p^{1/T}$（溫度 $T > 1$）取樣。因為 $p\propto e^{-E}$，所以 $p^{1/T}\propto e^{-E/T}$：所有能量差都除以 $T$。

手算

同樣的低谷，$T = 10$ 時能量差變成 2，接受機率變成 $e^{-2}\approx0.135$，很快就能跨過去。

代價是分佈也變了。高斯的 $p^{1/T}$ 是 $e^{-(x-\mu)^2/(2T\sigma^2)}$，也就是變異數變成 $T\sigma^2$，標準差變成 $\sqrt T$ 倍。

提高溫度讓峰變矮變寬，低谷被填高。（兩條曲線都把峰頂縮放到相同高度，方便比較形狀。）

加熱之後的分佈已經不是 $p$。實務上的做法是同時跑好幾個溫度的鏈、偶爾交換狀態（parallel tempering），或從高溫開始逐漸降溫（annealing）。第 14 講的 annealed Langevin 先在大雜訊（類似高溫）下跑，再逐步減少雜訊，是同一種想法。

Demo D7 · 兩個分得很開的峰：溫度越高越容易跨，但分佈也跟著變形

怎麼看：$T = 1$ 時，鏈幾乎只待在出發的那個峰，右上的峰明明佔 50%，卻幾乎沒有樣本。把 $T$ 調高，鏈開始在兩峰之間來回（看讀數的跨越次數），但每個峰都被拉寬，樣本不再代表原來的分佈。

溫度 T =

p^{1/T} 的等高線（越深越高）MH 樣本（20,000 步，從左下的峰出發）

一句話帶走：能量差決定跨越有多難；加熱把能量差除以 $T$，容易跨越，但分佈也變了，所以要逐步降溫或交換。

09在 CS236 哪裡用到

本頁的概念	出現的地方
不需要 $Z$ 的取樣	第 11 講：從 EBM 取樣、contrastive divergence 要的負樣本
Metropolis、接受率	第 11 講投影片 24 與 Demo E；第 12 講開頭「Metropolis MCMC 為什麼會收斂」
Detailed balance	第 12 講用它說明 Metropolis 與 Langevin 的穩態；補充課 F 的 MALA
Gibbs、block Gibbs	第 11 講的 Ising 去雜訊、RBM
混合慢、多峰陷阱	第 13 講的三個陷阱、第 14 講的 annealed Langevin
短鏈、沒跑到穩態	第 11 講 contrastive divergence 只跑幾步 MCMC 的偏差

常見誤解

誤解：「被拒絕的提議就不算樣本」。拒絕時要把目前的 $x$ 再記一次。只記錄被接受的點，會讓分佈偏向容易被接受的區域。

誤解：「接受率越高越好」。$\sigma$ 很小時接受率接近 100%，但鏈幾乎不動（Demo D5）。好的步長在「接受率」和「每步走多遠」之間取捨。

誤解：「穩態一定要滿足 detailed balance」。那只是充分條件。不可逆的鏈也可以有穩態（Demo D3 的繞圈鏈），甚至常常混合得更快。

誤解：「直方圖看起來很平滑，鏈就收斂了」。鏈可能困在一個峰裡，畫出一個漂亮但錯誤的直方圖（Demo D5 的小 $\sigma$、Demo D7 的 $T = 1$）。

誤解：「10,000 個 MCMC 樣本 = 10,000 個獨立樣本」。要用有效樣本數衡量；黏的鏈可能只相當於一百個。

自我測驗

1. 兩狀態鏈 $T = \begin{bmatrix}0.9&0.1\\0.3&0.7\end{bmatrix}$ 的穩態分佈是什麼？

解 $\pi_1\cdot0.1 = \pi_2\cdot0.3$（兩狀態時 global balance 就是 detailed balance），得 $\boldsymbol\pi = (0.75, 0.25)$。

2. 寫出 detailed balance，並證明它推得出穩態。

$\pi_iT_{ij} = \pi_jT_{ji}$；對 $i$ 加總：$\sum_i\pi_iT_{ij} = \pi_j\sum_iT_{ji} = \pi_j$，即 $\boldsymbol\pi T = \boldsymbol\pi$。

3. 對稱提議下的 Metropolis 接受率是什麼？為什麼不需要 $Z$？

$\min\{1, \tilde p(x')/\tilde p(x)\}$；$Z$ 同時出現在分子分母，消掉了。

4. 兩狀態「黏」鏈的 $\lambda_2 = 0.9$。10,000 個樣本的有效樣本數大約多少？

$10000\cdot\frac{1-0.9}{1+0.9}\approx526$。

5. Gibbs 在 $\rho = 0.99$ 的二維高斯上，$x_1$ 的一步自相關是多少？為什麼這麼慢？

$\rho^2\approx0.98$。每步只能沿座標軸移動，在很細長的橢圓裡要走很多小碎步才能從一端到另一端。

6. 為什麼從 $p^{1/T}$（$T > 1$）取樣比較容易在峰之間移動？代價是什麼？

能量差都除以 $T$，跨越低谷時的接受機率變大。代價是取到的已經不是 $p$，需要 tempering 或 annealing 把結果拉回原本的分佈。

本頁的概念	出現的地方
不需要 \(Z\) 的取樣	第 11 講：從 EBM 取樣、contrastive divergence 要的負樣本
Metropolis、接受率	第 11 講投影片 24 與 Demo E；第 12 講開頭「Metropolis MCMC 為什麼會收斂」
Detailed balance	第 12 講用它說明 Metropolis 與 Langevin 的穩態；補充課 F 的 MALA
Gibbs、block Gibbs	第 11 講的 Ising 去雜訊、RBM
混合慢、多峰陷阱	第 13 講的三個陷阱、第 14 講的 annealed Langevin
短鏈、沒跑到穩態	第 11 講 contrastive divergence 只跑幾步 MCMC 的偏差