作者：chen_h
微訊號 & QQ：862251340
微信公眾號：coderpai
我的部落格：請點選這裡
本文是整理網上的幾篇部落格和論文所得出來的，所有的原文連線都在文末。
在科學研究中，如何生成服從某個概率分佈的樣本是一個重要的問題。如果樣本維度很低，只有一兩維，我們可以用反切法，拒絕取樣和重要性取樣等方法。但是對於高位樣本，這些方法就不適用了。這時我們就可以使用一些“高檔”的演算法，比如Metropolis-Hasting演算法和Gibbs Sampling演算法。
Metropolis-Hasting演算法和Gibbs Sampling演算法是馬爾科夫鏈蒙特卡洛（Markov Chain Mento Carlo，MCMC）方法。
1. 馬爾科夫鏈蒙特卡洛（MCMC）方法
MCMC方法是用蒙特卡洛方法去體現馬爾科夫鏈的方法。在講MCMC之前，必須要先講一下馬爾科夫鏈，馬爾鏈的數學定義為：
也就是說當前狀態只和前一個狀態有關，而與其他狀態無關，Markov Chain 體現的是狀態空間的轉換關係，下一個狀態只和當前狀態有關。比如下圖就是一個馬爾科夫鏈的示意圖：
圖中轉換關係可以用一個概率轉換矩陣 p 表示，如下：
在上圖中，我們假設當前狀態為 u = [0.5, 0.2, 0.3]，那麼下一個矩陣的狀態就是 u*p = [0.18, 0.64, 0.18]，再下一個矩陣的狀態就是 u*p*p = [ 0.108, 0.316, 0.576]，依照這個轉換矩陣一直轉換下去，最後的系統就趨近於一個穩定狀態 [0.22, 0.41, 0.37]。而事實證明無論你從哪個點出發，經過很長的 Markov Chain 之後都會穩定到這一點。
>>> u
array([ 0.5, 0.2, 0.3])
>>> p
array([[ 0. , 1. , 0. ],
[ 0. , 0.1, 0.9],
[ 0.6, 0.4, 0. ]])
>>> for _ in xrange(100):
... u = np.dot(u,p)
... print u
複製程式碼再舉一個例子，社會學家經常把人按其經濟狀況分為3類：下層（lower-class）、中層（middle-class）、上層（upper-class）。我們用1,2,3分別代表這三個階層。社會學家們發現決定一個人的收入階層的最重要的因素就是其父母的結束階層。如果一個人的收入屬於下層類別，那麼他的孩子屬於下層收入的概率是 0.65，屬於中層收入的概率是 0.28，屬於上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉移概率如下：
同理，圖中的轉換關係可以用一個概率轉換矩陣 p 表示，如下：
假設當前這一代人處在下層、中層、上層的概率分佈向量是：
那麼他們的子女的分佈比例將是：
他們的孫子代的分佈比例將是：
以此類推，第 n 代子孫的收入分佈比例將是：
我們舉個具體的例子，假設初始概率分佈為：
則我們可以計算前 n 代人的分佈狀況如下：
>>> x
array([ 0.21, 0.68, 0.11])
>>> p
array([[ 0.65, 0.28, 0.07],
[ 0.15, 0.67, 0.18],
[ 0.12, 0.36, 0.52]])
>>> for _ in xrange(10):
... x = np.dot(x,p)
... print x
...
[ 0.2517 0.554 0.1943]
[ 0.270021 0.511604 0.218375]
[ 0.27845925 0.49699556 0.22454519]
[ 0.28249327 0.49179188 0.22571485]
[ 0.28447519 0.48985602 0.22566879]
[ 0.28546753 0.48909735 0.22543512]
[ 0.28597071 0.48878278 0.22524651]
[ 0.28622796 0.488645 0.22512704]
[ 0.28636017 0.48858171 0.22505812]
[ 0.28642834 0.48855152 0.22502014]複製程式碼我們發現從第7代人開始，這個分佈就穩定不變了，事實上，在這個問題中，從任意初始概率分佈開始都會收斂到這個穩定的結果。也就是說，收斂的行為和初始概率分佈 π0 無關。這說明這個收斂行為主要是由概率轉移矩陣 P 決定的。我們可以計算一下 P^n ：
我們發現，當 n 足夠大的時候，這個 P^n 矩陣的每一行都是穩定地收斂到 π=[0.286,0.489,0.225] 這個概率分佈。自然的，這個收斂現象並非是我們這個馬氏鏈都有的，而是絕大多數馬氏鏈的共同行為，關於馬氏鏈的收斂我們有如下漂亮的定理：
這個馬氏鏈的收斂定理非常重要，所有的MCMC（Markov Chain Monte Carlo）方法都是以這個定理作為理論基礎的。
對於給定的概率分佈 p(x) ，我們希望能與便捷的方式生成它對應的樣本。由於馬氏鏈能收斂到平穩分佈，於是一個很漂亮的想法就是：如果我們能構造一個轉移矩陣為 P 的馬氏鏈，使得該馬氏鏈的平穩分佈恰好是 p(x)，那麼我們從任何一個初始狀態 x(0) 出發沿著馬氏鏈轉移，得到一個轉移序列 x(0)，x(1)，x(2)，...，x(n)，x(n+1)，...，如果馬氏鏈在第 n 步已經收斂了，於是我們就得到了 π(x) 的樣本 x(n)，x(n+1)，...。
這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統的平穩性質， Metropolis 考慮了物理學中常見的波爾茲曼分佈的取樣問題，首次提出了基於馬氏鏈的蒙特卡羅方法，即Metropolis演算法，並在最早的計算機上程式設計實現。Metropolis 演算法是首個普適的取樣方法，並啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis的這篇論文被收錄在《統計學中的重大突破》中， Metropolis演算法也被遴選為二十世紀的十個最重要的演算法之一。
我們接下來介紹的MCMC 演算法是 Metropolis 演算法的一個改進變種，即常用的 Metropolis-Hastings 演算法。由上一節的例子和定理我們看到了，馬氏鏈的收斂性質主要由轉移矩陣 P 決定, 所以基於馬氏鏈做取樣的關鍵問題是如何構造轉移矩陣 P，使得平穩分佈恰好是我們要的分佈p(x)。如何能做到這一點呢？我們主要使用如下的定理。
其實這個定理是顯而易見的，因為細緻平穩條件的物理含義就是對於任何兩個狀態 i，j， 從 i 轉移出去到 j 而丟失的概率質量，恰好會被從 j 轉移回ii 的概率質量補充回來，所以狀態ii上的概率質量 π(i) 是穩定的，從而 π(x) 是馬氏鏈的平穩分佈。數學上的證明也很簡單，由細緻平穩條件可得：
由於 π 是方程 πP=π 的解，所以 π 是平穩分佈。
假設我們已經有一個轉移矩陣為 Q 的馬氏鏈（q(i,j)表示從狀態 i 轉移到狀態 j 的概率，也可以寫為 q(j|i) 或者q(i→j)），顯然，通常情況下：
也就是細緻平穩條件不成立，所以 p(x) 不太可能是這個馬氏鏈的平穩分佈。我們可否對馬氏鏈做一個改造，使得細緻平穩條件成立呢？譬如，我們引入一個 α(i,j), 我們希望：
取什麼樣的 α(i,j) 以上等式能成立呢？最簡單的，按照對稱性，我們可以取：
於是(*)式就成立了。所以有：
於是我們把原來具有轉移矩陣 Q 的一個很普通的馬氏鏈，改造為了具有轉移矩陣 Q′ 的馬氏鏈，而 Q′ 恰好滿足細緻平穩條件，由此馬氏鏈 Q′ 的平穩分佈就是 p(x)！
在改造 Q 的過程中引入的 α(i,j) 稱為接受率，物理意義可以理解為在原來的馬氏鏈上，從狀態 i 以 q(i,j) 的概率轉跳轉到狀態 j 的時候，我們以 α(i,j) 的概率接受這個轉移，於是得到新的馬氏鏈 Q′ 的轉移概率為 q(i,j)α(i,j) 。
假設我們已經有一個轉移矩陣 Q (對應元素為 q(i,j) ), 把以上的過程整理一下，我們就得到了如下的用於取樣概率分佈 p(x) 的演算法。
上述過程中 p(x),q(x|y) 說的都是離散的情形，事實上即便這兩個分佈是連續的，以上演算法仍然是有效，於是就得到更一般的連續概率分佈 p(x)的取樣演算法，而 q(x|y) 就是任意一個連續二元概率分佈對應的條件分佈。
以上的 MCMC 取樣演算法已經能很漂亮的工作了，不過它有一個小的問題：馬氏鏈 Q 在轉移的過程中的接受率 α(i,j) 可能偏小，這樣取樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉，這使得馬氏鏈遍歷所有的狀態空間要花費太長的時間，收斂到平穩分佈 p(x) 的速度太慢。有沒有辦法提升一些接受率呢?
假設 α(i,j)=0.1,α(j,i)=0.2, 此時滿足細緻平穩條件，於是
上式兩邊擴大5倍，我們改寫為
看，我們提高了接受率，而細緻平穩條件並沒有打破！這啟發我們可以把細緻平穩條件(**) 式中的 α(i,j),α(j,i) 同比例放大，使得兩數中最大的一個放大到1，這樣我們就提高了取樣中的跳轉接受率。所以我們可以取：
這個公式可以進一步簡化為下面的公式：
於是，經過對上述MCMC 取樣演算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 演算法。
對於分佈 p(x)，我們構造轉移矩陣 Q′ 使其滿足細緻平穩條件
此處 x 並不要求是一維的，對於高維空間的 p(x)，如果滿足細緻平穩條件
那麼以上的 Metropolis-Hastings 演算法一樣有效。
2. Gibbs Sampling演算法
對於高維的情形，由於接受率 α 的存在(通常 α<1), 以上 Metropolis-Hastings 演算法的效率不夠高。能否找到一個轉移矩陣 Q 使得接受率 α=1 呢？我們先看看二維的情形，假設有一個概率分佈 p(x,y), 考察 x座標相同的兩個點 A(x1,y1),B(x1,y2)，我們發現
所以得到
即
基於以上等式，我們發現，在 x=x1 這條平行於 y 軸的直線上，如果使用條件分佈 p(y|x1) 做為任何兩個點之間的轉移概率，那麼任何兩個點之間的轉移滿足細緻平穩條件。同樣的，如果我們在 y=y1y=y1 這條直線上任意取兩個點 A(x1,y1),C(x2,y1)，也有如下等式
於是我們可以如下構造平面上任意兩點之間的轉移概率矩陣 Q
有了如上的轉移矩陣 Q, 我們很容易驗證對平面上任意兩點 X,Y, 滿足細緻平穩條件
於是這個二維空間上的馬氏鏈將收斂到平穩分佈 p(x,y)。而這個演算法就稱為 Gibbs Sampling 演算法,是 Stuart Geman 和Donald Geman 這兩兄弟於1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個演算法在現代貝葉斯分析中佔據重要位置。
Gibbs Sampling 演算法中的馬氏鏈轉移
以上取樣過程中，如圖所示，馬氏鏈的轉移只是輪換的沿著座標軸 x 軸和 y 軸做轉移，於是得到樣本 (x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),⋯ 馬氏鏈收斂後，最終得到的樣本就是 p(x,y) 的樣本，而收斂之前的階段稱為 burn-in period。額外說明一下，我們看到教科書上的 Gibbs Sampling 演算法大都是座標軸輪換取樣的，但是這其實是不強制要求的。最一般的情形可以是，在 t 時刻，可以在 x 軸和 y 軸之間隨機的選一個座標軸，然後按條件概率做轉移，馬氏鏈也是一樣收斂的。輪換兩個座標軸只是一種方便的形式。
以上的過程我們很容易推廣到高維的情形，對於(
此時轉移矩陣 Q 由條件分佈 p(y|x1) 定義。上式只是說明了一根座標軸的情形和二維情形類似，很容易驗證對所有座標軸都有類似的結論。所以 n 維空間中對於概率分佈 p(x1,x2,⋯,xn) 可以如下定義轉移矩陣
如果當前狀態為 (x1,x2,⋯,xn)，馬氏鏈轉移的過程中，只能沿著座標軸做轉移。沿著 xi 這根座標軸做轉移的時候，轉移概率由條件概率 p(xi|x1,⋯,xi−1,xi+1,⋯,xn) 定義；
其它無法沿著單根座標軸進行的跳轉，轉移概率都設定為 0。
於是我們可以把Gibbs Smapling 演算法從取樣二維的 p(x,y) 推廣到取樣 n 維的 p(x1,x2,⋯,xn)
以上演算法收斂後，得到的就是概率分佈 p(x1,x2,⋯,xn) 的樣本，當然這些樣本並不獨立，但是我們此處要求的是取樣得到的樣本符合給定的概率分佈，並不要求獨立。同樣的，在以上演算法中，座標軸輪換取樣不是必須的，可以在座標軸輪換中引入隨機性，這時候轉移矩陣 Q 中任何兩個點的轉移概率中就會包含座標軸選擇的概率，而在通常的 Gibbs Sampling 演算法中，座標軸輪換是一個確定性的過程，也就是在給定時刻tt，在一根固定的座標軸上轉移的概率是1。
Reference：
Metropolis-Hastings 和 Gibbs sampling
隨機取樣方法整理與講解（MCMC、Gibbs Sampling等）
LDA-math-MCMC 和 Gibbs Sampling
Introduction to Monte Carlo Methods
An Introduction to MCMC for Machine Learning
CoderPai 是一個專注於演算法實戰的平臺，從基礎的演算法到人工智慧演算法都有設計。如果你對演算法實戰感興趣，請快快關注我們吧。加入AI實戰微信群，AI實戰QQ群，ACM演算法微信群，ACM演算法QQ群。詳情請關注 “CoderPai” 微訊號（coderpai）。