一、概述

1. 概述

在對一個概率分佈進行隨機抽樣，或者是求函式關於該概率分佈的數學期望時可以使用馬爾可夫鏈蒙特卡羅法（MCMC）。相比與拒絕取樣法和重要性取樣法，MCMC更適用於隨機變數是多元的、概率密度函式是非標準形式的、隨機變數各分量不獨立等情況。

對於多元隨機變數$x$，滿足$x\in \mathcal{X}$，其概率密度函式為$p(x)$，$f(x)$為定義在$x\in \mathcal{X}$的函式，目標是獲得概率分佈$p(x)$的樣本集合以及求函式$f(x)$的數學期望$E_{p(x)}[f(x)]$。

應用MCMC解決這個問題。基本想法是：在隨機變數$x$的狀態空間$S$上定義一個滿足遍歷定理的馬爾可夫鏈$X$，使其平穩分佈就是抽樣的目標分佈$p(x)$。然後在這個馬爾可夫鏈上隨機遊走，每個時刻得到一個樣本。

根據遍歷定理，當時間趨於無窮時，樣本的分佈趨於平穩分佈，樣本的函式均值趨近函式的數學期望。所以，當時間足夠長時（時刻大於某個正整數$m$），在之後的時間（時間小於等於某個正整數$n$，$n>m$）裡隨機遊走得到的樣本集合 { x m + 1 , x m + 2 , ⋯   , x n } \left \{x_{m+1},x_{m+2},\cdots ,x_{n}\right \} {xm+1​,xm+2​,⋯,xn​}就是目標概率分佈的抽樣結果，得到的函式均值（遍歷均值）就是要計算的數學期望值：

$$\hat{E}f=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

到時刻$m$為止的時間段稱為燃燒期。

2. 需要注意的幾個知識點

• 由於這個馬爾可夫鏈滿足遍歷定理，隨機遊走的初始點並不影響得到的結果，也就是說從不同的起始點出發，都會收斂到同一平穩分佈。

• MCMC的收斂性的判斷往往是經驗性的，比如，在馬爾可夫鏈上進行隨機遊走，檢驗遍歷均值是否收斂。具體的方法有：
  ①每隔一段時間取一次樣本，得到多個樣本以後，計算遍歷均值，當計算的均值穩定後，認為馬爾可夫鏈已經收斂。
  ②在馬爾可夫鏈上並行進行多個隨機遊走，比較各個隨機遊走的遍歷均值是否接近一致。

• MCMC中得到的樣本序列，相鄰的樣本點是相關的，而不是獨立的。因此，在需要獨立樣本時，可以在該樣本序列中再次進行隨機抽樣，比如每隔一段時間取一次樣本，將這樣得到的子樣本集合作為獨立樣本集合。

• 一般來說，MCMC比拒絕取樣法更容易實現，因為只需要定義馬爾可夫鏈，而不需要定義建議分佈。一般來說MCMC比拒絕取樣效率更高，因為沒有大量被拒絕的樣本，雖然燃燒期的成本也要拋棄。

3. 馬爾可夫鏈蒙特卡羅法的基本步驟

①首先，在隨機變數$x$的狀態空間上構造一個滿足遍歷定義的馬爾可夫鏈，使其平穩分佈為目標分佈$p(x)$;
②從狀態空間的某一點$x_0$出發，用構造的馬爾可夫鏈進行隨機遊走，產生樣本序列$x_0,x_1,\cdots ,x_t,\cdots$;
③應用馬爾可夫鏈的遍歷定理，確定正整數$m$和$n$（$m<n$），得到樣本集合 { x m + 1 , x m + 2 , ⋯   , x n } \left \{x_{m+1},x_{m+2},\cdots ,x_{n}\right \} {xm+1​,xm+2​,⋯,xn​}，求得$f(x)$的均值（遍歷均值）：

$$\hat{E}f=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

這裡有幾個重要問題：

①如何定義馬爾可夫鏈，保證MCMC的條件成立；
②如何確定收斂步數$m$，保證樣本抽樣的無偏性；
③如何確定迭代步數$n$，保證遍歷均值計算的精度。

二、Metropilis-Hastings演算法（MH演算法）

1. 基本原理

假設要抽樣的概率分佈為$p(x)$。MH演算法採用轉移核為$p(x,x^{{}^{\prime }})$的馬爾可夫鏈：

$$p(x,x^{{}^{\prime }})=q(x,x^{{}^{\prime }})\alpha (x,x^{{}^{\prime }})$$

