侵删
一、Bayes inference
可以简单而不负责任的讲,贝叶斯学派起源于一个公式:
图片来自http://www.ruanyifeng.com/blog/2011/08/bayesian_inference_part_one.html
这个公式当年被Thomas Bayes先生提出来用以解决“逆概”问题( Inverse probability,现在不称“逆概”,而叫做统计推断,inferential statistics)。
举个例子,如果告诉一个黑色袋子里有N个红球和M个篮球,求每次取出红球的概率,这种问题就叫“顺概”问题。而如果依次告诉前X次抽出球的颜色,求袋子里红球的频率,则这种问题为“逆概”问题。
事实上,贝叶斯学派的对头,概率学派,也有其解决“逆概”问题的方法,即从样本估计总体 (frequentist inference)。简单来讲,如果前十次抽出的球中,有三个为红色,七个为蓝色,则认为袋子中红球与蓝旗的比例为3:7(即红球频率0.3)。这属于“归纳(inductive)”的思维方法。
不同于概率学派,贝叶斯采取的是“演绎(deductive)”的思维方式。先假设袋子里红球的频率为0.5,称为先验概率p(A)。如果第一次取出的是红球,求后验概率p(A|B)(即第一次取出的是红球的情况下,这个袋子里红球的频率)。可见,贝叶斯方法可以通过这样不断地取样,继而利用调整因子p(B|A)/p(B)不断地调整后验概率,当后验概率趋于稳定时,则可用以估算袋子里红球的频率。
这个公式十分伟大,然而其所蕴含的思想却是十分简单而普通,亦为我们常人所频繁使用。举例来说,当我没出门的时候,推断天会下雨的概率为0.5(p(A)),即先验概率。而我出门后,发现天空中有积雨云,那么,利用这一与下雨相关率极高的事件(B),我则可以把天会下雨的概率上升到0.9(P(A|B)),即后验概率。维基上的一句话对此公式极富解释力,如下(稍有修改):
“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).”
然而这个公式在它被提出后的很长时间不能得到发展,究其原因在于,在统计学家眼里,世界上所有事物都是概率分布,而非确定状态。因此,贝叶斯公式通常面对的并不是一个概率到另一个概率的问题,而是一个概率分布到另一个概率分布,因此十分耗费计算能力,受限于时代的计算能力,贝叶斯公式一直在世界的角落徘徊,直到MCMC方法的出现。
二、MC与MC
两个MC,是完全不同的东西。前一个是指Markov Chain, 马尔科夫链,后一个为Monte Carlo,蒙特卡洛。先说后面那个。
Monte Carlo, 蒙特卡洛,原指摩纳哥大公国一著名赌场。
而赌博总是与概率 有关,因此蒙特卡洛在此指的是Monte Carlo 方法,随机模拟,即利用计算机进行随机抽样,产生大的数据结果,用以探究目标的概率分布概况,或是模拟一个概率分布。
随机模拟小时候玩过,拿计算机随机产生一个数。那个时候没有概念,以为随机的意思就概率相同。事实上,随机过程中蕴含的概率分布是可以有偏倚的。随机模拟中有一个重要问题就是模拟特定的概率分布,最简单的算是均匀分布 Uniform(0,1),可通过线性同余发生器。而其他一些著名分布,例如正态分布,指数分布,Gamma分布,t分布,都可以通过这个均匀分布转换得到。然而很多时候,现实中面对的总是一些十分复杂的分布,以至于蒙特卡洛方法长时间里得不到发展,直到前一个MC,即马尔科夫链的提出,承当了概率分布间的转换器。
Markov Chain,马尔科夫链,又称马氏链,指的是状态间转换的随机过程。http://www.52nlp.cn/lda-math-mcmc-%E5%92%8C-gibbs-sampling1中举了一个很好的例子。如果将社会阶层分为“上”,“中”,“下”三个,分别以1,2,3代表,则某个人社会阶层的代代传递轨迹则形成了一条马尔科夫链,比如说:1->1->2->3->2->1......
马尔科夫链有三个重要特质。
其一为“memorylessness”,即链的下一个状态只依赖于现在的状态,而与之前的状态都无关。这样就能很好的把贝叶斯公式融合进去,贝叶斯公式干的事情就是基于此刻状态,预测下一状态。
其二,马尔科夫链有概率转换矩阵P。如上面“社会阶层”的例子,则有P:
其三,马尔科夫链具有收敛性,且不管初始状态如何,其收敛后的结果相似。此话该如何理解?还举上面“社会阶层”的例子,假设初代社会的阶层分布已知,π0=[0.21,0.68,0.11],则我们可以计算前n代人的分布状况如下:
可见从第七代开始,社会阶层的分布则开始趋于稳定,即所谓的收敛。而且,不管初始分布如何变,其收敛后的分布(即平稳分布)总是相同的。
三、MCMC (Markov Chain Monte Carlo)
综上,我认为MCMC可以理解为一种“链式取样”的过程。
总体的无限性或不可测量性是世界之常态,因此,人类认识世界的一个主要方式就是取样。
最简单的取样应该就是随机取样了,袋子里有100个球,抽取10个球作为样本,推测袋子里的情况(统计推断)。我以前一提起抽样,就以为是随机抽样,现在才发现,原来抽样是个大学问(日后若是有机会,也想给自己总结一下”抽样“这个命题)。然而人类对总体的探求不总是想知道总体的整体状态,还有一种情况,比如说黑屋子里有1,000,000个人,想知道其中最高的有多高。对于这种情况,如果还是用随机抽样的话,效率则会显得捉襟见肘了。此时,则有其他的抽样方法可以发挥重要作用,MCMC就是其中一例。
这种链式取样每抽一个人就会评估这个人的身高是最高的概率(或然率),并以此为先验概率,接着抽取下一个人,利用调节因子(肯定跟这个人身高有关)求得后验概率,即后来抽取的这个人是这群人中身高最高的一个的概率。如果后验概率比先验概率高,则保留这次取样,若比先验概率低,则放弃这次取样,接着再循环这一步骤。
由于最高的身高是不会变的(收敛性?),因此这种趋向性的链式取样在解决这种问题时,总是会无限逼近这一最高值,显然比随机抽样的方式效率高出很多。
三、Metropolis-Hastings MCMC
是一种十分常见的MCMC方法。起源于Metropolis算法,即
1.随机定起始状态i,计算其概率f(i)。
2.选择一临近状态j,计算其概率f(j)。
3.如果f(j)/f(i) >=1,则接受状态j,替换掉原状态i。
如果f(j)/f(i) <1,则利用均匀分布 Uniform(0,1)产生一随机数,
如果该随机数小于f(j)/f(i),则接接受状态j,替换掉原状态i。
如果该随机数大于f(j)/f(i),则保持原状态j不变。
4.重复步骤2、3, 直到它到达一个均衡分布(蛇游到山顶不愿下来了)。
此算法有个前提就是在i状态提议j状态的概率和在j状态提议i状态的概率相等。如果不等,则利用Hastings corrections,即成Metropolis-Hastings MCMC。
可见,这个算法对马尔科夫链的抽样进行了一些限制,让链的方向总是走向概率较高的状态,一旦到顶,则盘踞其上,无法下来,从而探测平稳分布(又称后验概率分布)中概率最高值对应的状态。
四、Metropolis-coupled MCMC(MC)
如果只是一个山坡,想要到达山顶,则利用第三节中的Metropolis-Hastings MCMC方法会十分有效。然而现实情况下,大多后验概率分布更像是有很多小山坡的地貌。蛇游到其中一个山坡的顶上则不愿下来了,因此无法确定后验概率分布的最高峰。为了解决这个问题,则有了这里的Metropolis-coupled MCMC算法。
如果把Metropolis-Hastings MCMC算法比作一条蛇的话,那么它就是一条近视的蛇,一旦到达某一坡顶,就算附近有一座更高的峰它也看不见,因此不愿下来。而Metropolis-coupled MCMC 算法除了放了一条这样的近视蛇之外(cold chain),还放n条目光十分跳跃的蛇(heated chains),以实现峰顶与峰顶之间的跳跃。除此之外,这个算法每个迭代后还会有一个交换的过程,利用Metropolis式的比较法随机交换两条蛇的位置。
在算法的最后,只有cold chain的轨迹被输出来。
四、MrBayes
MrBayes 是一套生物信息软件,用以系统发育的贝叶斯推断。其核心算法就是Metropolis-coupled MCMC(MC)。简单来讲,就是采取链式取样的方式,在天文数字的系统树中求得似然率最高的那一棵树。
运行贝叶斯最简单命令行如下:
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的运行界面。
exe,输入aligment后的序列文件,nexus格式。
prset 和 lset用来设定进化模型。我一般会先用protest筛选进化模型。
mcmcp,设定后续mcmc算法的参数。下设一下参数:
samplefreq,取样频率(每算多少代取一次样);
printfreq,展示迭代结果频率(每算多少代展示一次结果);
diagnfreq,诊断频率,即隔多少迭代诊断一次马尔科夫链是否达到平稳分布。
mcmc,开始跑mcmc算法,ngen参数,设定总共的迭代代数。
mcmc的运行需要一定时间,会产出如下文件:
该mcmc的运行共分两个run,默认来说,每个run共设四条马尔科夫取样链,一冷,三热。当两个run的冷链结果开始稳定聚合,则可认为马尔科夫链开始收敛(我相信diagnfreq中的diagn就是利用此)。
最后,sump 和 sumt 分别用来统计收敛后的冷链取样结果。burnin,即舍去的样本。就本例来说,共跑2,000,000代,500代一取样,则共取4000个样,舍去2000个,则剩余2000个样本用来拟合目标分布。
PS: burn in原指产品工艺中的预烧实验,即一批样预烧直至商品损坏,从而判断这批商品的质量。比如说声场一批皮鞋,要知道其鞋底质量,则会从中抽样,用机器反复折叠,模拟人走路时的鞋底折叠,直至鞋底损坏,则可估计该批鞋的质量。可见,用来burn in的样本最终都是要损毁的,因此在这里取“舍去样本的个数”之意。
总结来说:
贝叶斯是利用后验概率作为依据的一种统计推断方法
而MCMC是一种在巨大样本中采取链式的“取样方法”,两者不可混淆。
当求后验概率的公式中出现 #连续函数的积分,或者#巨大样本的求和的情况时,这样的后验概率直接是解不出来的,这个时候,就用MCMC这种方法,从大样本中取样,依据取样结果来解后验概率的值。