一文读懂MCMC算法(马尔科夫链蒙特卡洛)

理解本文需要一些贝叶斯基础,小白可移步https://www.jianshu.com/p/c942f8783b97

为了理解MCMC,我们依然是从一个具体的事例出发:假设当我们来到了一个小人国,我们感兴趣的是小人国的国民的身高分布情况,即有多少人高1米,有多少人高0.5米,又有多少人像我们正常人一样高。一种解决这个问题的暴力方法是找遍这个小人国的所有人,然后都测量身高,但显然,这是一个愚公移山式的方法,在很多情况下都是不可能的。

所以,由于精力有限,我们只找到了10个小人国的人民,这十个人的高度分别是:
1.14,1.02,1.08,0.96,0.79,0.94,1.00,0.93,1.13,1.02

聪明的我们的直觉是,这大概符合一个正态分布,然后我们碰到了一个开挂了的长老,他说:“没错,就是正态分布,而且标准差sd=0.1,现在让我看看你们这些愚蠢的人类能不能知道这个正态分布的平均值μ是多少!”。

正态分布Normal distribution

一对分别名为马尔科夫和蒙特卡洛的名侦探组合就此登场,他们说:“首先,我们先随便猜一个平均值μ,比如μ(1) = 0.8好了。”

*我这里用μ(n)表示第n个提出的μ值,所以μ(1)是提出的第一个μ值,μ(2)是第二个提出的μ值。

接着他们要做的事,是要确定另一个值μ(2),通常我们要谨慎一点去选择一个和之前提出来的值μ(1)差别不大的值,比如μ(2) = 0.7。

接着的问题就是,我们需要判断:究竟是μ(1) 更符合实际情况,还是μ(2) 更符合实际情况呢?但要如何作出这个判断呢,这里就要用到前面的贝叶斯公式了。

判断哪个值更好,实际上是在问,基于目前观察到的数据,得到这个参数μ的可能性哪个更大?即:
已知D = {1.14,1.02,1.08,0.96,0.79,0.94,1.00,0.93,1.13,1.02}
p(μ(2)|D) 大于还是小于还是等于 p(μ(1)|D) ?

这不就是在问谁的后验概率更大嘛?

为了解决这件事,一种思路是我们要把p(μ(2)|D) 和 p(μ(1)|D) 都用前面的贝叶斯公式解出来。但你很快就会发现在这种情况下证据概率p(D)会很难算。

但如果我们转念一想,我们要做的是比较p(μ(2)|D) 和p(μ(1)|D) ,那么我们其实只要求p(μ(2)|D) / p(μ(1)|D) 就可以了,如果这个比值大于1,则说明μ2的后验概率更大,更符合实际情况

而实际上,
p(μ(2)|D) / p(μ(1)|D)
= (p(μ(2))p(D|μ(2)) / p(D)) / (p(μ(1))p(D|μ(1)) / p(D))
= p(μ(2))p(D|μ(2)) / p(μ(1))p(D|μ(1))

可以看到,由于分子和分母上的P(D)被相互抵消了,剩下来需要知道的值就只剩下μ(1)和μ(2)的先验概率,以及分别在μ=μ(1)和μ=μ(2)时得到数据D的似然概率了。

现在,我们面临的问题要比之前简化了一些。但实际上我们还需要处理似然概率的计算和先验概率的问题。

先说说似然概率p(D|μ(2)) 和p(D|μ(1)),此时的似然概率是完全可以算出来的。因为我们已经假设了数据D符合的是正态分布模型了,且已知sd=0.1(前面大师说的),所以当我们假设μ=μ(1)时,就确定了一个平均值为μ(1)和标准差为0.1的正态分布,也就确定了这个正态分布的概率密度函数f(x),接着基于我们的数据D计算x = 1.14,1.02,1.08...等值的概率密度,再将这些值相乘,便得到了似然概率*。

