理解本文需要一些贝叶斯基础,小白可移步//www.greatytc.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,现在让我看看你们这些愚蠢的人类能不能知道这个正态分布的平均值μ是多少!”。
一对分别名为马尔科夫和蒙特卡洛的名侦探组合就此登场,他们说:“首先,我们先随便猜一个平均值μ,比如μ(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个参数数值)
最终,将被采样的参数数值的频率分布作为对该参数的后验概率分布的近似。
看完以后,你是不是想说这么复杂的事,是人干的吗!?
废话,这种事当然是计算机来干啊,你还想手算不成?