侵刪

轉(zhuǎn)： 從 bayes 到mrbayes

一、Bayes inference

可以簡單而不負(fù)責(zé)任的講，貝葉斯學(xué)派起源于一個(gè)公式：

圖片來自http://www.ruanyifeng.com/blog/2011/08/bayesian_inference_part_one.html

這個(gè)公式當(dāng)年被Thomas Bayes先生提出來用以解決“逆概”問題（ Inverse probability，現(xiàn)在不稱“逆概”，而叫做統(tǒng)計(jì)推斷，inferential statistics）。

舉個(gè)例子，如果告訴一個(gè)黑色袋子里有N個(gè)紅球和M個(gè)籃球，求每次取出紅球的概率，這種問題就叫“順概”問題。而如果依次告訴前X次抽出球的顏色，求袋子里紅球的頻率，則這種問題為“逆概”問題。

事實(shí)上，貝葉斯學(xué)派的對(duì)頭，概率學(xué)派，也有其解決“逆概”問題的方法，即從樣本估計(jì)總體 (frequentist inference)。簡單來講，如果前十次抽出的球中，有三個(gè)為紅色，七個(gè)為藍(lán)色，則認(rèn)為袋子中紅球與藍(lán)旗的比例為3:7（即紅球頻率0.3）。這屬于“歸納（inductive）”的思維方法。

不同于概率學(xué)派，貝葉斯采取的是“演繹（deductive）”的思維方式。先假設(shè)袋子里紅球的頻率為0.5，稱為先驗(yàn)概率p(A)。如果第一次取出的是紅球，求后驗(yàn)概率p(A|B)（即第一次取出的是紅球的情況下，這個(gè)袋子里紅球的頻率）?？梢姡惾~斯方法可以通過這樣不斷地取樣，繼而利用調(diào)整因子p(B|A)/p(B)不斷地調(diào)整后驗(yàn)概率，當(dāng)后驗(yàn)概率趨于穩(wěn)定時(shí)，則可用以估算袋子里紅球的頻率。

這個(gè)公式十分偉大，然而其所蘊(yùn)含的思想?yún)s是十分簡單而普通，亦為我們常人所頻繁使用。舉例來說，當(dāng)我沒出門的時(shí)候，推斷天會(huì)下雨的概率為0.5（p(A)），即先驗(yàn)概率。而我出門后，發(fā)現(xiàn)天空中有積雨云，那么，利用這一與下雨相關(guān)率極高的事件(B)，我則可以把天會(huì)下雨的概率上升到0.9(P(A|B))，即后驗(yàn)概率。維基上的一句話對(duì)此公式極富解釋力，如下（稍有修改）：

“the bayesian approach combines the prior probability of a P(A) with the likelihood of the data (B) to produce a posterior probability distribution on P(A|B).”

然而這個(gè)公式在它被提出后的很長時(shí)間不能得到發(fā)展，究其原因在于，在統(tǒng)計(jì)學(xué)家眼里，世界上所有事物都是概率分布，而非確定狀態(tài)。因此，貝葉斯公式通常面對(duì)的并不是一個(gè)概率到另一個(gè)概率的問題，而是一個(gè)概率分布到另一個(gè)概率分布，因此十分耗費(fèi)計(jì)算能力，受限于時(shí)代的計(jì)算能力，貝葉斯公式一直在世界的角落徘徊，直到MCMC方法的出現(xiàn)。

二、MC與MC

兩個(gè)MC，是完全不同的東西。前一個(gè)是指Markov Chain, 馬爾科夫鏈，后一個(gè)為Monte Carlo，蒙特卡洛。先說后面那個(gè)。

Monte Carlo， 蒙特卡洛，原指摩納哥大公國一著名賭場。

而賭博總是與概率 有關(guān)，因此蒙特卡洛在此指的是Monte Carlo 方法，隨機(jī)模擬，即利用計(jì)算機(jī)進(jìn)行隨機(jī)抽樣，產(chǎn)生大的數(shù)據(jù)結(jié)果，用以探究目標(biāo)的概率分布概況，或是模擬一個(gè)概率分布。

隨機(jī)模擬小時(shí)候玩過，拿計(jì)算機(jī)隨機(jī)產(chǎn)生一個(gè)數(shù)。那個(gè)時(shí)候沒有概念，以為隨機(jī)的意思就概率相同。事實(shí)上，隨機(jī)過程中蘊(yùn)含的概率分布是可以有偏倚的。隨機(jī)模擬中有一個(gè)重要問題就是模擬特定的概率分布，最簡單的算是均勻分布 Uniform(0,1)，可通過線性同余發(fā)生器。而其他一些著名分布，例如正態(tài)分布，指數(shù)分布，Gamma分布，t分布，都可以通過這個(gè)均勻分布轉(zhuǎn)換得到。然而很多時(shí)候，現(xiàn)實(shí)中面對(duì)的總是一些十分復(fù)雜的分布，以至于蒙特卡洛方法長時(shí)間里得不到發(fā)展，直到前一個(gè)MC，即馬爾科夫鏈的提出，承當(dāng)了概率分布間的轉(zhuǎn)換器。

