2014年6月11日星期三

MCMC算法简介

Other interesting/ ML MAP Bayesian: https://engineering.purdue.edu/kak/Tutorials/
Gibbs Motif finding: http://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/gibbs.pdf

原文链接：百度贴吧  http://tieba.baidu.com/p/1681867134

       前天一个应用数学专业的朋友问我一个问题。这种类型的问题在工程上叫做求解反问题(inverse problem)。这种问题在他们看来是非常麻烦的。但是作为一个学统计的人而言，这种问题是再常见不过的了。也许是不同学科对问题的理解不同吧。我这里就引用这个例子介绍一下统计学的基本思想以及MCMC算法。
原始问题是这样的，有一个正问题，背景似乎是研究岩石断层什么的：

       幻灯片中的等式成立(由相关物理知识决定)，称为forward problem。其中theta和phi是可以控制的量，称为输入。左侧的量称为输出，因此是输入量的函数。问题是：有颜色的参数是未知的。
       我们需要通过一组输入和输出的数据来估计参数，这样的问题称为反问题。从统计学角度上看，这个问题属于参数估计问题。
----------------------------------------
1、统计学/贝叶斯统计
我们首先要用统计语言重新描述这个问题：
       假设我们有n次试验，对应于theta和phi的不同组合。我们将第i次试验表述为

       其中epsilon_i是第i次实验的随机误差。
       假定epsilon_i服从独立同分布的正态分布。不妨假设正态分布的均值为0，因为如果不为0均值也可以被截距项A吸收掉。设正态分布方差为sigma^2，未知。
       统计学中的一个核心概念称为似然函数（likelihood function），它是观测到的量的概率分布。这里我们观测到R，所以就要写出R的联合概率分布。注意到R_i也是独立的正态分布，因此联合概率密度为

       由于概率密度要满足积分为1的正则性条件。因此统计上常常省略密度前面的常数项，然后用正比于表示所缺常数由正则性条件确定。
       然后我们这里应用贝叶斯统计的观点，认为参数也是有一定分布的。这是因为在样本量有限的情况下参数没办法准确确定，因此用一个分布来描述是合理的。假设我们给所有参数赋予无信息先验，根据贝叶斯公式，就得到后验密度

       贝叶斯统计认为参数的一切信息都包含在后验分布当中。于是我们只要研究这个后验密度函数就可以了。
       然而不幸的是，这个后验密度有5维，而且还不知道前面的归一化常数，因此没办法直观把握。于是贝叶斯统计中有了一个重要的边缘化的思想：考虑每个参数的边缘分布，即将其他参数积分掉。尽管贝叶斯统计历史悠久，这个边缘化思想却是近30年来随着Markov Chain Monte Carlo (MCMC)的兴起才开始出现的。
       简单地说，就是通过计算形如：

的式子求出每个参数的边缘分布，然后统计推断就结束了。
________________________________
*边缘分布(对除theta外其他变量的积分，从而消除其他变量影响)
       在介绍MCMC算法之前。我们先讲一下计算这个边缘分布的难点。第一个难点是归一化常数不知道，第二个则是高维数值积分本身的困难。
       因为这两个困难，目前认为最有效的算法是Monte Carlo算法。算法的基本思想是，产生服从后验分布的样本，设想如果我们能够产生服从P(A,B_1,B_2,\phi,\sigma^2|R)的样本，那么它的第一维就服从P(A|R)，依此类推。如果我能产生N个样本，N足够大，我们就可以用样本的边缘分布近似总体的边缘分布，也就是说，取出样本的第一维，就可以得到A的近似后验分布。因此问题转化为了已知密度函数（缺归一化常数）的抽样问题。
_______________________________
2、MCMC算法
       MCMC不是一个算法，而是一类算法的总称。从数学上讲，算法的思想是产生一个Markov Chain，以目标分布为平稳分布。根据Markov Chain的理论，一个Markov Chain从任意初值出发，都会收敛到其平稳分布（如果存在）。注意这个收敛是依分布的。换句话说，如果我模拟了一条这样的马尔科夫链，去除掉前面一部分样本之后，就可以认为后面的样本来自于平稳分布。这样上面所述的抽样问题就解决了。
       MCMC与传统的独立抽样的关键区别在于：抽样的结果是相关的。理论上这对最终结果没有影响（数学上称为遍历理论）。然而实际上仍然有一系列的问题。其中最主要的是收敛诊断问题。和传统算法不同，MCMC的结果无法直接看出是否收敛。因此有的人认为使用MCMC需要慎重。但是根据我的经验，MCMC的表现基本还是可以相信的。
       2.1、Metropolis算法
       史上第一个MCMC算法诞生于1942年，现在称为Metropolis算法。
       这个算法首先要求一个状态空间上的对称转移函数。即设x,y是状态空间上任意两点，则从x转移到y的概率和从y转移到x的概率相等。
       Metropolis算法描述如下：