其中$q(x,x^{{}^{\prime }})$稱為建議分佈（proposal distribution），$\alpha (x,x^{{}^{\prime }})$稱為接受分佈（acceptance distribution）。

$q(x,x^{{}^{\prime }})$是另一個馬爾可夫鏈的轉移核，並且$q(x,x^{{}^{\prime }})$是不可約的，即其概率值恆不為$0$，同時也是一個容易抽樣的分佈。接受分佈$\alpha (x,x^{{}^{\prime }})$是：

$$\alpha (x,x^{{}^{\prime }})=min\left\{1,\frac{p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)}{p(x)q(x,x^{{}^{\prime }})}\right\}$$

這時，轉移核$p(x,x^{{}^{\prime }})$可以寫成：

$$p(x,x^{{}^{\prime }})=\{\begin{array}{c}q(x,x^{{}^{\prime }}),p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)\ge p(x)q(x,x^{{}^{\prime }})\\ q(x^{{}^{\prime }},x)\frac{p(x^{{}^{\prime }})}{p(x)},p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)<p(x)q(x,x^{{}^{\prime }})\end{array}$$

轉移核為$p(x,x^{{}^{\prime }})$的馬爾可夫鏈上的隨機遊走以以下方式進行。如果在時刻$t-1$處於狀態$x$，即$x_{t-1}=x$，則先按建議分佈$q(x,x^{{}^{\prime }})$抽樣產生一個候選狀態$x^{{}^{\prime }}$，然後按照接受分佈$\alpha (x,x^{{}^{\prime }})$抽樣決定是否接受狀態$x^{{}^{\prime }}$。以概率$\alpha (x,x^{{}^{\prime }})$接受$x^{{}^{\prime }}$，決定時刻$t$轉移到狀態$x^{{}^{\prime }}$，而以概率$1-\alpha (x,x^{{}^{\prime }})$拒絕$x^{{}^{\prime }}$，決定時刻$t$仍停留在狀態$x$。具體地，從區間$(0,1)$上的均勻分佈中抽取一個隨機數$u$，決定時刻$t$的狀態：

$$x_t=\{\begin{array}{c}{x}^{{}^{\prime }},u\le \alpha (x,x^{{}^{\prime }})\\ x,u>\alpha (x,x^{{}^{\prime }})\end{array}$$

可以證明，轉移核為$p(x,x^{{}^{\prime }})$的馬爾可夫鏈是可逆馬爾可夫鏈（滿足遍歷定理），其平穩分佈就是$p(x)$，即要抽樣的目標分佈。也就是說這是MCMC的一個具體實現。

2. 定理

由轉移核$p(x,x^{{}^{\prime }})=q(x,x^{{}^{\prime }})\alpha (x,x^{{}^{\prime }})$構成的馬爾可夫鏈是可逆的，即
$$p(x)p(x,x^{{}^{\prime }})=p(x^{{}^{\prime }})p(x^{{}^{\prime }},x)$$
並且$p(x)$是該馬爾可夫鏈的平穩分佈。

證明如下：

若$x=x^{{}^{\prime }}$，則上式顯然成立。
若$x\ne x^{{}^{\prime }}$，則：
$$\begin{gathered} p(x)p(x,x^{{}^{\prime }})=p(x)q(x,x^{{}^{\prime }})min\left\{1,\frac{p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)}{p(x)q(x,x^{{}^{\prime }})}\right\} \\ =min\left\{p(x)q(x,x^{{}^{\prime }}),p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)\right\} \\ =p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)min\left\{\frac{p(x)q(x,x^{{}^{\prime }})}{p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)},1\right\} \\ =p(x^{{}^{\prime }})p(x^{{}^{\prime }},x) \end{gathered}$$
由$p(x)p(x,x^{{}^{\prime }})=p(x^{{}^{\prime }})p(x^{{}^{\prime }},x)$知：
$$\begin{gathered} \int p(x)p(x,x^{{}^{\prime }})\mathrm{d}x \\ =\int p(x^{{}^{\prime }})p(x^{{}^{\prime }},x)\mathrm{d}x \\ =p(x^{{}^{\prime }})\int p(x^{{}^{\prime }},x)\mathrm{d}x \\ =p(x^{{}^{\prime }}) \end{gathered}$$
所以$p(x)$是該馬爾可夫鏈的平穩分佈。

3. 建議分佈

