MCMC及R实现

1贝叶斯集锦（3）：从MC、MC到MCMC2013-07-3123:03:39#####一份草稿贝叶斯计算基础一、从MC、MC到MCMC斯坦福统计学教授PersiDiaconis是一位传奇式的人物。Diaconis14岁就成了一名魔术师,为了看懂数学家Feller的概率论著作，24岁时进入大学读书。他向《科学美国人》投稿介绍他的洗牌方法，在《科学美国人》上常年开设数学游戏专栏的著名数学科普作家马丁•加德纳给他写了推荐信去哈佛大学，当时哈佛的统计学家Mosteller正在研究魔术，于是Diaconis成了Mosteller的学生。（对他这段传奇经历有兴趣的读者可以看一看统计学史话《女士品茶》）。下面要讲的这个故事，是Diaconis在他的文章TheMarkovChainMonteCarloRevolution中给出的破译犯人密码的例子。一天，一位研究犯罪心理学的心理医生来到斯坦福拜访Diaconis。他带来了一个囚犯所写的密码信息。他希望Diaconis帮助他把这个密码中的信息找出来。这个密码里的每个符号应该对应着某个字母，但是如何把这些字母准确地找出来呢？Diaconis和他的学生Marc采用了一种叫做MCMC（马尔科夫链蒙特卡洛）的方法解决了这个问题。这个MCMC方法就是这一节我们所要讨论的内容。(1)贝叶斯推断的计算问题在上节我们看到，贝叶斯统计学是利用后验分布对θ进行推断。这种推断的计算很多情况下要用积分计算来完成。比如，我们要计算θ的函数g(θ)的期望：E(g(θ∣x))=∫g(θ)fθ∣x(θ∣x)dθ其中函数f表示后验分布。当g(θ)=θ时，得到的就是关于θ的点估计。2但是对很多贝叶斯推断问题来说，有时候后验分布过于复杂，使得积分没有显示结果，数值方法也很难应用；有时候需要计算多重积分（比如后验分布是多元分布时）。这些都会带来计算上的很大困难。这也是在很长的时期内，贝叶斯统计得不到快速发展的一个原因。1990年代MCMC（MarkovChainMonteCarlo，马尔科夫链蒙特卡洛）计算方法引入到贝叶斯统计学之后，一举解决了这个计算的难题。可以说，近年来贝叶斯统计的蓬勃发展，特别是在各个学科的广泛应用和MCMC方法的使用有着极其密切的关系。（2）蒙特卡洛方法（MonteCarlo）蒙特卡洛方法是一种随机模拟方法，随机模拟的思想由来已久（参见下面的蒲丰投针的例子），但是由于难于取得随机数，随机模拟的方法一直发展缓慢。而蒙特卡洛方法的出现得益于现代电子计算机的诞生，在1944年由Metropolis和Ulam提出于二战时美国原子弹研究的曼哈顿工程之中。蒙特卡洛这个名字是由Metropolis起的，借用了那个著名的赌场的名字，因为赌博总是和概率相关。例：蒲丰投针1777年，法国学者蒲丰提出的一种用随机试验来求圆周率的方法set.seed(1234)l-0.8#针的长度n-1e+06#重复100000次.u1-runif(n)#取随机数x-1/2*u1#x是针中心到最近的线的距离u2-runif(n)y-l/2*sin(u2*2*pi)z-as.numeric(x=y)#相交的充要条件是x=y。pi.e-n*l/sum(z)#π的估计式pi.e##[1]3.141蒙特卡洛模拟的过程是随机的，但解决的不仅有随机的问题也可以有确定性的问题。蒙特卡洛可以视为一种思想的泛称，只要在解决问题分人过程中，利用大量的随机样本，然后对这些样本结果进行概率分析从而得到问题求解的方法，都可以称之为蒙特卡洛方法。蒙特卡洛方法的基础是（利用计算机）生成指定分布的随机数的抽样。在统计上发展了一些方法来解决这个问题。一般来说，从均匀分布U（0,1）的样本较容易生成，通过线性同余发生器可以生成伪随机数序列，这些序列的各种统计指标和均匀分布U(0,1)的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。常见的分布的随机数都可以基于均匀分布的样本生成。对定积分的计算是蒙特卡洛方法的一种重要的应用，而很多统计问题，比如计算概率、矩（期望、方差等）都可以归结为定积分的计算比如说我们想计算θ的期望，期望的计算公式为：E(θ)=∫θp(θ)dθ怎么样通过抽样的方法来解决这个问题呢？3在θ的分布中抽取独立样本θ1,θ2,...,θn，n足够大用样本均值θˉ=1n∑ni=1θi作为E(θ)当n趋近无穷时，利用大数定律，样本均值会收敛到E(θ)（这称为仿真一致性）从上面这个简单例子，蒙特卡洛方法的基本思路是：（a）针对实际问题建立一个简单易行的概率统计模型，使问题所求的解为该模型的概率分布或者数字特征，比如：某个事件的概率或者是某个随机变量的期望值。（b）对模型中的随机变量建立抽样方法，在计算机上进行模拟试验，得到足够的随机抽样，并对相关事件进行统计（c）对试验结果进行分析，给出所求解的估计及其精度（方差）的估计上面例子中所用的这种模拟方法是把有限区间上的积分看做这个区间上均匀分布随机变量的期望值，叫做逆分布函数法（inverse-CDFmethod，也称作逆变换法）。这种方法简便易行，但对于无穷区间的积分是不适用的。所以还有一些其它的抽样方法，比如importancesampling;rejectionsampling;slicesampler等，不再一一详述。对于复杂或者高维的分布，利用蒙特卡洛方法生成随机样本是比较困难的，因此需要发展一些更复杂的随机模拟技术。（3）马尔科夫链（MarkovChain）和MCMC在蒙特卡洛模拟中，我们在后验分布中抽取样本，当这些样本独立时，利用大数定律样本均值会收敛到期望值。如果得到的样本是不独立的，那么就要借助于马尔科夫链进行抽样。MCMC方法就是为了这个目的而诞生的。马尔科夫链又称为马尔科夫过程是一种离散的随机过程，用俄罗斯数学家安德烈.马尔科夫的名字命名。随机过程可以看做一个随时间变化的随机变量序列，在物理上用来表示在某个空间中物体的位置随机变化的轨迹；在贝叶斯统计中，物体的轨迹可以看做蒙特卡洛算法的结果，而“空间”就是支持后验分布p(θ∣x)的样本空间。马尔科夫链的定义如下，设θ(t)是一个随机过程，如果它满足下面的性质：p(θ(t+h)=xt+h∣θ(s)=xs,s≤t)=p(θ(t+h)=xt+h∣θ(t)=xt,∨h0这个定义又称为马尔科夫性质，对一个马尔科夫链来说，未来状态只与当前t时刻有关，而与t时刻之前的历史状态无关（条件独立）。在时刻t，一个遍历的马尔科夫链从状态i可以向链的任一状态（包括状态i自身）转移。从状态i到状态j的转移概率记作p(i,j)或者p(i\rightarrowj),p(j\midi)。所有这些转移概率构成了状态转移矩阵，记为P。4马尔科夫链的一个很重要的性质是平稳分布。简单的说，主要统计性质不随时间而变的马尔科夫链就可以认为是平稳的。数学上有马氏链收敛定理，当步长n足够大时，一个非周期且任意状态联通的马氏链可以收敛到一个平稳分布π(x)。这个定理就是所有的MCMC方法的理论基础。我们对一个马氏链进行状态转移。设从初始分布π0出发，假设到第n步这个链可以收敛到平稳分布π0，那个这个过程可以表示为：X0π0(x),X1π1(x),...,Xnπn(x)=π(x),Xn+1π(x),Xn+2π(x),...在利用马氏链进行抽样时，在收敛之前的一段时间，比如上面的前n-1次迭代，各个状态的边际分布还不能认为是稳定分布，所以在进行估计的时候，应该把前面的这n-1个迭代值去掉。这个过程称为“burn-in”。MCMC方法就是构造合适的马尔科夫链进行抽样而使用蒙特卡洛方法进行积分计算。既然马尔科夫链可以收敛到平稳分布。我们可以建立一个以π为平稳分布的马尔科夫链，对这个链运行足够长时间之后，可以达到平稳状态。此时马尔科夫链的值就相当于在分布π(x)中抽取样本。利用马尔科夫链进行随机模拟的方法就是MCMC。这种方法最早是由Metropolis研究物理学中粒子系统的平稳性质想到的，他得到了第一个MCMC的算法：Metropolis算法。由此，一系列的MCMC方法相继诞生。二、M-H算法MCMC方法最早由Metropolis（1954）给出，后来Metropolis的算法由Hastings改进，合称为M-H算法。M-H算法是MCMC的基础方法。由M-H算法演化出了许多新的抽样方法，包括目前在MCMC中最常用的Gibbs抽样也可以看做M-H算法的一个特例。我们来介绍基本的M-H算法。马氏链的收敛性质主要由转移矩阵P决定,所以基于马氏链做抽样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布π(x)这里用到了马尔科夫链的另一个性质，如果具有转移矩阵P和分布π(x)的马氏链对所有的状态i,j满足下面的等式：π(i)p(i,j)=π(j)p(j,i)这个等式称为细致平衡方程。满足细致平衡方程的分布π(x)是平稳的。所以我们希望抽样的马尔科夫链是平稳的，可以把细致平衡方程作为出发点。但是一般情况下，任给的分布分布π(x)不一定是平稳的。那么为了让细致平衡方程成立，我们可以引入一个函数α(i,j)满足≤α(i,j)≤1,使得π(i)p(i,j)α(i,j)=π(j)p(j,i)α(j,i)按照对称性，我们可以取：α(i,j)=π(j)p(j,i)α(j,i)=π(j)p(i,j)记5q(i,j)=p(i,j)α(i,ij)q(j,i)=p(j,i)α(j,i)这样就得到了新的以q(i,j)为转移概率的马尔科夫链，且具有平稳分布π(x)这里的α(i,j)叫做接受率，也就是在原来的马氏链上以概率α(i,j)接受从状态i到状态j的转移概率为p(i,j)的状态转移。进一步的，如果接受率过小，会导致马氏链过多地拒绝状态转移，这样马氏链的收敛速度会很慢。因此，我们可以考虑放大α(i,j)，使π(i)p(i,j)α(i,j)=π(j)p(j,i)α(j,i)中的两个方程中大的\alpha$取到1，小的同比例放大。这样拒绝率就可以表示为：α(i,j)=min(1,π(j)p(j,i)/π(i)p(i,j))这就是M-H算法。我们的讨论都是针对离散状态的，对连续状态的马尔科夫链依然有相同的结果。在连续状态中表示状态转移概率的项用条件概率密度代替，称为转移核。M-H算法的步骤：(1)构造合适的提议分布（Proposaldistribution)g(⋅∣Xt)在分布g中产生X0；(2)迭代下面的步骤：在g(⋅∣Xt)中生成新样本Y；从均匀分布U(0,1)抽取随机数U；如果U满足U≤f(Y)g(Xt∣Y)/f(Xt)g(Y∣Xt),则令Xt+1=Y（转移到新状态），否则Xt+1=Xt（状态不变）。其中f是目标分布，也就是我们需要进行抽样的后验分布；增加t值，进行下一步迭代。例：对自由度为v的t分布进行抽样，t分布的密度值由R函数dt来计算N-3000#抽样个数v-4#自由度x-c()x0-0#初始值x[1]-x0k-0#k表示拒绝转移的次数u-runif(N)#抽取均匀分布随机数for(iin2:N){y-rnorm(1,x[i-1],2)if(u[i]=(dt(y,v)/dt(x[i-1],v)))x[i]-yelse{x[i]-x[i-1]k-k+1}}k##[1]13616回到我们前面的故事。因为在文本中前后字母有依赖的关系，所以Marc利用前后字母的依赖关系，来描述整个密码文本中出现这些字符的概率模型。他用《战争与和平》作为标准的文本，统计一个字母到另一个字母的一步转移概率。然后他利用了M-H算法，在假设所有对应关系出现的可能性相等的前提下(也就是无信息先验的情况)，首先随机给出了字符和字母的对应关系。利用前边得到的转移概率，计算这种对应关系出现的概率p1。随机抽取两个字符，互换它们的对应关系，此时对应概率p2。如果P2P1，接受新的对应关系；否则，抛一枚以P2/P1的概率出现正面的硬币，如果出现正面，则接受新的对应关系，否