摘要：算法中的馬氏鏈轉(zhuǎn)移以上采樣過程中，如圖所示，馬氏鏈的轉(zhuǎn)移只是輪換的沿著坐標(biāo)軸軸和軸做轉(zhuǎn)移，于是得到樣本馬氏鏈?zhǔn)諗亢螅罱K得到的樣本就是的樣本，而收斂之前的階段稱為。

作者：chen_h
微信號 & QQ：862251340
微信公眾號：coderpai
簡書地址：https://www.jianshu.com/p/278...

本文是整理網(wǎng)上的幾篇博客和論文所得出來的，所有的原文連接都在文末。

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

Metropolis-Hasting算法和Gibbs Sampling算法是馬爾科夫鏈蒙特卡洛（Markov Chain Mento Carlo，MCMC）方法。

MCMC方法是用蒙特卡洛方法去體現(xiàn)馬爾科夫鏈的方法。在講MCMC之前，必須要先講一下馬爾科夫鏈，馬爾鏈的數(shù)學(xué)定義為：

也就是說當(dāng)前狀態(tài)只和前一個狀態(tài)有關(guān)，而與其他狀態(tài)無關(guān)，Markov Chain 體現(xiàn)的是狀態(tài)空間的轉(zhuǎn)換關(guān)系，下一個狀態(tài)只和當(dāng)前狀態(tài)有關(guān)。比如下圖就是一個馬爾科夫鏈的示意圖：

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

在上圖中，我們假設(shè)當(dāng)前狀態(tài)為 u = [0.5, 0.2, 0.3]，那么下一個矩陣的狀態(tài)就是 u*p = [0.18, 0.64, 0.18]，再下一個矩陣的狀態(tài)就是 u*p*p = [ 0.108, 0.316, 0.576]，依照這個轉(zhuǎn)換矩陣一直轉(zhuǎn)換下去，最后的系統(tǒng)就趨近于一個穩(wěn)定狀態(tài) [0.22, 0.41, 0.37]。而事實證明無論你從哪個點出發(fā)，經(jīng)過很長的 Markov Chain 之后都會穩(wěn)定到這一點。

>>> 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

再舉一個例子，社會學(xué)家經(jīng)常把人按其經(jīng)濟狀況分為3類：下層（lower-class）、中層（middle-class）、上層（upper-class）。我們用1,2,3分別代表這三個階層。社會學(xué)家們發(fā)現(xiàn)決定一個人的收入階層的最重要的因素就是其父母的結(jié)束階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是 0.65，屬于中層收入的概率是 0.28，屬于上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉(zhuǎn)移概率如下：

同理，圖中的轉(zhuǎn)換關(guān)系可以用一個概率轉(zhuǎn)換矩陣 p 表示，如下：

假設(shè)當(dāng)前這一代人處在下層、中層、上層的概率分布向量是：

那么他們的子女的分布比例將是：

他們的孫子代的分布比例將是：

以此類推，第 n 代子孫的收入分布比例將是：

我們舉個具體的例子，假設(shè)初始概率分布為：

則我們可以計算前 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]

我們發(fā)現(xiàn)從第7代人開始，這個分布就穩(wěn)定不變了，事實上，在這個問題中，從任意初始概率分布開始都會收斂到這個穩(wěn)定的結(jié)果。也就是說，收斂的行為和初始概率分布 π0 無關(guān)。這說明這個收斂行為主要是由概率轉(zhuǎn)移矩陣 P 決定的。我們可以計算一下 P^n ：

我們發(fā)現(xiàn)，當(dāng) n 足夠大的時候，這個 P^n 矩陣的每一行都是穩(wěn)定地收斂到 π=[0.286,0.489,0.225] 這個概率分布。自然的，這個收斂現(xiàn)象并非是我們這個馬氏鏈都有的，而是絕大多數(shù)馬氏鏈的共同行為，關(guān)于馬氏鏈的收斂我們有如下漂亮的定理：

