在科學(xué)研究中，如何生成服從某個(gè)概率分布的樣本是一個(gè)重要的問題。 如果樣本維度很低，只有一兩維，我們可以用反切法、拒絕采樣和重要性采樣等方法。 但是對(duì)于高維樣本，這些方法就不適用了。這時(shí)我們就要使用一些 “高檔” 的算法，比如下面要介紹的 Metropolis-Hasting 算法和 Gibbs sampling 算法。

      Metropolis-Hasting 算法和 Gibbs sampling 算法是馬爾科夫鏈蒙特卡洛（Markov Chain Mento Carlo, MCMC）方法。我們先介紹 MCMC 方法。

1. 馬爾科夫蒙特卡洛方法

      MCMC 方法是用蒙特卡洛方法去體現(xiàn)馬爾科夫鏈的方法。馬爾科夫鏈?zhǔn)菭顟B(tài)空間的轉(zhuǎn)換關(guān)系，下一個(gè)狀態(tài)只和當(dāng)前的狀態(tài)有關(guān)。比如下圖就是一個(gè)馬爾科夫鏈的示意圖。

圖中轉(zhuǎn)移關(guān)系可以用一個(gè)概率轉(zhuǎn)換矩陣 p 表示，

如果當(dāng)前狀態(tài)分布為 (u(x)), 那么下一個(gè)矩陣的狀態(tài)就是 ( u(x)p ), 再下一個(gè)就是(u(x)p^2),... 最后會(huì)收斂到一個(gè)平穩(wěn)分布 (\pi)。這個(gè)平穩(wěn)分布 (\pi) 只和概率轉(zhuǎn)移矩陣 p 有關(guān)，而和初始狀態(tài)分布 u 是什么沒有關(guān)系。

      如何判斷一個(gè)馬爾科夫鏈?zhǔn)欠衲苁諗康狡椒€(wěn)分布，以及如何判斷一個(gè)狀態(tài)分布是不是一個(gè)馬爾科夫鏈的平穩(wěn)分布呢？我們有下面定理。

細(xì)致平衡條件: 已知各態(tài)歷經(jīng)的的馬爾科夫鏈有概率轉(zhuǎn)移矩陣  ，以及已知狀態(tài)分布  。如果對(duì)于任意兩個(gè)狀態(tài) i 和 j，下面公式成立，則馬爾科夫鏈能夠收斂到  。

這里的各態(tài)歷經(jīng)是指任意兩個(gè)狀態(tài)之間可以通過有限步到達(dá)。

怎么證明細(xì)致平衡條件呢？我也不知道啊。

      MCMC 方法的基本原理是利用細(xì)致平衡條件構(gòu)建一個(gè)概率轉(zhuǎn)移矩陣，使得目標(biāo)概率就是概率轉(zhuǎn)移矩陣的平穩(wěn)分布。 Metropolis-Hasting 和 Gibbs sampling 算法本質(zhì)上是構(gòu)建概率轉(zhuǎn)移矩陣的不同方法。

2. Metropolis-Hastings 算法

      Metropolis-Hastings 算法先提出一個(gè)可能不符合條件的概率轉(zhuǎn)移矩陣 q, 然后再進(jìn)行調(diào)整。比如我們提出的 q 是均勻概率，即從任意狀態(tài)到任意狀態(tài)的概率是相等的。顯然在絕大部分情況下，q 的穩(wěn)定概率不是目標(biāo)概率 (\pi)，即不滿足細(xì)致平衡條件。

如何讓這個(gè)不等式轉(zhuǎn)變成等式呢？根據(jù)對(duì)稱性，我們?nèi)菀椎玫较旅娴牡仁健?br>

