☰

圓州率
🌐
☰

Feature Image

Markov Chain Monte Carlo (MCMC)

數學, 機率論, 馬可夫鏈 (Markov Chain)
介紹 Markov Chain Monte Carlo (MCMC)，一種適用於複雜分佈上採樣的模型。

介紹

若要對分佈採樣或計算期望值問題，可以透過 accept-reject method 等 Monte Carlo 方法；但對於複雜、多變量、變數不獨立的分佈，可以透過 Markov Chain Monte Carlo (MCMC) 方法。

假設一隨機變量 $x \in X$ 帶有 pdf $p (x)$，如何對 $p (x)$ 採樣？MCMC 方法是找出一個 Markov chain $X = \{ X_0, X_1, \cdots, X_t, \cdots \}$ 定義在狀態空間 $S$ 上且滿足 ergodic theorem，則 $X$ 的穩定分佈就會與採樣分佈 $p$ 相同，接著在此 Markov chain 隨機採樣，當時間趨於無限大時，採樣就會收斂至穩定分佈。

因此核心問題就是如何建構此 Markov chain，一種方案是找出 transition function 並建構 reversible Markov chain 使其滿足 ergodic theorem，因此起始點就會與結果無關，最終都會收斂至唯一的穩定分佈，一般常用的方案有 Metropolis–Hastings algorithm 或 Gibbs sampling。

MCMC 在相鄰兩次的採樣是相關的，若需要獨立樣本可以考慮每隔數個週期後再接受結果，並且相對於 accept-reject method，並不需要找到已知且容易採樣的分佈進行拒絕，因為在高維度不容易尋找，一般而言 MCMC 方法的效率會高於 accept-reject method。

步驟

  1. 建構滿足 ergodic theorem 的 Markov chain，並使其穩定分佈與採樣目標 $p (x)$ 相同
  2. 選定狀態空間中任意起始點 $x_0$，並透過 Markov chain 執行隨機漫步以獲取 $\{ x_0, x_1, \cdots, x_t, \cdots \}$
  3. 選定兩整數 $m < n$，選取 $\{ x_{m + 1}, x_{m + 2}, \cdots, x_n \}$ 作為採樣結果

然而有幾個點需要注意

  1. 如何定義 Markov chain 使其滿足 ergodic theorem
  2. 如何選取收斂步數 $m$ 使結果不偏
  3. 如何選取迭代步數 $n$ 使結果精確

MCMC 與機器學習

MCMC 能用於機器學習中的模型學習與推理，使其在機器學習中扮演重要角色，尤其在 Bayesian learning；考慮觀測到的資料為隨機變量 $y \in Y$ 而模型參數為 $x \in X$，Bayesian learning 會選取能使得 posterior 機率最大的參數，其中 posterior