這個馬氏鏈的收斂定理非常重要，所有的MCMC（Markov Chain Monte Carlo）方法都是以這個定理作為理論基礎(chǔ)的。

對于給定的概率分布 p(x) ，我們希望能與便捷的方式生成它對應(yīng)的樣本。由于馬氏鏈能收斂到平穩(wěn)分布，于是一個很漂亮的想法就是：如果我們能構(gòu)造一個轉(zhuǎn)移矩陣為 P 的馬氏鏈，使得該馬氏鏈的平穩(wěn)分布恰好是 p(x)，那么我們從任何一個初始狀態(tài) x(0) 出發(fā)沿著馬氏鏈轉(zhuǎn)移，得到一個轉(zhuǎn)移序列 x(0)，x(1)，x(2)，...，x(n)，x(n+1)，...，如果馬氏鏈在第 n 步已經(jīng)收斂了，于是我們就得到了 π(x) 的樣本 x(n)，x(n+1)，...。

這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統(tǒng)的平穩(wěn)性質(zhì)， Metropolis 考慮了物理學(xué)中常見的波爾茲曼分布的采樣問題，首次提出了基于馬氏鏈的蒙特卡羅方法，即Metropolis算法，并在最早的計算機上編程實現(xiàn)。Metropolis 算法是首個普適的采樣方法，并啟發(fā)了一系列 MCMC方法，所以人們把它視為隨機模擬技術(shù)騰飛的起點。 Metropolis的這篇論文被收錄在《統(tǒng)計學(xué)中的重大突破》中， Metropolis算法也被遴選為二十世紀(jì)的十個最重要的算法之一。

我們接下來介紹的MCMC 算法是 Metropolis 算法的一個改進變種，即常用的 Metropolis-Hastings 算法。由上一節(jié)的例子和定理我們看到了，馬氏鏈的收斂性質(zhì)主要由轉(zhuǎn)移矩陣 P 決定, 所以基于馬氏鏈做采樣的關(guān)鍵問題是如何構(gòu)造轉(zhuǎn)移矩陣 P，使得平穩(wěn)分布恰好是我們要的分布p(x)。如何能做到這一點呢？我們主要使用如下的定理。

其實這個定理是顯而易見的，因為細(xì)致平穩(wěn)條件的物理含義就是對于任何兩個狀態(tài) i，j， 從 i 轉(zhuǎn)移出去到 j 而丟失的概率質(zhì)量，恰好會被從 j 轉(zhuǎn)移回ii 的概率質(zhì)量補充回來，所以狀態(tài)ii上的概率質(zhì)量 π(i) 是穩(wěn)定的，從而 π(x) 是馬氏鏈的平穩(wěn)分布。數(shù)學(xué)上的證明也很簡單，由細(xì)致平穩(wěn)條件可得：

由于 π 是方程 πP=π 的解，所以 π 是平穩(wěn)分布。

假設(shè)我們已經(jīng)有一個轉(zhuǎn)移矩陣為 Q 的馬氏鏈（q(i,j)表示從狀態(tài) i 轉(zhuǎn)移到狀態(tài) j 的概率，也可以寫為 q(j|i) 或者q(i→j)），顯然，通常情況下：

也就是細(xì)致平穩(wěn)條件不成立，所以 p(x) 不太可能是這個馬氏鏈的平穩(wěn)分布。我們可否對馬氏鏈做一個改造，使得細(xì)致平穩(wěn)條件成立呢？譬如，我們引入一個 α(i,j), 我們希望：

取什么樣的 α(i,j) 以上等式能成立呢？最簡單的，按照對稱性，我們可以取：

于是(*)式就成立了。所以有：

于是我們把原來具有轉(zhuǎn)移矩陣 Q 的一個很普通的馬氏鏈，改造為了具有轉(zhuǎn)移矩陣 Q′ 的馬氏鏈，而 Q′ 恰好滿足細(xì)致平穩(wěn)條件，由此馬氏鏈 Q′ 的平穩(wěn)分布就是 p(x)！