Markov Chain，馬爾科夫鏈，又稱馬氏鏈，指的是狀態(tài)間轉(zhuǎn)換的隨機(jī)過程。http://www.52nlp.cn/lda-math-mcmc-%E5%92%8C-gibbs-sampling1中舉了一個(gè)很好的例子。如果將社會(huì)階層分為“上”，“中”，“下”三個(gè)，分別以1,2,3代表，則某個(gè)人社會(huì)階層的代代傳遞軌跡則形成了一條馬爾科夫鏈，比如說：1->1->2->3->2->1......

馬爾科夫鏈有三個(gè)重要特質(zhì)。

其一為“memorylessness”，即鏈的下一個(gè)狀態(tài)只依賴于現(xiàn)在的狀態(tài)，而與之前的狀態(tài)都無關(guān)。這樣就能很好的把貝葉斯公式融合進(jìn)去，貝葉斯公式干的事情就是基于此刻狀態(tài)，預(yù)測下一狀態(tài)。

其二，馬爾科夫鏈有概率轉(zhuǎn)換矩陣P。如上面“社會(huì)階層”的例子，則有P：

其三，馬爾科夫鏈具有收斂性，且不管初始狀態(tài)如何，其收斂后的結(jié)果相似。此話該如何理解？還舉上面“社會(huì)階層”的例子，假設(shè)初代社會(huì)的階層分布已知，π0=[0.21,0.68,0.11]，則我們可以計(jì)算前n代人的分布狀況如下：

可見從第七代開始，社會(huì)階層的分布則開始趨于穩(wěn)定，即所謂的收斂。而且，不管初始分布如何變，其收斂后的分布（即平穩(wěn)分布）總是相同的。

三、MCMC (Markov Chain Monte Carlo)

綜上，我認(rèn)為MCMC可以理解為一種“鏈?zhǔn)饺印钡倪^程。

總體的無限性或不可測量性是世界之常態(tài)，因此，人類認(rèn)識(shí)世界的一個(gè)主要方式就是取樣。

最簡單的取樣應(yīng)該就是隨機(jī)取樣了，袋子里有100個(gè)球，抽取10個(gè)球作為樣本，推測袋子里的情況（統(tǒng)計(jì)推斷）。我以前一提起抽樣，就以為是隨機(jī)抽樣，現(xiàn)在才發(fā)現(xiàn)，原來抽樣是個(gè)大學(xué)問（日后若是有機(jī)會(huì)，也想給自己總結(jié)一下”抽樣“這個(gè)命題）。然而人類對(duì)總體的探求不總是想知道總體的整體狀態(tài)，還有一種情況，比如說黑屋子里有1,000,000個(gè)人，想知道其中最高的有多高。對(duì)于這種情況，如果還是用隨機(jī)抽樣的話，效率則會(huì)顯得捉襟見肘了。此時(shí)，則有其他的抽樣方法可以發(fā)揮重要作用，MCMC就是其中一例。

這種鏈?zhǔn)饺用砍橐粋€(gè)人就會(huì)評(píng)估這個(gè)人的身高是最高的概率（或然率），并以此為先驗(yàn)概率，接著抽取下一個(gè)人，利用調(diào)節(jié)因子（肯定跟這個(gè)人身高有關(guān)）求得后驗(yàn)概率，即后來抽取的這個(gè)人是這群人中身高最高的一個(gè)的概率。如果后驗(yàn)概率比先驗(yàn)概率高，則保留這次取樣，若比先驗(yàn)概率低，則放棄這次取樣，接著再循環(huán)這一步驟。

由于最高的身高是不會(huì)變的（收斂性？），因此這種趨向性的鏈?zhǔn)饺釉诮鉀Q這種問題時(shí)，總是會(huì)無限逼近這一最高值，顯然比隨機(jī)抽樣的方式效率高出很多。

三、Metropolis-Hastings MCMC

是一種十分常見的MCMC方法。起源于Metropolis算法，即

1.隨機(jī)定起始狀態(tài)i，計(jì)算其概率f(i)。

2.選擇一臨近狀態(tài)j，計(jì)算其概率f(j)。

3.如果f(j)/f(i) >=1，則接受狀態(tài)j，替換掉原狀態(tài)i。

如果f(j)/f(i) <1，則利用均勻分布 Uniform(0,1)產(chǎn)生一隨機(jī)數(shù)，

如果該隨機(jī)數(shù)小于f(j)/f(i)，則接接受狀態(tài)j，替換掉原狀態(tài)i。

如果該隨機(jī)數(shù)大于f(j)/f(i)，則保持原狀態(tài)j不變。

4.重復(fù)步驟2、3， 直到它到達(dá)一個(gè)均衡分布（蛇游到山頂不愿下來了）。

此算法有個(gè)前提就是在i狀態(tài)提議j狀態(tài)的概率和在j狀態(tài)提議i狀態(tài)的概率相等。如果不等，則利用Hastings corrections，即成Metropolis-Hastings MCMC。