建議分佈$q(x,x^{{}^{\prime }})$有多種可能的形式，這裡介紹兩種常用形式。

• 第一種形式

假設建議分佈是對稱的，即對任意的$x$和$x^{{}^{\prime }}$有：

$$q(x,x^{{}^{\prime }})=q(x^{{}^{\prime }},x)$$

這樣的建議分佈稱為Metropolis選擇，也是MH演算法最初採用的建議分佈。這時，接受分佈$\alpha (x,x^{{}^{\prime }})$簡化為：

$$\alpha (x,x^{{}^{\prime }})=min\left\{1,\frac{p(x^{{}^{\prime }})}{p(x)}\right\}$$

Metropolis特例：

①$q(x,x^{{}^{\prime }})=p(x^{{}^{\prime }}|x)$，定義為多元正態分佈，其均值是$x$，其協方差矩陣是常數矩陣（因為協方差矩陣是常數矩陣，所以$q(x,x^{{}^{\prime }})$對稱）。
②$q(x,x^{{}^{\prime }})=q(|x=x^{{}^{\prime }}|)$，這時演算法稱為隨機遊走Metropolis演算法。例如：
$$q(x,x^{{}^{\prime }})\propto exp\left\{-\frac{(x^{{}^{\prime }}-x{)}^2}{2}\right\}$$

Metropolis選擇的特點是當$x^{{}^{\prime }}$與$x$接近時，$q(x,x^{{}^{\prime }})$的概率值高，否則$q(x,x^{{}^{\prime }})$的概率值低。狀態轉移在附近點的可能性更大。

• 第二種形式

第二種形式稱為獨立抽樣。假設$q(x,x^{{}^{\prime }})$與當前狀態$x$無關，即$q(x,x^{{}^{\prime }})=q(x^{{}^{\prime }})$。建議分佈的計算按照$q(x^{{}^{\prime }})$獨立抽樣進行。此時，接受分佈可以寫成：

$$\alpha (x,x^{{}^{\prime }})=min\left\{1,\frac{w(x^{{}^{\prime }})}{w(x)}\right\}$$

其中$w(x^{{}^{\prime }})=\frac{p(x^{{}^{\prime }})}{q(x^{{}^{\prime }})},w(x)=\frac{p(x)}{q(x)}$。

獨立抽樣實現簡單，但可能收斂速度慢，通常選擇接近目標狀態分佈$p(x)$的分佈作為建議分佈$q(x)$。

4. 滿條件分佈

MCMC的目標分佈通常是多元聯合概率分佈$p(x)=p(x_1,x_2,\cdots ,x_k)$，其中$x=(x_1,x_2,\cdots ,x_k{)}^T$為$k$維隨機變數。如果條件概率分佈$p(x_I|x_{-I})$中所有$k$個變數全部出現，其中 x I = { x i , i ∈ I } , x − I = { x i , i ∉ I } , I ⊂ K = { 1 , 2 , ⋯   , k } x_{I}=\left \{x_{i},i\in I\right \},x_{-I}=\left \{x_{i},i\notin I\right \},I\subset K=\left \{1,2,\cdots ,k\right \} xI​={xi​,i∈I},x−I​={xi​,i∈/​I},I⊂K={1,2,⋯,k}，那麼稱這種條件概率分佈為滿條件分佈（full conditional distribution）。

滿條件分佈有以下性質：對任意的$x,x^{{}^{\prime }}\in \mathcal{X}$和任意的$I\subset K$有：

$$p(x_I|x_{-I})=\frac{p(x)}{\int p(x)\mathrm{d}{x}_I}\propto p(x)$$

而且，對任意的$x,x^{{}^{\prime }}\in \mathcal{X}$和任意的$I\subset K$有：

$$\frac{p(x_I^{{}^{\prime }}|x_{-I}^{{}^{\prime }})}{p(x_I|x_{-I})}=\frac{p(x^{{}^{\prime }})}{p(x)}$$

MH演算法中可以利用上述性質，簡化運算，提高效率。具體地，通過滿條件分佈概率的比$\frac{p(x_I^{{}^{\prime }}|x_{-I}^{{}^{\prime }})}{p(x_I|x_{-I})}$計算聯合概率的比$\frac{p(x^{{}^{\prime }})}{p(x)}$，而前者更容易計算。

5. 基本步驟