在改造 Q 的過程中引入的 α(i,j) 稱為接受率，物理意義可以理解為在原來的馬氏鏈上，從狀態(tài) i 以 q(i,j) 的概率轉(zhuǎn)跳轉(zhuǎn)到狀態(tài) j 的時候，我們以 α(i,j) 的概率接受這個轉(zhuǎn)移，于是得到新的馬氏鏈 Q′ 的轉(zhuǎn)移概率為 q(i,j)α(i,j) 。

假設(shè)我們已經(jīng)有一個轉(zhuǎn)移矩陣 Q (對應(yīng)元素為 q(i,j) ), 把以上的過程整理一下，我們就得到了如下的用于采樣概率分布 p(x) 的算法。

上述過程中 p(x),q(x|y) 說的都是離散的情形，事實上即便這兩個分布是連續(xù)的，以上算法仍然是有效，于是就得到更一般的連續(xù)概率分布 p(x)的采樣算法，而 q(x|y) 就是任意一個連續(xù)二元概率分布對應(yīng)的條件分布。

以上的 MCMC 采樣算法已經(jīng)能很漂亮的工作了，不過它有一個小的問題：馬氏鏈 Q 在轉(zhuǎn)移的過程中的接受率 α(i,j) 可能偏小，這樣采樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉(zhuǎn)，這使得馬氏鏈遍歷所有的狀態(tài)空間要花費太長的時間，收斂到平穩(wěn)分布 p(x) 的速度太慢。有沒有辦法提升一些接受率呢?

假設(shè) α(i,j)=0.1,α(j,i)=0.2, 此時滿足細(xì)致平穩(wěn)條件，于是

上式兩邊擴大5倍，我們改寫為

看，我們提高了接受率，而細(xì)致平穩(wěn)條件并沒有打破！這啟發(fā)我們可以把細(xì)致平穩(wěn)條件(**) 式中的 α(i,j),α(j,i) 同比例放大，使得兩數(shù)中最大的一個放大到1，這樣我們就提高了采樣中的跳轉(zhuǎn)接受率。所以我們可以取：

這個公式可以進一步簡化為下面的公式：

于是，經(jīng)過對上述MCMC 采樣算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。

對于分布 p(x)，我們構(gòu)造轉(zhuǎn)移矩陣 Q′ 使其滿足細(xì)致平穩(wěn)條件

此處 x 并不要求是一維的，對于高維空間的 p(x)，如果滿足細(xì)致平穩(wěn)條件

那么以上的 Metropolis-Hastings 算法一樣有效。

對于高維的情形，由于接受率 α 的存在(通常 α<1), 以上 Metropolis-Hastings 算法的效率不夠高。能否找到一個轉(zhuǎn)移矩陣 Q 使得接受率 α=1 呢？我們先看看二維的情形，假設(shè)有一個概率分布 p(x,y), 考察 x坐標(biāo)相同的兩個點 A(x1,y1),B(x1,y2)，我們發(fā)現(xiàn)

所以得到

即

基于以上等式，我們發(fā)現(xiàn)，在 x=x1 這條平行于 y 軸的直線上，如果使用條件分布 p(y|x1) 做為任何兩個點之間的轉(zhuǎn)移概率，那么任何兩個點之間的轉(zhuǎn)移滿足細(xì)致平穩(wěn)條件。同樣的，如果我們在 y=y1y=y1 這條直線上任意取兩個點 A(x1,y1),C(x2,y1)，也有如下等式

于是我們可以如下構(gòu)造平面上任意兩點之間的轉(zhuǎn)移概率矩陣 Q

有了如上的轉(zhuǎn)移矩陣 Q, 我們很容易驗證對平面上任意兩點 X,Y, 滿足細(xì)致平穩(wěn)條件