** 可以这样理解这一似然概率的计算:如果我们此时假定的正态分布与数据实际对应的正态分布越接近,就自然可能有越多的数据落在高概率密度函数的区域(即分布的平均值附近),如此,作为概率密度函数的连乘的似然概率自然也会更高。相比之下,如果你现在确定的正态分布的平均值为1500,标准差为1,那么它在x = 1.14的概率密度(概率密度的具体数值不等于概率,但是两者的数学意味是接近的)就会高度趋近于0,将这样一个数作为因子去计算似然概率,似然概率也显然将会比较低。

说完了似然概率,就轮到先验概率p(μ(2))和p(μ(1))的问题了。先验概率要怎么去算呢?答案是不用算!我们自己来定。但是先验概率毫无疑问对MCMC算法是有影响的,就像我们之前从之前贝爷的故事里看到的那样,后验概率是受到先验概率影响的。之所以一枚90%击中率的硬币几乎不能预测一个人是否得病,是因为得这种病本身的先验概率就超级低。 一个你需要记住的简单事实是,我们设定的先验概率越是背离真实的情况,就需要越多的数据去将先验概率修正,让后验概率符合实际的情况。 从这个意义上说, 贝叶斯推理不是无中生有,而是要先对数据背后的结果有一个信念(belief)的基础上,根据所见的数据,不断地修正原本的信念,使之逐渐接近数据背后对应的真实情况(贝叶斯公式本身就带有学习、更新的意味,所以学界现在还有种说法是我们的大脑是贝叶斯式的)

当我们看到数据的时候,通过观察数据,我们最开始会猜想μ=1的概率比较高,因此我们可以假定μ的先验概率是服从平均值为1,sd为0.5的概率分布,有了这样的先验概率分布,我们就可以计算得到当μ=μ1,μ2时分别对应的概率密度了。

搞定了先验概率和似然概率,就可以计算前面的公式p(μ(2))p(D|μ(2)) / p(μ(1))p(D|μ(1))了。当这个比值大于等于1时,我们就接受μ(2),而如果这个比值小于1,我们就以这个比值为概率接受μ(2)。比如比值为0.5时,我们只有50%的概率接受μ(2)。当不接受的时候,我们得重新寻找一个μ(2),再进行同样的后验概率比较。

反复进行这样的步骤之后,我们可以想象,我们自然会更大程度地访问那些后验概率更高的μ值。我们访问不同的μ值的频率分布,就是关于μ值的后验概率分布(的近似)。至于这背后具体的数学推导,我们就不谈了。但注意,参数的近似后验分布并不是我们想要拟合的模型「即最开始的问题:小人国的国民的身高分布情况」。还记得我们最开始假设小人国的身高分布情况符合正态分布,且我们已经得知这个正态分布的标准差sd=1,而MCMC最终会告诉我们根据现有的数据,我们推断小人国身高分布的平均值μ,符合某个概率分布(比如平均值为1,sd为5),如果我们觉得合适,我们可以将μ的后验分布的平均值作为μ的最可能值。即,「小人国的国民的身高分布情况最有可能符合μ=1,sd=1的概率分布」。

最后总结一下MCMC算法:
(1)确定参数值的先验分布。
(2)先确定第一个访问(或者说采样)的参数的数值,作为当前参数数值
(3)根据当前访问的参数的数值,以一定的方式(比如Metropolis sampler)提出下一个待考虑访问的参数的数值。
(4)以比值的形式,比较当前参数数值和待考虑访问的参数数值的后验概率,计算后验概率涉及到先验概率和后验概率的概率密度。根据这个比值的大小,接受或拒绝该待考虑采样的参数数值。接受后则将该参数数值视为当前参数数值。
(5)重复(3)和(4),直到符合某种终止条件(比如说访问了10000个参数数值)

最终,将被采样的参数数值的频率分布作为对该参数的后验概率分布的近似。

看完以后,你是不是想说这么复杂的事,是人干的吗!?

废话,这种事当然是计算机来干啊,你还想手算不成?

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

推荐阅读更多精彩内容