可見，這個(gè)算法對(duì)馬爾科夫鏈的抽樣進(jìn)行了一些限制，讓鏈的方向總是走向概率較高的狀態(tài)，一旦到頂，則盤踞其上，無法下來，從而探測平穩(wěn)分布（又稱后驗(yàn)概率分布）中概率最高值對(duì)應(yīng)的狀態(tài)。

四、Metropolis-coupled MCMC（MC）

如果只是一個(gè)山坡，想要到達(dá)山頂，則利用第三節(jié)中的Metropolis-Hastings MCMC方法會(huì)十分有效。然而現(xiàn)實(shí)情況下，大多后驗(yàn)概率分布更像是有很多小山坡的地貌。蛇游到其中一個(gè)山坡的頂上則不愿下來了，因此無法確定后驗(yàn)概率分布的最高峰。為了解決這個(gè)問題，則有了這里的Metropolis-coupled MCMC算法。

如果把Metropolis-Hastings MCMC算法比作一條蛇的話，那么它就是一條近視的蛇，一旦到達(dá)某一坡頂，就算附近有一座更高的峰它也看不見，因此不愿下來。而Metropolis-coupled MCMC 算法除了放了一條這樣的近視蛇之外（cold chain），還放n條目光十分跳躍的蛇（heated chains），以實(shí)現(xiàn)峰頂與峰頂之間的跳躍。除此之外，這個(gè)算法每個(gè)迭代后還會(huì)有一個(gè)交換的過程，利用Metropolis式的比較法隨機(jī)交換兩條蛇的位置。

在算法的最后，只有cold chain的軌跡被輸出來。

四、MrBayes

MrBayes 是一套生物信息軟件，用以系統(tǒng)發(fā)育的貝葉斯推斷。其核心算法就是Metropolis-coupled MCMC（MC）。簡單來講，就是采取鏈?zhǔn)饺拥姆绞?，在天文?shù)字的系統(tǒng)樹中求得似然率最高的那一棵樹。

運(yùn)行貝葉斯最簡單命令行如下：

mb

exe <filename>

prset aamodelpr=fixed(jones)

lset rates=invgamma

mcmcp samplefreq=500 printfreq=500 diagnfreq=500

mcmc ngen=2000000

sump <filename> burnin=2000

sumt <filename> burnin=2000

解釋如下：

mb，打開MrBayes的運(yùn)行界面。

exe，輸入aligment后的序列文件，nexus格式。

prset 和 lset用來設(shè)定進(jìn)化模型。我一般會(huì)先用protest篩選進(jìn)化模型。

mcmcp，設(shè)定后續(xù)mcmc算法的參數(shù)。下設(shè)一下參數(shù)：

samplefreq，取樣頻率（每算多少代取一次樣）；

printfreq，展示迭代結(jié)果頻率（每算多少代展示一次結(jié)果）；

diagnfreq，診斷頻率，即隔多少迭代診斷一次馬爾科夫鏈?zhǔn)欠襁_(dá)到平穩(wěn)分布。

mcmc，開始跑mcmc算法，ngen參數(shù)，設(shè)定總共的迭代代數(shù)。

mcmc的運(yùn)行需要一定時(shí)間，會(huì)產(chǎn)出如下文件：

該mcmc的運(yùn)行共分兩個(gè)run，默認(rèn)來說，每個(gè)run共設(shè)四條馬爾科夫取樣鏈，一冷，三熱。當(dāng)兩個(gè)run的冷鏈結(jié)果開始穩(wěn)定聚合，則可認(rèn)為馬爾科夫鏈開始收斂（我相信diagnfreq中的diagn就是利用此）。

最后，sump 和 sumt 分別用來統(tǒng)計(jì)收斂后的冷鏈取樣結(jié)果。burnin，即舍去的樣本。就本例來說，共跑2，000,000代，500代一取樣，則共取4000個(gè)樣，舍去2000個(gè)，則剩余2000個(gè)樣本用來擬合目標(biāo)分布。

PS: burn in原指產(chǎn)品工藝中的預(yù)燒實(shí)驗(yàn)，即一批樣預(yù)燒直至商品損壞，從而判斷這批商品的質(zhì)量。比如說聲場一批皮鞋，要知道其鞋底質(zhì)量，則會(huì)從中抽樣，用機(jī)器反復(fù)折疊，模擬人走路時(shí)的鞋底折疊，直至鞋底損壞，則可估計(jì)該批鞋的質(zhì)量?？梢?，用來burn in的樣本最終都是要損毀的，因此在這里取“舍去樣本的個(gè)數(shù)”之意。

總結(jié)來說：

貝葉斯是利用后驗(yàn)概率作為依據(jù)的一種統(tǒng)計(jì)推斷方法

而MCMC是一種在巨大樣本中采取鏈?zhǔn)降摹叭臃椒ā保瑑烧卟豢苫煜?/p>

當(dāng)求后驗(yàn)概率的公式中出現(xiàn) #連續(xù)函數(shù)的積分，或者#巨大樣本的求和的情況時(shí)，這樣的后驗(yàn)概率直接是解不出來的，這個(gè)時(shí)候，就用MCMC這種方法，從大樣本中取樣，依據(jù)取樣結(jié)果來解后驗(yàn)概率的值。