①任意選擇一個初始值$x_0$；
②對$i=1,2,\cdots ,n$迴圈執行：
  (a)設狀態$x_{i-1}=x$，按照建議分佈$q(x,x^{{}^{\prime }})$隨機抽取一個候選狀態$x^{{}^{\prime }}$；
  (b)計算接收概率：

$$\alpha (x,x^{{}^{\prime }})=min\left\{1,\frac{p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)}{p(x)q(x,x^{{}^{\prime }})}\right\}$$

  ©從區間$(0,1)$中按均勻分佈隨機抽取一個數$u$。若$u\le \alpha (x,x^{{}^{\prime }})$，則狀態$x_i=x^{{}^{\prime }}$，否則，狀態$x_i=x$；

③得到樣本集合 { x m + 1 , x m + 2 , ⋯   , x n } \left \{x_{m+1},x_{m+2},\cdots ,x_{n}\right \} {xm+1​,xm+2​,⋯,xn​}，計算函式樣本均值$f_{mn}$：

$$f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)$$

6. 單分量MH演算法

在MH演算法中，通常需要對多元變數分佈進行抽樣，有時對多元變數的抽樣是困難的。可以對多元變數的每一變數的條件分佈依次分別進行抽樣，從而實現對整個多元變數的一次抽樣，這就是單分量MH（single-component Metropolis-Hastings）演算法。

假設馬爾可夫鏈的狀態有$k$維隨機變數表示：

$$x=(x_1,x_2,\cdots ,x_k{)}^T$$

為了生成容量為$n$的樣本集合$\left\{x^{(1)},x^{(2)},\cdots ,x^{(n)}\right\}$，單分量MH演算法由下面的$k$步迭代實現MH演算法的一次迭代。

設在第$i-1$次迭代結束時分量$x_j$的取值為$x_j^{(i-1)}$，在第$i$次迭代的第$j$步，對分量$x_j$根據MH演算法更新，得到其新的取值$x_j^{(i)}$。首先，由建議分佈$q(x_j^{(i-1)},x_j|x_{-j}^{(i)})$抽樣產生分量$x_j$的候選值$x_j^{{}^{\prime }(i)}$，這裡$x_{-j}^{(i)}$表示在第$i$次迭代的第$j-1$步後的$x^i$除去$x_j^{i-1}$的所有值，即：

$$x_{-j}^{(i)}=(x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\cdots ,x_k^{(i-1)}{)}^T$$

其中分量$1,2,\cdots ,j-1$已經更新。然後，按照接受概率：

$$\alpha (x_j^{(i-1)},x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})=min\left\{1,\frac{p(x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})q(x_j^{{}^{\prime }(i)},x_j^{(i-1)}|x_{-j}^{(i)})}{p(x_j^{(i-1)}|x_{-j}^{(i)})q(x_j^{(i-1)},x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})}\right\}$$

抽樣決定是否接受候選值$x_j^{{}^{\prime }(i)}$。如果$x_j^{{}^{\prime }(i)}$被接受，則令$x_j^{(i)}=x_j^{{}^{\prime }(i)}$；否則$x_j^{(i)}=x_j^{(i-1)}$。其餘分量在第$j$步不改變。馬爾可夫鏈的轉移概率為：

$$p(x_j^{(i-1)},x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})=\alpha (x_j^{(i-1)},x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})q(x_j^{(i-1)},x_j^{{}^{\prime }(i)}|x_{-j}^{(i)})$$

三、吉布斯抽樣

吉布斯抽樣可以認為是MH演算法的特殊情況，但是更容易實現，因此被廣泛使用。

1. 基本原理

吉布斯抽樣（Gabbs sampling） 用於多元變數聯合分佈的抽樣和估計。其基本做法是，從聯合概率分佈定義滿條件概率分佈，依次對滿條件概率分佈進行抽樣，得到樣本的序列。可以證明這樣的抽樣過程是在一個馬爾可夫鏈上的隨機遊走，每一個樣本對應著馬爾可夫鏈的狀態，平穩分佈就是目標的聯合分佈。整體成為一個MCMC，燃燒期之後的樣本就是聯合分佈的隨機樣本。

假設多元變數的聯合概率分佈為$p(x)=p(x_1,x_2,\cdots ,x_k)$。吉布斯抽樣從一個初始樣本$x^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^T$出發，不斷進行迭代，每一次迭代得到聯合概率分佈的一個樣本$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T$。最終得到樣本序列$\left\{x^{(0)},x^{(1)},\cdots ,x^{(n)}\right\}$。

