3.1 隨機模擬

隨機模擬(或者統計模擬)方法有一個很酷的別名是蒙特卡羅方法(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。

隨機模擬與計算機

現代的統計模擬方法最早由數學家烏拉姆提出，被Metropolis命名為蒙特卡羅方法，蒙特卡羅是著名的賭場，賭博總是和統計密切關聯的，所以這個命名風趣而貼切，很快被大家廣泛接受。被不過據說費米之前就已經在實驗中使用了，但是沒有發表。說起蒙特卡羅方法的源頭，可以追溯到18世紀，布豐當年用于計算[Math Processing Error]π的著名的投針實驗就是蒙特卡羅模擬實驗。統計采樣的方法其實數學家們很早就知道，但是在計算機出現以前，隨機數生成的成本很高，所以該方法也沒有實用價值。隨著計算機技術在二十世紀后半葉的迅猛發展，隨機模擬技術很快進入實用階段。對那些用確定算法不可行或不可能解決的問題，蒙特卡羅方法常常為人們帶來希望。

蒙特卡羅方法

統計模擬中有一個重要的問題就是給定一個概率分布[Math Processing Error]p(x)，我們如何在計算機中生成它的樣本。一般而言均勻分布?[Math Processing Error]Uniform(0,1)的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成[Math Processing Error][0,1]之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布?[Math Processing Error]Uniform(0,1)?的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。

生成一個概率分布的樣本

而我們常見的概率分布，無論是連續的還是離散的分布，都可以基于[Math Processing Error]Uniform(0,1)?的樣本生成。例如正態分布可以通過著名的 Box-Muller 變換得到

[Box-Muller 變換] ?如果隨機變量?[Math Processing Error]U1,U2?獨立且[Math Processing Error]U1,U2～Uniform[0,1]，
[Math Processing Error]Z0=?2ln?U1cos(2πU2)Z1=?2ln?U1sin(2πU2)
則?[Math Processing Error]Z0,Z1?獨立且服從標準正態分布。

其它幾個著名的連續分布，包括指數分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通過類似的數學變換得到；離散的分布通過均勻分布更加容易生成。更多的統計分布如何通過均勻分布的變換生成出來，大家可以參考統計計算的書，其中 Sheldon M. Ross 的《統計模擬》是寫得非常通俗易懂的一本。

不過我們并不是總是這么幸運的，當[Math Processing Error]p(x)的形式很復雜，或者[Math Processing Error]p(x)?是個高維的分布的時候，樣本的生成就可能很困難了。 譬如有如下的情況

• [Math Processing Error]p(x)=p~(x)∫p~(x)dx,而?[Math Processing Error]p~(x)?我們是可以計算的，但是底下的積分式無法顯式計算。
• [Math Processing Error]p(x,y)?是一個二維的分布函數，這個函數本身計算很困難，但是條件分布?[Math Processing Error]p(x|y),p(y|x)的計算相對簡單;如果?[Math Processing Error]p(x)?是高維的，這種情形就更加明顯。

此時就需要使用一些更加復雜的隨機模擬的方法來生成樣本。而本節中將要重點介紹的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一種，這兩個方法在現代貝葉斯分析中被廣泛使用。要了解這兩個算法，我們首先要對馬氏鏈的平穩分布的性質有基本的認識。

3.2 馬氏鏈及其平穩分布

馬氏鏈的數學定義很簡單
[Math Processing Error]P(Xt+1=x|Xt,Xt?1,?)=P(Xt+1=x|Xt)
也就是狀態轉移的概率只依賴于前一個狀態。

我們先來看馬氏鏈的一個具體的例子。社會學家經常把人按其經濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1,2,3 分別代表這三個階層。社會學家們發現決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是 0.65, 屬于中層收入的概率是 0.28, 屬于上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉移概率如下

使用矩陣的表示方式，轉移概率矩陣記為
[Math Processing Error]P=[0.650.280.070.150.670.180.120.360.52]

假設當前這一代人處在下層、中層、上層的人的比例是概率分布向量?[Math Processing Error]π0=[π0(1),π0(2),π0(3)]，那么他們的子女的分布比例將是?[Math Processing Error]π1=π0P, 他們的孫子代的分布比例將是?[Math Processing Error]π2=π1P=π0P2, ……, 第[Math Processing Error]n代子孫的收入分布比例將是?[Math Processing Error]πn=πn?1P=π0Pn。

假設初始概率分布為[Math Processing Error]π0=[0.21,0.68,0.11]，則我們可以計算前[Math Processing Error]n代人的分布狀況如下

我們發現從第7代人開始，這個分布就穩定不變了，這個是偶然的嗎？我們換一個初始概率分布[Math Processing Error]π0=[0.75,0.15,0.1]?試試看，繼續計算前[Math Processing Error]n代人的分布狀況如下

我們發現，到第9代人的時候, 分布又收斂了。最為奇特的是，兩次給定不同的初始概率分布，最終都收斂到概率分布?[Math Processing Error]π=[0.286,0.489,0.225]，也就是說收斂的行為和初始概率分布?[Math Processing Error]π0?無關。這說明這個收斂行為主要是由概率轉移矩陣[Math Processing Error]P決定的。我們計算一下?[Math Processing Error]Pn

[Math Processing Error]P20=P21=?=P100=?=[0.2860.4890.2250.2860.4890.2250.2860.4890.225]

我們發現，當?[Math Processing Error]n?足夠大的時候，這個[Math Processing Error]Pn矩陣的每一行都是穩定地收斂到[Math Processing Error]π=[0.286,0.489,0.225]?這個概率分布。自然的，這個收斂現象并非是我們這個馬氏鏈獨有的，而是絕大多數馬氏鏈的共同行為，關于馬氏鏈的收斂我們有如下漂亮的定理：

馬氏鏈定理：?如果一個非周期馬氏鏈具有轉移概率矩陣[Math Processing Error]P,且它的任何兩個狀態是連通的，那么?[Math Processing Error]limn→∞Pijn?存在且與[Math Processing Error]i無關，記?[Math Processing Error]limn→∞Pijn=π(j), 我們有

[Math Processing Error]limn→∞Pn=[π(1)π(2)?π(j)?π(1)π(2)?π(j)??????π(1)π(2)?π(j)??????]

?[Math Processing Error]π(j)=∑i=0∞π(i)Pij

?[Math Processing Error]π?是方程?[Math Processing Error]πP=π?的唯一非負解

其中,
[Math Processing Error]π=[π(1),π(2),?,π(j),?],∑i=0∞πi=1
[Math Processing Error]π稱為馬氏鏈的平穩分布。

這個馬氏鏈的收斂定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以這個定理作為理論基礎的。?定理的證明相對復雜，一般的隨機過程課本中也不給證明，所以我們就不用糾結它的證明了，直接用這個定理的結論就好了。我們對這個定理的內容做一些解釋說明：

該定理中馬氏鏈的狀態不要求有限，可以是有無窮多個的；

定理中的“非周期“這個概念我們不打算解釋了，因為我們遇到的絕大多數馬氏鏈都是非周期的；

兩個狀態[Math Processing Error]i,j是連通并非指[Math Processing Error]i?可以直接一步轉移到[Math Processing Error]j([Math Processing Error]Pij>0),而是指?[Math Processing Error]i?可以通過有限的[Math Processing Error]n步轉移到達[Math Processing Error]j([Math Processing Error]Pijn>0)。馬氏鏈的任何兩個狀態是連通的含義是指存在一個[Math Processing Error]n, 使得矩陣[Math Processing Error]Pn?中的任何一個元素的數值都大于零。

我們用?[Math Processing Error]Xi?表示在馬氏鏈上跳轉第[Math Processing Error]i步后所處的狀態，如果?[Math Processing Error]limn→∞Pijn=π(j)?存在，很容易證明以上定理的第二個結論。由于
[Math Processing Error]P(Xn+1=j)=∑i=0∞P(Xn=i)P(Xn+1=j|Xn=i)=∑i=0∞P(Xn=i)Pij
上式兩邊取極限就得到?[Math Processing Error]π(j)=∑i=0∞π(i)Pij

從初始概率分布?[Math Processing Error]π0?出發，我們在馬氏鏈上做狀態轉移，記[Math Processing Error]Xi的概率分布為[Math Processing Error]πi, 則有
[Math Processing Error]X0～π0(x)Xi～πi(x),πi(x)=πi?1(x)P=π0(x)Pn
由馬氏鏈收斂的定理, 概率分布[Math Processing Error]πi(x)將收斂到平穩分布?[Math Processing Error]π(x)。假設到第[Math Processing Error]n步的時候馬氏鏈收斂，則有
[Math Processing Error]X0～π0(x)X1～π1(x)?Xn～πn(x)=π(x)Xn+1～π(x)Xn+2～π(x)?
所以?[Math Processing Error]Xn,Xn+1,Xn+2,?～π(x)?都是同分布的隨機變量，當然他們并不獨立。如果我們從一個具體的初始狀態?[Math Processing Error]x0?開始,沿著馬氏鏈按照概率轉移矩陣做跳轉，那么我們得到一個轉移序列?[Math Processing Error]x0,x1,x2,?xn,xn+1?,?由于馬氏鏈的收斂行為，[Math Processing Error]xn,xn+1,??都將是平穩分布?[Math Processing Error]π(x)?的樣本。

3.3 Markov Chain Monte Carlo

對于給定的概率分布[Math Processing Error]p(x),我們希望能有便捷的方式生成它對應的樣本。由于馬氏鏈能收斂到平穩分布， 于是一個很的漂亮想法是：如果我們能構造一個轉移矩陣為[Math Processing Error]P的馬氏鏈，使得該馬氏鏈的平穩分布恰好是[Math Processing Error]p(x), 那么我們從任何一個初始狀態[Math Processing Error]x0出發沿著馬氏鏈轉移, 得到一個轉移序列?[Math Processing Error]x0,x1,x2,?xn,xn+1?,， 如果馬氏鏈在第[Math Processing Error]n步已經收斂了，于是我們就得到了?[Math Processing Error]π(x)?的樣本[Math Processing Error]xn,xn+1?。

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

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

定理：[細致平穩條件]?如果非周期馬氏鏈的轉移矩陣[Math Processing Error]P和分布[Math Processing Error]π(x)?滿足
[Math Processing Error](1)π(i)Pij=π(j)Pjifor alli,j
則?[Math Processing Error]π(x)?是馬氏鏈的平穩分布，上式被稱為細致平穩條件(detailed balance condition)。

其實這個定理是顯而易見的，因為細致平穩條件的物理含義就是對于任何兩個狀態[Math Processing Error]i,j, 從?[Math Processing Error]i?轉移出去到[Math Processing Error]j?而丟失的概率質量，恰好會被從?[Math Processing Error]j?轉移回[Math Processing Error]i?的概率質量補充回來，所以狀態[Math Processing Error]i上的概率質量[Math Processing Error]π(i)是穩定的，從而[Math Processing Error]π(x)是馬氏鏈的平穩分布。數學上的證明也很簡單，由細致平穩條件可得
[Math Processing Error]∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=π(j)∑i=1∞Pji=π(j)?πP=π
由于[Math Processing Error]π?是方程?[Math Processing Error]πP=π的解，所以[Math Processing Error]π是平穩分布。

假設我們已經有一個轉移矩陣為[Math Processing Error]Q馬氏鏈([Math Processing Error]q(i,j)表示從狀態?[Math Processing Error]i轉移到狀態[Math Processing Error]j的概率，也可以寫為?[Math Processing Error]q(j|i)或者[Math Processing Error]q(i→j)), 顯然，通常情況下?[Math Processing Error]p(i)q(i,j)≠p(j)q(j,i)?也就是細致平穩條件不成立，所以?[Math Processing Error]p(x)?不太可能是這個馬氏鏈的平穩分布。我們可否對馬氏鏈做一個改造，使得細致平穩條件成立呢？譬如，我們引入一個?[Math Processing Error]α(i,j), 我們希望
[Math Processing Error](2)p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)(?)
取什么樣的?[Math Processing Error]α(i,j)?以上等式能成立呢？最簡單的，按照對稱性，我們可以取
[Math Processing Error]α(i,j)=p(j)q(j,i)，α(j,i)=p(i)q(i,j)
于是(*)式就成立了。所以有
[Math Processing Error](3)p(i)q(i,j)α(i,j)?Q′(i,j)=p(j)q(j,i)α(j,i)?Q′(j,i)(??)
于是我們把原來具有轉移矩陣[Math Processing Error]Q的一個很普通的馬氏鏈，改造為了具有轉移矩陣[Math Processing Error]Q′的馬氏鏈，而?[Math Processing Error]Q′恰好滿足細致平穩條件，由此馬氏鏈[Math Processing Error]Q′的平穩分布就是[Math Processing Error]p(x)！

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