這時(shí)整個(gè)概率轉(zhuǎn)移矩陣滿足細(xì)致平衡條件。從 i 狀態(tài)轉(zhuǎn)到 j 狀態(tài)的概率是 (q(j|i) \pi(j)q(i|j))，實(shí)現(xiàn)這個(gè)轉(zhuǎn)移概率的方式是 i 狀態(tài)以 q(j|i) 概率跳轉(zhuǎn)到 j 狀態(tài)，然后以 (\pi(j)q(i|j)) 接受跳轉(zhuǎn) (拒絕跳轉(zhuǎn)就退回 i 狀態(tài))。這樣整個(gè) Metropolis-Hasting 算法的框架就建立起來了。

      這個(gè)原始的 Metropoli-Hasting 算法的有一個(gè)小問題。 跳轉(zhuǎn)接受概率 (\pi(j)q(i|j)) 和 (\pi(i)q(j|i)) 的值很小，算法進(jìn)行過程充斥著跳轉(zhuǎn)拒絕。為了改進(jìn)這點(diǎn)，Metropoli-Hasting 算法的方法是公式兩邊同時(shí)乘以一個(gè)系數(shù)，使得 (\pi(j)q(i|j)) 和 (\pi(i)q(j|i)) 中大的一項(xiàng) scale 到 1，得到下面的公式。

這個(gè)公式可以進(jìn)一步簡(jiǎn)化為下面的公式

      根據(jù)上面的推導(dǎo)，我們?nèi)菀椎玫?Metropolis-Hasting 算法的流程。

3. Gibbs sampling 算法

      Gibbs sampling 算法是 Metropolis-Hasting 算法的一個(gè)特例。很雞賊的一個(gè)特例。m 維的一個(gè)樣本跳轉(zhuǎn)到另一個(gè)樣本的過程，可以拆解為 m 個(gè)子過程，每一個(gè)子過程對(duì)應(yīng)一個(gè)維度。這時(shí)概率轉(zhuǎn)移矩陣是 m 個(gè)子概率轉(zhuǎn)移矩陣之積，即 (p = \prod_{i=k}^{m} p_k )

其中 (p_k) 表示第 k 維的變化概率。在 (p_k) 中，兩個(gè)狀態(tài)之間只有 k 維不同，其跳轉(zhuǎn)概率如下所示；不然為 0。

\begin{eqnarray}
p_k(\pmb{x}{\dashv k, k=v_2}|\pmb{x}{\dashv k, k=v_1}) = \frac{\pi(\pmb{x}{\dashv k, k=v_2})}{\sum{v}\pi(\pmb{x}_{\dashv k, k=v})} \nonumber
\end{eqnarray}

其中 (\pmb{x}_{\dashv k, k=v_2}) 表示樣本第 k 維數(shù)據(jù)為 (v_2)，其它維度固定。這時(shí)候我們發(fā)現(xiàn)如下公式

\begin{eqnarray}
&& \pi(\pmb{x}{\dashv k, k=v_1}) p(\pmb{x}{\dashv k, k=v_2}|\pmb{x}{\dashv k, k=v_1}) \nonumber \
&=& \pi(\pmb{x}{\dashv k, k=v_1}) \frac{\pi(\pmb{x}{\dashv k, k=v_2})}{\sum{v}\pi(\pmb{x}{\dashv k, k=v})} \nonumber \
&=& \pi(\pmb{x}{\dashv k, k=v_2}) \frac{\pi(\pmb{x}{\dashv k, k=v_1})}{\sum{v}\pi(\pmb{x}{\dashv k, k=v})} \nonumber \
&=& \pi(\pmb{x}{\dashv k, k=v_2}) p(\pmb{x}{\dashv k, k=v_1}|\pmb{x}{\dashv k, k=v_2}) \nonumber \
\end{eqnarray}
即 (p_k) 和 (\pi) 滿足細(xì)致平衡條件的等式。那么 (p_k) 就是我們要構(gòu)建的概率轉(zhuǎn)移矩陣嘛？答案是否定的。因?yàn)橥暾募?xì)致平衡條件需要各態(tài)歷經(jīng)。在概率轉(zhuǎn)移矩陣 (p_k) 下, 只有 k 維數(shù)據(jù)子啊變化，因此一個(gè)狀態(tài)永遠(yuǎn)不能到達(dá)和它第 k-1 維數(shù)據(jù)不同的狀態(tài)。

      最終我們構(gòu)建的概率轉(zhuǎn)移矩陣是 m 個(gè)子概率轉(zhuǎn)移矩陣之積

我們很容易證明 (p) 依然滿足細(xì)致平衡條件中的等式，同時(shí)還滿足各態(tài)歷經(jīng)。根據(jù)上述推導(dǎo)，我們得到 Gibbs sampling 的算法過程。

4. 總結(jié)

      Metropolist-Hasting 和 Gibbs sampling 是有效的 MCMC 算法，能夠解決高維空間的采樣問題。