原文傳送門：也談MCMC方法與Gibbs抽樣?
MCMC，即傳說中的Markov Chain Mento Carlo方法姑子。其主要用于統(tǒng)計推理中進行模擬抽樣，尤其在貝葉斯推理中有著非常廣泛的應(yīng)用。如算法模型的后驗參數(shù)估計問題怪与，很多情況下其后驗概率分布沒有確定性的解析解，或者解析解計算起來非常復(fù)雜缅疟，便可以通過MCMC模擬抽樣分别，根據(jù)大數(shù)定律，參數(shù)的期望便可以通過對抽樣樣本的求均值來評估存淫。
山人第一次見到MCMC兄還是在研究僧階段耘斩，那時候以Latent Direichlet Allocation(LDA)為代表的Blei先生的一系列主題模型算法還很火，甚至你還能看見Andrew Ng的身影桅咆。于是導(dǎo)師欣然的把其另一篇層次主題模型的論文括授，Hierarchical LDA(hLDA)甩給我們，拍著我們的肩膀岩饼，語重心長的說荚虚，好好干，會很有前景的籍茧。于是我的MCMC初體驗是這樣的：
What the hell? 于是直到現(xiàn)在還對MCMC念念不忘版述。好吧，是耿耿于懷寞冯。最近又看見Quora上有人討論MCMC和Gibbs抽樣院水，再看時腊徙，發(fā)現(xiàn)雖然有一兩年未看，腦部神經(jīng)元還是不停的工作檬某，現(xiàn)在理解起來竟然清晰許多撬腾。 MCMC是Markov Chain和Mento Carlo兩個概念的組合，我們不妨分而治之恢恼，先看看各自的含義民傻。
I-Markov Chain
即馬爾科夫鏈，這哥么大家肯定不會陌生场斑，還記得Hidden Markov Model么（Baum-Welch算法會推導(dǎo)了么:( )馬爾科夫鏈的一個重要屬性就是無記憶性漓踢。其表示的隨機過程，在一個狀態(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隨機過程青责，其定義有主要有兩點挺据，即狀態(tài)空間和轉(zhuǎn)移概率矩陣。如下圖所示脖隶，一個簡單的馬爾科夫鏈隨機過程扁耐，包含三個狀態(tài)：
其狀態(tài)之間的轉(zhuǎn)移概率矩陣如下：
假設(shè)在狀態(tài)Πi時，你在Bull Market 狀態(tài)产阱，且當(dāng)前概率分布為[0,1,0]婉称。在下一個Πi+1狀態(tài)時的概率分布為
Πi+1=Πi.P(2)
則結(jié)果為Πi+1=[.15.8.05]。如此類推构蹬，下一個狀態(tài)分布則為：
Πi+1=Πi.P2(3)
如此下去王暗，最終發(fā)現(xiàn)我們會得到一個穩(wěn)定的狀態(tài)，此時
Π=Π.P(4)
即狀態(tài)分布變得穩(wěn)定(Stationary)庄敛，不會再隨著狀態(tài)轉(zhuǎn)移概率的變化而變化俗壹。且我們發(fā)現(xiàn)，即使我們的初始狀態(tài)分布矩陣不是[0,1,0]而是另外一個值铐姚，如[0.4,0.3,0.3]時策肝，最終經(jīng)過多次轉(zhuǎn)移，也會達到最終的穩(wěn)定(Stationary)狀態(tài)隐绵，且穩(wěn)定狀態(tài)的分布是一致的之众，即最終的Stationary狀態(tài)與初始分布矩陣沒有關(guān)系，只與狀態(tài)轉(zhuǎn)移矩陣有關(guān)依许。那末是不是所有的狀態(tài)轉(zhuǎn)移矩陣都能最終達到穩(wěn)定狀態(tài)呢棺禾？答案自然不是，還是需要馬氏鏈定理的保證峭跳，簡單說就是
如果一個非周期馬氏鏈具有概率轉(zhuǎn)移矩陣P膘婶，且它的任何兩個狀態(tài)都是聯(lián)通的缺前，那么如果limn?>∞Pnij=π(j)存在且僅與j有關(guān)，那么這樣的一個穩(wěn)定分布就是存在的悬襟。
這里還有一點山人剛開始時也是非常模糊衅码。就是很多算法中提到，當(dāng)經(jīng)過了burn-in階段脊岳，狀態(tài)分布穩(wěn)定以后開始取樣計算概率分布逝段，當(dāng)時就想，既然都穩(wěn)定了割捅，π都保持不變了奶躯，取的樣本不都一樣么？其實這里所說的狀態(tài)穩(wěn)定是指滿足了某一個概率分布亿驾，即穩(wěn)定后抽樣出的樣本都是同分布的嘹黔。而在穩(wěn)定之前則可能不同的樣本是產(chǎn)生自不同的概率分布。
II-Monte Carlo
說完了馬爾科夫再來說說蒙特卡洛方法吧莫瞬，其名子來源于摩納哥的蒙特卡洛賭場儡蔓，是一種通過模擬抽樣求積分的方法。一個經(jīng)典的應(yīng)用便是計算圓周率乏悄。這個名叫“hit and miss"的實驗過程為：假設(shè)有一個單位長度為1的正方形區(qū)域浙值，再以正方形的中心為圓心恳不，單位長度為半徑畫一個正方形的內(nèi)切圓檩小。有一個隨機數(shù)發(fā)射器隨機的往正方形區(qū)域里發(fā)射。當(dāng)經(jīng)過N多次以后烟勋，圓周率可以估算為(hawaii.edu)：π=4NhitNshot
大學(xué)微積分中我們學(xué)過常見函數(shù)求積分的方法规求，如I=∫∞θg(θ)p(θ)dθ，p是θ的概率密度函數(shù)卵惦，求其在g上的積分阻肿。但在實際應(yīng)用中，函數(shù)g往往是不可積的沮尿，且θ可能是高緯向量丛塌，使得我們很難求得其解析解。在大數(shù)定律和中心極限定理的保證下畜疾，蒙特卡洛方法則通過模擬抽樣的方法為求其近似解提供了一條途徑赴邻。我們可以通過從概率密度函數(shù)p中抽樣出θ，最終MC近似的解為：I′=∑Mi=1g(θi)啡捶。
應(yīng)用到貝葉斯推理中姥敛，如果我們能夠通過抽樣的方式從參數(shù)變量的聯(lián)合分布中抽取到足夠多的樣本數(shù)據(jù)，我們便可以通過貝葉斯參數(shù)估計等方法求得其近似值瞎暑。但往往參數(shù)的聯(lián)合分布各個變量并非獨立彤敛，且很復(fù)雜与帆。尤其如LDA等主題生成模型里，要對聯(lián)合分布抽樣幾乎是不可能的墨榄。有么有可能通過某種控制變量法玄糟，對條件概率進行抽樣，借用馬爾科夫鏈中條件概率轉(zhuǎn)移矩陣達到穩(wěn)定狀態(tài)后的概率分布就是其變量的聯(lián)合分布下的樣本點呢袄秩？
III-MCMC類方法
于是茶凳，為了避免構(gòu)造一個復(fù)雜繁瑣的聯(lián)合分布函數(shù)來進行蒙特卡洛抽樣播揪，MCMC類方法神兵天降。通過構(gòu)造一個狀態(tài)轉(zhuǎn)移概率矩陣箱沦，那末當(dāng)其到達穩(wěn)定狀態(tài)時，分布便是所求的聯(lián)合概率分布雇庙。而聯(lián)合分布函數(shù)的樣本點則是每一次狀態(tài)轉(zhuǎn)移時自然產(chǎn)生的谓形。這么牛掰的想法當(dāng)然不是山人想到的疆前，一個叫著Metropolis的哥么在1953年研究粒子系統(tǒng)的平穩(wěn)性質(zhì)便提出來了。而目前我們常用的一個叫著Metropolis-Hastings算法便是在其基礎(chǔ)上的一個改進竹椒。
1 細致平穩(wěn)條件
我們在前面提到了童太，我們可以通過構(gòu)造一個狀態(tài)轉(zhuǎn)移概率矩陣，使得其平穩(wěn)狀態(tài)下的概率分布就是我們想要的分布胸完。但不是隨意構(gòu)造一個狀態(tài)轉(zhuǎn)移概率矩陣就能滿足的书释。那需要什么樣的條件呢赊窥？細致平穩(wěn)條件就是這樣一個充分條件。如果非周期馬氏鏈的轉(zhuǎn)移概率矩陣P和分布π(x)滿足：
π(i)Pij=π(j)Pjiforalli,j(5)
則π(x)就是該馬氏鏈的平穩(wěn)分布扯再。那自然不是所有的概率矩陣和分布都滿足等式（5）中的條件址遇，我們可以通過對馬氏鏈進行一個小小的改造：
π(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時，我們以α(i,j)的概率接受這個跳轉(zhuǎn)十籍。
一般的MCMC采樣算法的接受率通過和一個Uniform[0,1]分布采樣的值u作比較唇礁，如果接受率大于這個值，則接受這次轉(zhuǎn)移围俘，從i轉(zhuǎn)移到j(luò)狀態(tài)琢融，反之則保持原i狀態(tài)。但是我們在實際應(yīng)用中使用這個方法時發(fā)現(xiàn)漾抬，很多情況下接受率普遍很低宿亡，導(dǎo)致馬氏鏈狀態(tài)轉(zhuǎn)移緩慢纳令，最終收斂的速度非常慢。為了解決這個問題坤按，我們還是采用類似等式（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)多，且維度比較高時骇扇，MH算法的接受率仍然差強人意。要是每次都接受該多好啊少孝。那什么樣的情況下，我從i到j(luò)時稍走，每次都能接受呢？（即接受率為1）粱胜。最終發(fā)現(xiàn)狐树，我們每次可以沿著垂直于某個變量維度的軸走。即通過迭代的方法冗恨，每一次只對一個變量進行采樣味赃。舉一個二維空間的例子，假設(shè)一個概率分布p(x,y)傲武，來看x坐標(biāo)相同的兩個點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)
下圖給出了一個更直觀的表示：
即喘垂，從A到B和從B到A的轉(zhuǎn)移是直接滿足細致平穩(wěn)條件的。因此我們不需要等式(6)中的接受率來幫忙得院，即接受率為1.圖中假設(shè)初始狀態(tài)為A章贞，則從A到下一個概率轉(zhuǎn)移矩陣分別為：
Q(A?>B)=p(yB|x1)
Q(A?>C)=p(xC|y1)
Q(A?>D)=0(12)
因此類似于曼哈頓距離的方法，狀態(tài)轉(zhuǎn)移總是沿著橫平豎直的街區(qū)進行蜕径。這邊是Gibbs抽樣算法的核心思想。下圖給出了一個Gibbs抽樣的直觀圖丧荐。
3 收斂條件的判斷
我們都知道當(dāng)概率狀態(tài)轉(zhuǎn)移穩(wěn)定時虹统，其分布便是所要求的聯(lián)合概率分布。但我們不可能通過如等式(2),(3)的方法來每轉(zhuǎn)換一步就求其概率分布渡冻，比較是否改變忧便。主要原因有二，其一是不可把所有變量間的轉(zhuǎn)移概率都找到珠增，其二矩陣計算耗時耗力蒂教。常見的方法便是通過burn-in的方法，多跑幾次凝垛。也有通過計算當(dāng)前狀態(tài)下的聯(lián)合分布可能性函數(shù)，然后根據(jù)Autocorrelation Function(ACF)的變化速率來判斷迭代是否收斂炭分。
So long, and thanks for all the fish.
參考
[1]PRML讀書會第十一章 Sampling Methods
[2]LDA-math-MCMC 和 Gibbs Sampling
[5]What-are-Markov-Chain-Monte-Carlo-methods-in-laymans-terms