馬氏鏈轉移和接受概率

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

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

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

假設?[Math Processing Error]α(i,j)=0.1,α(j,i)=0.2, 此時滿足細致平穩條件，于是
[Math Processing Error]p(i)q(i,j)×0.1=p(j)q(j,i)×0.2
上式兩邊擴大5倍，我們改寫為
[Math Processing Error]p(i)q(i,j)×0.5=p(j)q(j,i)×1
看，我們提高了接受率，而細致平穩條件并沒有打破！這啟發我們可以把細致平穩條件(**) 式中的[Math Processing Error]α(i,j),α(j,i)?同比例放大，使得兩數中最大的一個放大到1，這樣我們就提高了采樣中的跳轉接受率。所以我們可以取
[Math Processing Error]α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}
于是，經過對上述MCMC 采樣算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。

對于分布?[Math Processing Error]p(x),我們構造轉移矩陣?[Math Processing Error]Q′?使其滿足細致平穩條件
[Math Processing Error]p(x)Q′(x→y)=p(y)Q′(y→x)
此處?[Math Processing Error]x?并不要求是一維的，對于高維空間的?[Math Processing Error]p(x)，如果滿足細致平穩條件
[Math Processing Error]p(x)Q′(x→y)=p(y)Q′(y→x)
那么以上的 Metropolis-Hastings 算法一樣有效。