$$\begin{align*} p (x | y) = \frac{p (x) p (y | x)}{\int_X P (y | x') p (x') dx'} \end{align*}$$

在其中有 3 種能應用 MCMC 的場景

  1. Normalization:

    $$\begin{align*} \int_X p (y | x') p (x') dx' \end{align*}$$

  2. Marginalization: 若存在隱藏變數 $z \in Z$，則需計算

    $$\begin{align*} p (x | y) = \int_Z p (x, z | y) dz \end{align*}$$

  3. Expectation: 若需計算 $f (x)$ 在分佈上的期望值

    $$\begin{align*} E_p [f (X)] = \int_X f (x) p (x | y) dx \end{align*}$$

上述任一情況如果模型複雜，則計算將會非常困難，但能用 MCMC 做一個通用且有效率的選項

Metropolis–Hasting Algorithm

這邊介紹 MCMC 方法中常用的一個方案 Metropolis–Hasting algorithm。給定 pdf $p (x)$ 作為採樣目標，Metropolis–Hasting algorithm 建構一個 Markov chain 與 transition function $p (x, x')$

$$\begin{align*} p (x, x') = q (x, x') \alpha (x, x') \end{align*}$$

其中 $q (x, x')$ 稱為 proposal distribution 提出狀態改變的可能，$\alpha (x, x')$ 稱為 acceptance distribution，負責接受是否維持當前狀態或是改變狀態。Proposal distribution $q (x, x')$ 是其他 irreducible Markov chain 的 transition function，一般會選擇容易採樣的分佈；acceptance distribution $\alpha (x, x')$ 寫作

$$\begin{align*} \alpha (x, x') = \min \left\lbrace 1, \frac{p (x') q (x', x)}{p (x) q (x, x')} \right\rbrace \end{align*}$$

使得

$$\begin{align*} p (x, x') = \begin{cases} q (x, x') & \text{if } p (x') q(x', x) \geq p (x) q (x, x') \\ q (x, x') \frac{p (x')}{p (x)} & \text{if } p (x') q(x', x) > p (x) q (x, x') \end{cases} \end{align*}$$

而在 Markov chain $p (x, x')$ 上的隨機漫步，是依照下列方式：在時間 $t - 1$ 時為狀態 $x$，設 $x_{t - 1} = x$，接著根據 $q (x, x')$ 生成出下一時間點候選的狀態 $x'$，接著有 $\alpha (x, x')$ 的機率接受此提案，若通過則設 $x_t = x'$；反之則拒絕此提案，維持原狀態 $x_t = x$。

透過此方式建構的 $p (x, x')$ 是 reversible Markov chain，即滿足 ergodic theorem，使得 $p (x, x')$ 的穩定分佈即為採樣目標 $p (x)$。

演算法

在 Markov chain $p (x, x')$ 上的隨機漫步 Metropolis–Hasting algorithm，是依照下列方式

  1. 選定起始點 $x_0$
  2. 對 $i = \conti{n}$ 執行
    1. 設定狀態 $x_{i - 1} = x$，依據 $q (x, x')$ 生成出下一時間點候選的狀態 $x'$
    2. 計算接受機率 $\alpha (x, x')$
    3. 採樣 $u \sim \unif (0, 1)$，若 $u \leq \alpha (x, x')$ 設 $x_i = x'$；反之設 $x_i = x$
  3. 獲取採樣集合 $\{ x_{m + 1}, x_{m + 2}, \cdots, x_n \}$

Proposal Distribution 選取

Proposal Distribution $q (x, x')$ 一般選取容易採樣的分佈，且建議選取兩種類型的

  1. 對稱 (symmetic) 的，即 $q (x, x') = q (x', x)$
  2. 獨立 (independent) 的，即 $q (x, x') = q (x')$

若選取對稱分佈，即 $q (x, x') = q (x', x)$，則 acceptance distribution 能被簡化成

$$\begin{align*} \alpha (x, x') = \left\lbrace 1, \frac{p (x')}{p (x)} \right\rbrace \end{align*}$$

此選取稱為 Metropolis selection；在這中有個特例，設定 $q (x, x') = q ( |x - x'|)$，這種設定被稱為 random walk Metropolis algorithm，這種選取傾向於選擇靠近 $x$ 的 $x'$，例如

$$\begin{align*} q (x, x') \propto \exp \left[ - \frac{(x - x')^2}{2} \right] \end{align*}$$

若選擇獨立分佈，即 $q (x, x') = q (x')$，則

$$\begin{align*} \alpha (x, x') = \left\lbrace 1, \frac{p (x') / q(x')}{p (x) / q (x)} \right\rbrace \end{align*}$$

這種抽樣方式簡單易實現，但收斂速度較慢，且一般 $q (x)$ 會選取較靠近目標分佈 $p (x)$

Full Conditional Distribution

若採樣目標為 $k$ 維度變量，$x = (\contia{x}{k})^T$，若選 $I \subset K = \{ \conti{k} \}$，則 full conditional distribution 可以被寫為 $p (x_{I} | x_{-I})$，其中 $x_{I} = \{ x_i | i \in I \}$，$x_{-I} = \{ x_i | i \not\in I \}$。

Full conditional distribution 擁有以下性質

$$\begin{align*} p (x_I | x_{-I}) = \frac{p (x)}{\int p (x) dx_I} \propto p (x) \end{align*}$$

而任意兩 $x, x'$ 的 Full conditional distribution 擁有以下性質

$$\begin{align*} \frac{p (x_I' | x_{-I}') }{p (x_I | x_{-I}) } = \frac{p (x')}{p (x)} \end{align*}$$

這項性質被用於改進 Metropolis–Hasting Algorithm，其中計算 $\frac{p (x')}{p (x)}$ 會以 $\frac{p (x_I' | x_{-I}') }{p (x_I | x_{-I}) }$ 取代。

Single-Component Metropolis–Hastings Algorithm

在 Metropolis–Hastings Algorithm 需對多變量分佈採樣，有時這種採樣非常困難，因此轉而對單一變數的條件分佈採樣，並遍歷所有變數以此簡化問題，這種方式被稱為 single-component Metropolis–Hastings algorithm

設採樣目標為 $k$ 維度分佈 $x = (\contia{x}{k})^T$，其中 $x_j$ 是 $x$ 的第 $j$ 個變數。設 Markov chain 在時間 $i$ 採樣結果為

$$\begin{align*} x^{(i)} = (x_1^{(i)}, x_2^{(i)}, \cdots, x_k^{(i)})^T \end{align*}$$

其中 $x_{j}^{(i)}$ 為時間 $i$ 下採樣 $x$ 的第 $j$ 個變數。

為了獲取採樣結果 $\{ x^{(1)}, x^{(2)}, \cdots, x^{(n)} \}$，在第 $i - 1$ 迭代結束時，第 $j$ 個變數為 $x_j^{(i - 1)}$，接著在第 $i$ 次迭代的第 $j$ 步會生成 $x_j'^{(i)}$，其生成採樣自 proposed distribution $q (x_j^{(i - 1)}, x_j'^{(i)} | x_{-j}^{(i)})$，其中 $x_{-j}^{(i)}$ 為移除第 $j$ 項的 $x^{(i)}$，且前 $j - 1$ 項為更新後結果，即

$$\begin{align*} x_{-j}^{(i)} = (x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_{k}^{(i - 1)})^T \end{align*}$$

接著決定是否接受此次更新，即

$$\begin{align*} \alpha (x_{j}^{(i - 1)}, x_{j}'^{(i)} | x_{-j}^{(i)} ) = \min \left\lbrace 1, \frac{p (x_j'^{(i)} | x_{-j}^{(i)}) q (x_j'^{(i)}, x_j^{(i - 1)} | x_{-j}^{(i)})}{p (x_j^{(i - 1)} | x_{-j}^{(i)}) q (x_j^{(i - 1)}, x_j^{(i)} | x_{-j}^{(i)})} \right\rbrace \end{align*}$$

若接受此採樣，則 $x_j^{(i)} = x_j'^{(i)}$；反之則拒絕此採樣，維持 $x_j^{(i)} = x_{j}^{(i - 1)}$

Gibbs Sampling

Gibbs sampling 是一種 MCMC 常用的方式，是 single-component Metropolis–Hastings algorithm 的一種特例，但更容易實現，常用於估計多變量參數。其概念是從 joint probability distribution 定義 full conditional distribution，並在其中採樣。

給定 $k$ 維度多變量分佈的 pdf $p (\contia{x}{k})$，從起始點 $x^{(0)} = (x_1^{(0)}, x_2^{(0)}, \cdots, x_k^{(0)})^T$，並在第 $i$ 次迭代的第 $j$ 步驟採樣 $x_j^{(i)}$ 自 $p (x_j | x_{-j}^{(i)})$，其中

$$\begin{align*} p (x_j | x_{- j}^{(i)}) = p (x_j | x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_{k}^{(i - 1)}) \end{align*}$$

其中定義 proposed distribution 是 full conditional distribution，即

$$\begin{align*} q (x, x') = p (x_j' | x_{-j}) \end{align*}$$

透過這個轉換使得 $\alpha (x, x') = 1$，即接受所有採樣，這也是 Gibbs sampling 與 single-component Metropolis–Hastings algorithm 最大的差異，Gibbs sampling 接受所有的採樣，使得每次採樣會保持移動，而非拒絕採樣使得鄰近採樣不會移動。

Gibbs sampling 適用於 full conditional distribution 容易採樣，因此能提高採樣效率；single-component Metropolis–Hastings algorithm 即便 full conditional distribution 不易採樣也適用，但效率較低。

演算法

  • 目標：給定 pdf $p (x)$ 與函數 $f (x)$，回傳採樣結果 $x_{m + 1}, x_{m + 2}, \cdots, x_n$ 與樣本平均 $f_{mn}$
  • 輸入：pdf $p (x)$、函數 $f (x)$、收斂次數 $m$、迭代次數 $n$
  • 輸出：採樣結果 $x_{m + 1}, x_{m + 2}, \cdots, x_n$ 與樣本平均 $f_{mn}$

  1. 設定起始點 $x^{(0)} = (x_1^{(0)}, ,x_2^{(0)}, \cdots, x_k^{(0)})^T$

  2. 設 $i = \conti{n}$

    1. 採樣 $x_1^{(i)}$ 自 $p (x_1 | x_2^{(i - 1)}, x_3^{(i - 1)}, \cdots, x_k^{(i - 1)})$
    2. …
    3. 採樣 $x_j^{(i)}$ 自 $p (x_j | x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_k^{(i - 1)})$
    4. …
    5. 採樣 $x_k^{(i)}$ 自 $p (x_k | x_1^{(i)}, x_2^{(i)}, \cdots, x_{k - 1}^{(i)})$

  3. 獲取採樣集合 $\{ x^{(m + 1)}, x^{(m + 2)}, \cdots, x^{(n)} \}$

  4. 計算樣本平均

    $$\begin{align*} f_{mn} = \frac{1}{n - m} \sum_{i = m + 1}^{n} f (x^{(i)}) \end{align*}$$

範例

用 Gibbs sampling 對二元常態分佈取樣，設

$$\begin{align*} X = (X_1, X_2)^T \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix} \right) \end{align*}$$

則 full conditional distribution 為

$$\begin{align*} X_1 | X_2 & \sim N (\rho X_1, (1 - \rho^2)) \\ X_2 | X_1 & \sim N (\rho X_2, (1 - \rho^2)) \end{align*}$$

接著可透過 Gibbs sampling 採樣，設起始狀態 $x^{(0)} = (x_1^{(0)}, x_2^{(0)})^T$，接著迭代

  1. 第 1 次迭代:
    1. $x_1^{(1)} \sim N (\rho x_2^{(0)}, (1 - \rho^2))$
    2. $x_2^{(1)} \sim N (\rho x_1^{(1)}, (1 - \rho^2))$
  2. …
  3. 第 $i$ 次迭代:
    1. $x_1^{(i)} \sim N (\rho x_2^{(i - 1)}, (1 - \rho^2))$
    2. $x_2^{(i)} \sim N (\rho x_1^{(i)}, (1 - \rho^2))$
  4. …
  5. 第 $i$ 次迭代:
    1. $x_1^{(n)} \sim N (\rho x_2^{(n - 1)}, (1 - \rho^2))$
    2. $x_2^{(n)} \sim N (\rho x_1^{(n)}, (1 - \rho^2))$

選定收斂次數 $m$ 與迭代次數 $n$ 後，$\{ x^{(m + 1)}, x^{(m + 2)}, \cdots, x^{(n)} \}$ 即為採樣結果

採樣計算

在 Gibbs sampling 須對 full conditional distribution 進行多次採樣，因此能透過分佈的特徵來改進或加速採樣效率。以 Bayesian learning 作為範例，設 $y$ 為觀測資料，$x = (\alpha, \theta, z)$ 為參數，其中 $\alpha$ 為超參數 (hyperparameters)，$\theta$ 為模型參數，$z$ 是未觀測資料。Bayesian learning 目的在於最大化 posterior 機率，如下

$$\begin{align*} p (x | y) = p (\alpha, \theta, z | y) \propto p (z, y | \theta) p (\theta | \alpha) p (\alpha) \end{align*}$$

其中 $p (z, y | \theta)$ 是完整資料的分佈，$p (\theta | \alpha)$ 是 prior 分佈，$p (\alpha)$ 是 hyperparameter 分佈。而 Gibbs sampling 估計 $p (x | y)$ 擁有以下特徵

$$\begin{align*} p (\alpha_i | \alpha_{-i}, \theta, z, y) & \propto p (\theta | \alpha) p (\alpha) \\ p (\theta_j | \theta_{-j}, \alpha, \theta, z, y) & \propto p (z, y | \theta) p (\theta | \alpha) \\ p (z_k | z_{-k}, \alpha, \theta, y) & \propto p (z, y | \theta) \end{align*}$$

full conditional probability 是多個 conditional probability 相乘，每個機率分佈僅與少數變數相關，因此對 full conditional probability 可以僅透過數個 conditional probability 的採樣來執行，以此降低運算複雜度

參考資料

  1. LI, Hang. Machine Learning Methods. Springer Nature, 2023.
  2. Lange, Kenneth, J. Chambers, and W. Eddy. Numerical analysis for statisticians. Vol. 1. New York: springer, 2010.