原文傳送門：也談MCMC方法與Gibbs抽樣?
MCMC，即傳說中的Markov Chain Mento Carlo方法。其主要用于統(tǒng)計(jì)推理中進(jìn)行模擬抽樣，尤其在貝葉斯推理中有著非常廣泛的應(yīng)用。如算法模型的后驗(yàn)參數(shù)估計(jì)問題，很多情況下其后驗(yàn)概率分布沒有確定性的解析解，或者解析解計(jì)算起來非常復(fù)雜，便可以通過MCMC模擬抽樣，根據(jù)大數(shù)定律，參數(shù)的期望便可以通過對抽樣樣本的求均值來評估。
山人第一次見到MCMC兄還是在研究僧階段，那時(shí)候以Latent Direichlet Allocation(LDA)為代表的Blei先生的一系列主題模型算法還很火，甚至你還能看見Andrew Ng的身影。于是導(dǎo)師欣然的把其另一篇層次主題模型的論文，Hierarchical LDA(hLDA)甩給我們，拍著我們的肩膀，語重心長的說，好好干，會(huì)很有前景的。于是我的MCMC初體驗(yàn)是這樣的：
What the hell? 于是直到現(xiàn)在還對MCMC念念不忘。好吧，是耿耿于懷。最近又看見Quora上有人討論MCMC和Gibbs抽樣，再看時(shí)，發(fā)現(xiàn)雖然有一兩年未看，腦部神經(jīng)元還是不停的工作，現(xiàn)在理解起來竟然清晰許多。 MCMC是Markov Chain和Mento Carlo兩個(gè)概念的組合，我們不妨分而治之，先看看各自的含義。
I-Markov Chain
即馬爾科夫鏈，這哥么大家肯定不會(huì)陌生，還記得Hidden Markov Model么（Baum-Welch算法會(huì)推導(dǎo)了么:( )馬爾科夫鏈的一個(gè)重要屬性就是無記憶性。其表示的隨機(jī)過程，在一個(gè)狀態(tài)空間里游走且未來的狀態(tài)只與當(dāng)前的狀態(tài)有關(guān)，而與之前的狀態(tài)均無關(guān)。這種無記憶性便稱之為馬爾科夫性。
p(xt+1|xt,xt?1...x1)=p(xt+1|xt)(1)
馬爾科夫鏈?zhǔn)且环N隨機(jī)過程，其定義有主要有兩點(diǎn)，即狀態(tài)空間和轉(zhuǎn)移概率矩陣。如下圖所示，一個(gè)簡單的馬爾科夫鏈隨機(jī)過程，包含三個(gè)狀態(tài)：
其狀態(tài)之間的轉(zhuǎn)移概率矩陣如下：
假設(shè)在狀態(tài)Πi時(shí)，你在Bull Market 狀態(tài)，且當(dāng)前概率分布為[0,1,0]。在下一個(gè)Πi+1狀態(tài)時(shí)的概率分布為
Πi+1=Πi.P(2)
則結(jié)果為Πi+1=[.15.8.05]。如此類推，下一個(gè)狀態(tài)分布則為：
Πi+1=Πi.P2(3)
如此下去，最終發(fā)現(xiàn)我們會(huì)得到一個(gè)穩(wěn)定的狀態(tài)，此時(shí)
Π=Π.P(4)
即狀態(tài)分布變得穩(wěn)定(Stationary)，不會(huì)再隨著狀態(tài)轉(zhuǎn)移概率的變化而變化。且我們發(fā)現(xiàn)，即使我們的初始狀態(tài)分布矩陣不是[0,1,0]而是另外一個(gè)值，如[0.4,0.3,0.3]時(shí)，最終經(jīng)過多次轉(zhuǎn)移，也會(huì)達(dá)到最終的穩(wěn)定(Stationary)狀態(tài)，且穩(wěn)定狀態(tài)的分布是一致的，即最終的Stationary狀態(tài)與初始分布矩陣沒有關(guān)系，只與狀態(tài)轉(zhuǎn)移矩陣有關(guān)。那末是不是所有的狀態(tài)轉(zhuǎn)移矩陣都能最終達(dá)到穩(wěn)定狀態(tài)呢？答案自然不是，還是需要馬氏鏈定理的保證，簡單說就是
如果一個(gè)非周期馬氏鏈具有概率轉(zhuǎn)移矩陣P，且它的任何兩個(gè)狀態(tài)都是聯(lián)通的，那么如果limn?>∞Pnij=π(j)存在且僅與j有關(guān)，那么這樣的一個(gè)穩(wěn)定分布就是存在的。
這里還有一點(diǎn)山人剛開始時(shí)也是非常模糊。就是很多算法中提到，當(dāng)經(jīng)過了burn-in階段，狀態(tài)分布穩(wěn)定以后開始取樣計(jì)算概率分布，當(dāng)時(shí)就想，既然都穩(wěn)定了，π都保持不變了，取的樣本不都一樣么？其實(shí)這里所說的狀態(tài)穩(wěn)定是指滿足了某一個(gè)概率分布，即穩(wěn)定后抽樣出的樣本都是同分布的。而在穩(wěn)定之前則可能不同的樣本是產(chǎn)生自不同的概率分布。
II-Monte Carlo
說完了馬爾科夫再來說說蒙特卡洛方法吧，其名子來源于摩納哥的蒙特卡洛賭場，是一種通過模擬抽樣求積分的方法。一個(gè)經(jīng)典的應(yīng)用便是計(jì)算圓周率。這個(gè)名叫“hit and miss"的實(shí)驗(yàn)過程為：假設(shè)有一個(gè)單位長度為1的正方形區(qū)域，再以正方形的中心為圓心，單位長度為半徑畫一個(gè)正方形的內(nèi)切圓。有一個(gè)隨機(jī)數(shù)發(fā)射器隨機(jī)的往正方形區(qū)域里發(fā)射。當(dāng)經(jīng)過N多次以后，圓周率可以估算為(hawaii.edu)：π=4NhitNshot
大學(xué)微積分中我們學(xué)過常見函數(shù)求積分的方法，如I=∫∞θg(θ)p(θ)dθ，p是θ的概率密度函數(shù)，求其在g上的積分。但在實(shí)際應(yīng)用中，函數(shù)g往往是不可積的，且θ可能是高緯向量，使得我們很難求得其解析解。在大數(shù)定律和中心極限定理的保證下，蒙特卡洛方法則通過模擬抽樣的方法為求其近似解提供了一條途徑。我們可以通過從概率密度函數(shù)p中抽樣出θ，最終MC近似的解為：I′=∑Mi=1g(θi)。
應(yīng)用到貝葉斯推理中，如果我們能夠通過抽樣的方式從參數(shù)變量的聯(lián)合分布中抽取到足夠多的樣本數(shù)據(jù)，我們便可以通過貝葉斯參數(shù)估計(jì)等方法求得其近似值。但往往參數(shù)的聯(lián)合分布各個(gè)變量并非獨(dú)立，且很復(fù)雜。尤其如LDA等主題生成模型里，要對聯(lián)合分布抽樣幾乎是不可能的。有么有可能通過某種控制變量法，對條件概率進(jìn)行抽樣，借用馬爾科夫鏈中條件概率轉(zhuǎn)移矩陣達(dá)到穩(wěn)定狀態(tài)后的概率分布就是其變量的聯(lián)合分布下的樣本點(diǎn)呢？
III-MCMC類方法
于是，為了避免構(gòu)造一個(gè)復(fù)雜繁瑣的聯(lián)合分布函數(shù)來進(jìn)行蒙特卡洛抽樣，MCMC類方法神兵天降。通過構(gòu)造一個(gè)狀態(tài)轉(zhuǎn)移概率矩陣，那末當(dāng)其到達(dá)穩(wěn)定狀態(tài)時(shí)，分布便是所求的聯(lián)合概率分布。而聯(lián)合分布函數(shù)的樣本點(diǎn)則是每一次狀態(tài)轉(zhuǎn)移時(shí)自然產(chǎn)生的。這么牛掰的想法當(dāng)然不是山人想到的，一個(gè)叫著Metropolis的哥么在1953年研究粒子系統(tǒng)的平穩(wěn)性質(zhì)便提出來了。而目前我們常用的一個(gè)叫著Metropolis-Hastings算法便是在其基礎(chǔ)上的一個(gè)改進(jìn)。
1 細(xì)致平穩(wěn)條件
我們在前面提到了，我們可以通過構(gòu)造一個(gè)狀態(tài)轉(zhuǎn)移概率矩陣，使得其平穩(wěn)狀態(tài)下的概率分布就是我們想要的分布。但不是隨意構(gòu)造一個(gè)狀態(tài)轉(zhuǎn)移概率矩陣就能滿足的。那需要什么樣的條件呢？細(xì)致平穩(wěn)條件就是這樣一個(gè)充分條件。如果非周期馬氏鏈的轉(zhuǎn)移概率矩陣P和分布π(x)滿足：
π(i)Pij=π(j)Pjiforalli,j(5)
則π(x)就是該馬氏鏈的平穩(wěn)分布。那自然不是所有的概率矩陣和分布都滿足等式（5）中的條件，我們可以通過對馬氏鏈進(jìn)行一個(gè)小小的改造：
π(i)Pijα(i,j)=π(j)Pjiα(j,i)(6)
于是新得到的馬氏鏈為P′(j,i)：
π(i)p(i,j)α(i,j)P′(i,j)=π(j)p(j,i)α(j,i)P′(j,i)(??)(7)
而只要通過對稱性，取α(i,j)為π(j)p(j,i)，取α(j,i)為π(i)p(i,j)即可。此處的α(.)稱之為接受率。其可以理解為，在原來的馬氏鏈上，從狀態(tài)i以p(i,j)的概率跳轉(zhuǎn)到狀態(tài)j時(shí)，我們以α(i,j)的概率接受這個(gè)跳轉(zhuǎn)。
一般的MCMC采樣算法的接受率通過和一個(gè)Uniform[0,1]分布采樣的值u作比較，如果接受率大于這個(gè)值，則接受這次轉(zhuǎn)移，從i轉(zhuǎn)移到j(luò)狀態(tài)，反之則保持原i狀態(tài)。但是我們在實(shí)際應(yīng)用中使用這個(gè)方法時(shí)發(fā)現(xiàn)，很多情況下接受率普遍很低，導(dǎo)致馬氏鏈狀態(tài)轉(zhuǎn)移緩慢，最終收斂的速度非常慢。為了解決這個(gè)問題，我們還是采用類似等式（6）的方法，分子分母的接受率同步增大。
α(x,y)α(y,x)=π(y)p(y,x)π(x)q(x,y)
我們可以把跳轉(zhuǎn)之后的狀態(tài)α(y,x)接受率為1，則我們可以得到下面的接受率公式(注意接受率取值范圍只能是[0,1])：
α(i,j)=min{π(j)p(j,i)π(i)p(i,j),1}(8)
按照式（8）的接受率，便是我們的Metropolis-Hastings算法。
2 Gibbs抽樣
當(dāng)變量狀態(tài)多，且維度比較高時(shí)，MH算法的接受率仍然差強(qiáng)人意。要是每次都接受該多好啊。那什么樣的情況下，我從i到j(luò)時(shí)，每次都能接受呢？（即接受率為1）。最終發(fā)現(xiàn)，我們每次可以沿著垂直于某個(gè)變量維度的軸走。即通過迭代的方法，每一次只對一個(gè)變量進(jìn)行采樣。舉一個(gè)二維空間的例子，假設(shè)一個(gè)概率分布p(x,y)，來看x坐標(biāo)相同的兩個(gè)點(diǎn)A(x1,y1)和B(x1,y2)，通過簡單的聯(lián)合概率和條件概率的關(guān)系我們可以得到：
p(x1,y1)p(y2|x1)=p(x1)p(y1|x1)p(y2|x1)(9)
p(x1,y2)p(y1|x1)=p(x1)p(y2|x1)p(y1|x1)(10)
很明顯，等式(9),(10)右邊是相等的，如(11)所示:
p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)(11)
下圖給出了一個(gè)更直觀的表示：
即，從A到B和從B到A的轉(zhuǎn)移是直接滿足細(xì)致平穩(wěn)條件的。因此我們不需要等式(6)中的接受率來幫忙，即接受率為1.圖中假設(shè)初始狀態(tài)為A，則從A到下一個(gè)概率轉(zhuǎn)移矩陣分別為：
Q(A?>B)=p(yB|x1)
Q(A?>C)=p(xC|y1)
Q(A?>D)=0(12)
因此類似于曼哈頓距離的方法，狀態(tài)轉(zhuǎn)移總是沿著橫平豎直的街區(qū)進(jìn)行。這邊是Gibbs抽樣算法的核心思想。下圖給出了一個(gè)Gibbs抽樣的直觀圖。
3 收斂條件的判斷
我們都知道當(dāng)概率狀態(tài)轉(zhuǎn)移穩(wěn)定時(shí)，其分布便是所要求的聯(lián)合概率分布。但我們不可能通過如等式(2),(3)的方法來每轉(zhuǎn)換一步就求其概率分布，比較是否改變。主要原因有二，其一是不可把所有變量間的轉(zhuǎn)移概率都找到，其二矩陣計(jì)算耗時(shí)耗力。常見的方法便是通過burn-in的方法，多跑幾次。也有通過計(jì)算當(dāng)前狀態(tài)下的聯(lián)合分布可能性函數(shù)，然后根據(jù)Autocorrelation Function(ACF)的變化速率來判斷迭代是否收斂。
So long, and thanks for all the fish.
參考
[1]PRML讀書會(huì)第十一章 Sampling Methods
[2]LDA-math-MCMC 和 Gibbs Sampling
[5]What-are-Markov-Chain-Monte-Carlo-methods-in-laymans-terms