3.2 Gibbs Sampling

對于高維的情形，由于接受率?[Math Processing Error]α的存在(通常?[Math Processing Error]α<1), 以上 Metropolis-Hastings 算法的效率不夠高。能否找到一個轉移矩陣Q使得接受率?[Math Processing Error]α=1?呢？我們先看看二維的情形，假設有一個概率分布[Math Processing Error]p(x,y), 考察[Math Processing Error]x坐標相同的兩個點[Math Processing Error]A(x1,y1),B(x1,y2)，我們發現
[Math Processing Error]p(x1,y1)p(y2|x1)=p(x1)p(y1|x1)p(y2|x1)p(x1,y2)p(y1|x1)=p(x1)p(y2|x1)p(y1|x1)
所以得到
[Math Processing Error](4)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)(???)
即
[Math Processing Error]p(A)p(y2|x1)=p(B)p(y1|x1)
基于以上等式，我們發現，在?[Math Processing Error]x=x1?這條平行于?[Math Processing Error]y軸的直線上，如果使用條件分布?[Math Processing Error]p(y|x1)做為任何兩個點之間的轉移概率，那么任何兩個點之間的轉移滿足細致平穩條件。同樣的，如果我們在[Math Processing Error]y=y1?這條直線上任意取兩個點?[Math Processing Error]A(x1,y1),C(x2,y1),也有如下等式
[Math Processing Error]p(A)p(x2|y1)=p(C)p(x1|y1).