于是這個二維空間上的馬氏鏈將收斂到平穩(wěn)分布 p(x,y)。而這個算法就稱為 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 這兩兄弟于1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個算法在現(xiàn)代貝葉斯分析中占據(jù)重要位置。

Gibbs Sampling 算法中的馬氏鏈轉(zhuǎn)移

以上采樣過程中，如圖所示，馬氏鏈的轉(zhuǎn)移只是輪換的沿著坐標(biāo)軸 x 軸和 y 軸做轉(zhuǎn)移，于是得到樣本 (x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),? 馬氏鏈?zhǔn)諗亢螅罱K得到的樣本就是 p(x,y) 的樣本，而收斂之前的階段稱為 burn-in period。額外說明一下，我們看到教科書上的 Gibbs Sampling 算法大都是坐標(biāo)軸輪換采樣的，但是這其實是不強制要求的。最一般的情形可以是，在 t 時刻，可以在 x 軸和 y 軸之間隨機的選一個坐標(biāo)軸，然后按條件概率做轉(zhuǎn)移，馬氏鏈也是一樣收斂的。輪換兩個坐標(biāo)軸只是一種方便的形式。

以上的過程我們很容易推廣到高維的情形，對于(*) 式，如果 x1 變?yōu)槎嗑S情形 x1，可以看出推導(dǎo)過程不變，所以細(xì)致平穩(wěn)條件同樣是成立的

此時轉(zhuǎn)移矩陣 Q 由條件分布 p(y|x1) 定義。上式只是說明了一根坐標(biāo)軸的情形和二維情形類似，很容易驗證對所有坐標(biāo)軸都有類似的結(jié)論。所以 n 維空間中對于概率分布 p(x1,x2,?,xn) 可以如下定義轉(zhuǎn)移矩陣

如果當(dāng)前狀態(tài)為 (x1,x2,?,xn)，馬氏鏈轉(zhuǎn)移的過程中，只能沿著坐標(biāo)軸做轉(zhuǎn)移。沿著 xi 這根坐標(biāo)軸做轉(zhuǎn)移的時候，轉(zhuǎn)移概率由條件概率 p(xi|x1,?,xi?1,xi+1,?,xn) 定義；

其它無法沿著單根坐標(biāo)軸進行的跳轉(zhuǎn)，轉(zhuǎn)移概率都設(shè)置為 0。

于是我們可以把Gibbs Smapling 算法從采樣二維的 p(x,y) 推廣到采樣 n 維的 p(x1,x2,?,xn)

以上算法收斂后，得到的就是概率分布 p(x1,x2,?,xn) 的樣本，當(dāng)然這些樣本并不獨立，但是我們此處要求的是采樣得到的樣本符合給定的概率分布，并不要求獨立。同樣的，在以上算法中，坐標(biāo)軸輪換采樣不是必須的，可以在坐標(biāo)軸輪換中引入隨機性，這時候轉(zhuǎn)移矩陣 Q 中任何兩個點的轉(zhuǎn)移概率中就會包含坐標(biāo)軸選擇的概率，而在通常的 Gibbs Sampling 算法中，坐標(biāo)軸輪換是一個確定性的過程，也就是在給定時刻tt，在一根固定的坐標(biāo)軸上轉(zhuǎn)移的概率是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

作者：chen_h
微信號 & QQ：862251340
簡書地址：https://www.jianshu.com/p/278...

CoderPai 是一個專注于算法實戰(zhàn)的平臺，從基礎(chǔ)的算法到人工智能算法都有設(shè)計。如果你對算法實戰(zhàn)感興趣，請快快關(guān)注我們吧。加入AI實戰(zhàn)微信群，AI實戰(zhàn)QQ群，ACM算法微信群，ACM算法QQ群。長按或者掃描如下二維碼，關(guān)注 “CoderPai” 微信號（coderpai）