设x是上一次的状态，通过对称转移函数得到新状态x',然后抽取随机变量U服从(0,1)上的均与分布，如果
1）U<=q(x')/q(x),则markov chain转移到x';
2）U>q(x')/q(x),则markov chain保持不变，仍为x，
其中q是抽样的目标的分布。
容易证明这样构造的markov chain以q为平稳分布。
       这个算法有一个推广，称为Hasting's Generalization。如果转移函数不是对称的，只需要将q(x')/q(x)改为{q(x')p(x'->x)}/{q(x)p(x->x')}，其中p(x'->x)是从x'转移到x的概率。
       2.2、Gibbs抽样
       Metropolis是一个包打天下的算法。但是和其他一切适用范围广的算法一样，这个算法收敛是很慢的。统计中最常见的MCMC算法称为Gibbs抽样，它可以看做是一种特殊的Metropolis算法。
       Metropolis算法的主要问题，是算法效率严重依赖于转移函数的选取。如果选得不好，算法表现会很差。Gibbs抽样则不存在这样的问题，但是它要求被抽样的函数具有一定的形式。Gibbs抽样主要针对高维抽样问题。这种问题如果直接用Metropolis会变得很难。Gibbs抽样的想法是：既然n维抽样很难，我们就把他化成n个一维的抽样。
      这里为了简单起见，我们举一个二维的例子，更高的维度可以直接推广：
设目标分布为q(x,y),我们取初值(x_0,y_0)，Gibbs抽样一次迭代为
1)抽取x_1服从q(x|y_0)，
2)抽取y_1服从q(y|x_1)，
这样就得到了(x_1,y_1)。一直重复下去就得到需要的链。
       很显然这个算法要求两个条件分布q(x|y)和q(y|x)具有简单的形式，否则这两个一维抽样也是不可以实现的。当然如果一维抽样不能直接实现，还可以通过嵌入Metropolis算法，用一个接受/拒绝的过程代替直接从条件分布中抽样。
和Metropolis算法相比Gibbs抽样少了拒绝的步骤，这是一般情况下Gibbs抽样效果更好的原因。
——————————————
       在介绍了两种基本算法之后，我们回到原来的问题。由于这个问题是多维的，我们首先考虑Gibbs抽样。
 I、让我们来看看每个参数的条件分布是什么：
       1、首先，A在给定B_1,B_2,phi,sigma^2下的分布是什么呢？
        答案是正态分布，因为分布的形式为exp{-“A的二次函数”}，这是正态分布的核（回顾一下我们看分布的时候是不关心前面的常数系数的，所以看到这样的函数形式，就可以确定它是正态分布了）。事实上，(A,B_1,B_2)这个三维向量在给定phi和sigma^2的条件下服从联合正态分布（核的形式为exp{-“A,B_1,B_2的二次型”}）。经过配方（由矩阵运算可以迅速得到），可以算出这个正态分布的均值和协方差矩阵。
       2、然后，sigma^2在给定A,B_1,B_2,phi下是什么分布呢？
       我们注意到如果把sigma^2换成1/sigma^2就具有gamma分布的形式。所以这个分布被称为inverse gamma分布。
       3、最后我们来看phi在给定A,B_1,B_2,sigma^2下的分布。
       这里补充一下：原始物理问题有个约束条件，phi取值于0°到45°之间。非常不幸的是，phi的分布不属于任何已知的分布类型。因此对于phi的抽样我们只能借助metropolis算法来实现了。
当然这里要怎么选取转移函数也是一个可以讨论的问题。但是我考虑到phi的取值范围比较小，就使用一种最笨的抽样方式来解决了。
       我使用独立的[0°,45°]的均匀分布作为转移。也就是说，不管上一次的结果是什么，下一次的新值都从这个均匀分布中独立的产生。这个转移函数显然是对称的。然后用metropolis定义的概率进行接受或拒绝就行了。这个转移函数的缺点是不能利用上次成功的经验，但是优点是两个点之间相关性会减小（如果接受）。由于这个问题本身的困难程度不大，只要抽样的次数足够多（如抽10W次），就可以比较精确得得到后验分布了。
    II、再讲一下抽样完毕以后如何做分析：
       首先要砍掉前面没有达到平稳的点，称为burn-in。对于这种简单问题砍掉5000个点已经绰绰有余了。然后对于剩下的点，一维维的画直方图，或者做密度估计，都可以。一维问题总是很容易处理的，通过画图你一定会对你感兴趣的参数获得认识了。
END
--------------------------讨论-------------------
     1、请问 怎么证明这个算法的正确性 比如说 他的不可约 平稳分布和极限分布是否相同等等
       一个Markov chain收敛条件是不可约、非周期、正常返。这三个条件是容易成立的(只要模型可识别)。至于极限分布是不是想要的分布，只要验证目标分布是这个迭代算法的不动点就行了。对于大多数的问题这种验证都是平凡的，因此在统计文献中一般都省略此步骤。
      2、在运用了MCMC后，你文中最后得到的结果应该是待估计参数的边缘后验分布吧，但实际运用中需要估计参数的确切值，这个确切的值是多少？比如A1的值估计出来为多少？
       贝叶斯统计认为，参数的所有信息都包含在后验分布当中。同时，通过有限的样本是无法准确估计参数的，因此用后验分布来表示参数的不确定性。
       如果需要点估计，可以根据实际情况选取后验均值或众数（posterior mode）。但是这样的方法会忽略参数的不确定性。因此我建议在后续计算中选用可以保留这种不确定性的方法。最近应用数学中有个很火的领域叫uncertainty quantification，与此相关。
................
—————MCMC资源————原文链接：http://asc.2dark.org/node/18
google：MCMC tutorial
David MacKay's book(electronic version availiable):
http://www.inference.phy.cam.ac.uk/mackay/itila/
发帖者 时间：
标签：

没有评论:

发表评论

订阅： 博文评论 (Atom)