平面上馬氏鏈轉移矩陣的構造

于是我們可以如下構造平面上任意兩點之間的轉移概率矩陣Q
[Math Processing Error]Q(A→B)=p(yB|x1)如果xA=xB=x1Q(A→C)=p(xC|y1)如果yA=yC=y1Q(A→D)=0其它

有了如上的轉移矩陣 Q, 我們很容易驗證對平面上任意兩點?[Math Processing Error]X,Y, 滿足細致平穩條件
[Math Processing Error]p(X)Q(X→Y)=p(Y)Q(Y→X)
于是這個二維空間上的馬氏鏈將收斂到平穩分布?[Math Processing Error]p(x,y)。而這個算法就稱為 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 這兩兄弟于1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個算法在現代貝葉斯分析中占據重要位置。

Gibbs Sampling 算法中的馬氏鏈轉移

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

以上的過程我們很容易推廣到高維的情形，對于(***) 式，如果[Math Processing Error]x1變為多維情形[Math Processing Error]x1，可以看出推導過程不變，所以細致平穩條件同樣是成立的
[Math Processing Error](5)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)
此時轉移矩陣 Q 由條件分布?[Math Processing Error]p(y|x1)?定義。上式只是說明了一根坐標軸的情形，和二維情形類似，很容易驗證對所有坐標軸都有類似的結論。所以[Math Processing Error]n維空間中對于概率分布?[Math Processing Error]p(x1,x2,?,xn)?可以如下定義轉移矩陣

如果當前狀態為[Math Processing Error](x1,x2,?,xn)，馬氏鏈轉移的過程中，只能沿著坐標軸做轉移。沿著?[Math Processing Error]xi?這根坐標軸做轉移的時候，轉移概率由條件概率?[Math Processing Error]p(xi|x1,?,xi?1,xi+1,?,xn)?定義；

?其它無法沿著單根坐標軸進行的跳轉，轉移概率都設置為 0。

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

以上算法收斂后，得到的就是概率分布[Math Processing Error]p(x1,x2,?,xn)的樣本，當然這些樣本并不獨立，但是我們此處要求的是采樣得到的樣本符合給定的概率分布，并不要求獨立。同樣的，在以上算法中，坐標軸輪換采樣不是必須的，可以在坐標軸輪換中引入隨機性，這時候轉移矩陣?[Math Processing Error]Q?中任何兩個點的轉移概率中就會包含坐標軸選擇的概率，而在通常的 Gibbs Sampling 算法中，坐標軸輪換是一個確定性的過程，也就是在給定時刻[Math Processing Error]t，在一根固定的坐標軸上轉移的概率是1。

轉載地址：https://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

這篇鏈接http://www.myexception.cn/cloud/2100158.html也值得學習。

總結

以上是生活随笔為你收集整理的(转载)机器学习知识点(十三)吉布斯采样法(Gibbs Sampling)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。