MCMC及R实现汇编.docVIP

贝叶斯集锦（3）：从MC、MC到MCMC 2013-07-31 23:03:39 #####一份草稿贝叶斯计算基础一、从MC、MC到MCMC斯坦福统计学教授Persi Diaconis是一位传奇式的人物。Diaconis14岁就成了一名魔术师,为了看懂数学家Feller的概率论著作，24岁时进入大学读书。他向《科学美国人》投稿介绍他的洗牌方法，在《科学美国人》上常年开设数学游戏专栏的著名数学科普作家马丁?加德纳给他写了推荐信去哈佛大学，当时哈佛的统计学家Mosteller 正在研究魔术，于是Diaconis成了Mosteller的学生。（对他这段传奇经历有兴趣的读者可以看一看统计学史话《女士品茶》）。 下面要讲的这个故事，是Diaconis 在他的文章The Markov Chain Monte Carlo Revolution中给出的破译犯人密码的例子。 一天，一位研究犯罪心理学的心理医生来到斯坦福拜访Diaconis。他带来了一个囚犯所写的密码信息。他希望Diaconis帮助他把这个密码中的信息找出来。 这个密码里的每个符号应该对应着某个字母，但是如何把这些字母准确地找出来呢？Diaconis和他的学生Marc采用了一种叫做MCMC（马尔科夫链蒙特卡洛）的方法解决了这个问题。这个MCMC方法就是这一节我们所要讨论的内容。 1 贝叶斯推断的计算问题在上节我们看到，贝叶斯统计学是利用后验分布对θ进行推断。这种推断的计算很多情况下要用积分计算来完成。比如，我们要计算θ的函数g θ 的期望：E g θ∣x ∫g θ fθ∣x θ∣x dθ其中函数f表示后验分布。当g θ θ时，得到的就是关于θ的点估计。但是对很多贝叶斯推断问题来说，有时候后验分布过于复杂，使得积分没有显示结果，数值方法也很难应用；有时候需要计算多重积分（比如后验分布是多元分布时）。这些都会带来计算上的很大困难。这也是在很长的时期内，贝叶斯统计得不到快速发展的一个原因。1990年代MCMC（Markov Chain Monte Carlo ，马尔科夫链蒙特卡洛）计算方法引入到贝叶斯统计学之后，一举解决了这个计算的难题。可以说，近年来贝叶斯统计的蓬勃发展，特别是在各个学科的广泛应用和MCMC方法的使用有着极其密切的关系。（2）蒙特卡洛方法（Monte Carlo）蒙特卡洛方法是一种随机模拟方法，随机模拟的思想由来已久（参见下面的蒲丰投针的例子），但是由于难于取得随机数，随机模拟的方法一直发展缓慢。而蒙特卡洛方法的出现得益于现代电子计算机的诞生，在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θ怎么样通过抽样的方法来解决这个问题呢？ 在θ的分布中抽取独立样本θ1,θ2,...,θn，n足够大 用样本均值θˉ 1n∑ni 1θi作为E θ 当n趋近无穷时，利用大数定律，样本均值会收敛到E θ （这称为仿真一致性） 从上面这个简单例子，蒙特卡洛方法的基本思