Mrbayes

侵删

转: 从 bayes 到mrbayes

一、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这种方法,从大样本中取样,依据取样结果来解后验概率的值。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,378评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,356评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,702评论 0 342
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,259评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,263评论 5 371
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,036评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,349评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,979评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,469评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,938评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,059评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,703评论 4 323
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,257评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,262评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,485评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,501评论 2 354
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,792评论 2 345

推荐阅读更多精彩内容