在每次迭代中，依次對$k$個隨機變數中的一個變數進行隨機抽樣。如果在第$i$次迭代中，對第$j$個變數進行隨機抽樣，那麼抽樣的分佈就是滿條件概率分佈$p(x_j,x_{-j}^{(i)})$。

設在第$i-1$步得到樣本$(x_1^{(i-1)},x_2^{(i-1)},\cdots ,x_k^{(i-1)}{)}^T$，在第$i$步，首先對第一個變數按照以下滿條件概率分佈隨機抽樣：

$$p(x_1|x_2^{(i-1)},\cdots ,x_k^{(i-1)})$$

得到$x_1^{(i)}$，之後依次對第$j$個變數按照以下滿條件概率分佈隨機抽樣：

$$p(x_j|x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\cdots ,x_k^{(i-1)}),j=2,\cdots ,k-1$$

得到$x_j^{(i)}$，最後對第$k$個變數按照以下滿條件概率分佈隨機抽樣：

$$p(x_k|x_1^{(i)},\cdots ,x_{k-1}^{(i)})$$

得到$x_k^{(i)}$，於是得到樣本$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T$。

2. 吉布斯抽樣與單分量MH演算法的關係

吉布斯抽樣是單分量MH演算法的特殊情況。定義建議分佈是當前變數$x_j$，$j=1,2,\cdots ,k$的滿條件概率分佈：

$$q(x,x^{{}^{\prime }})=p(x_j^{{}^{\prime }}|x_{-j})$$

這時，接受概率$\alpha =1$，

$$\begin{gathered} \alpha (x,x^{{}^{\prime }})=min\left\{1,\frac{p(x^{{}^{\prime }})q(x^{{}^{\prime }},x)}{p(x)q(x,x^{{}^{\prime }})}\right\} \\ =min\left\{1,\frac{p(x_{-j}^{{}^{\prime }})p(x_j^{{}^{\prime }}|x_{-j}^{{}^{\prime }})p(x_j|x_{-j}^{{}^{\prime }})}{p(x_{-j})p(x_j|x_{-j})p(x_j^{{}^{\prime }}|x_{-j})}\right\} \\ =1 \end{gathered}$$

這裡用到了$p(x_{-j})=p(x_{-j}^{{}^{\prime }})$和$p(\cdot |x_{-j})=p(\cdot |x_{-j}^{{}^{\prime }})$。

轉移核就是滿條件概率分佈：

$$p(x,x^{{}^{\prime }})=p(x_j^{{}^{\prime }}|x_{-j})$$

也就是說按照單分量的滿條件概率分佈$p(x_j^{{}^{\prime }}|x_{-j})$進行隨機抽樣，就能實現單分量MH演算法。吉布斯抽樣對每次抽樣的結果都接受，沒有拒絕，這一點和一般的MH演算法不同。

這裡，假設滿條件概率分佈$p(x_j^{{}^{\prime }}|x_{-j})$不為$0$，即馬爾可夫鏈是不可約的。

3. 基本步驟

①初始化。給出初始樣本$x^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^T$；
②對$i=1,2,\cdots ,n$迴圈執行：
  (1)由滿條件分佈$p(x_1|x_2^{(i-1)},\cdots ,x_k^{(i-1)})$抽取$x_1^{(i)}$；
   $⋮$
  (j)由滿條件分佈$p(x_j|x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\cdots ,x_k^{(i-1)})$抽取$x_j^{(i)}$；
   $⋮$
  (k)由滿條件分佈$p(x_k|x_1^{(i)},\cdots ,x_{k-1}^{(i)})$抽取$x_k^{(i)}$。
得到第$i$次迭代值$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T$。
③得到樣本集合$\left\{x^{(m+1)},x^{(m+2)},\cdots ,x^{(n)}\right\}$，計算函式樣本均值$f_{mn}$：

$$f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x^{(i)})$$

4. 對比單分量MH演算法

單分量MH演算法與吉布斯抽樣的不同之處在於，在前者演算法中，抽樣會在樣本點之間移動，但其間可能在某一些樣本點上停留（由於取樣被拒絕）；而在後者演算法中，抽樣點會在樣本點之間持續移動。

吉布斯抽樣適合於滿條件概率分佈容易抽樣的情況，而單分量MH演算法適合於滿條件概率分佈不容易抽樣的情況，這是使用容易抽樣的條件分佈作建議分佈。

參考資料

ref:李航《統計學習方法》