随机数产生,MCMC方法,EM算法,bootstrap方法

1、统计计算文献报告摘 要:本文主要围绕随机数的产生算法、MCMC算法、EM算法、Bootstrap算法四种统计计算方法，介绍各方法的算法原理；并针对每个算法给出了实例。关键词：随机数，MCMC方法，EM算法，bootstrap方法1 随机数的产生统计模拟的最基本问题是随机数生成，是按统计模型生成数据的基础。随机数生成可分成两类:0,1区间上均匀随机数生成和非均匀随机数生成，是随机数生成问题的两个基本研究领域1。1.1 0,1区间上均匀随机数生成0,1上均匀随机数生成是随机数生成之基础，非均匀随机数是由0,1上均匀随机数经相应运算而产生的。0,1上均匀随机数的质量决定了非均匀随机数的质量

2、。如何快速地生成0,1上高质量均匀随机数是普遍关注的问题。生成方法分为两类:物理方法和数学方法2。1.1.1 物理方法基本原理为：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。把具有随机性质的物理过程变换为随机数，使用物理随机随机数发生器，在计算机上可以得到真正的随机数，其随机性和均匀性都是很好的、而且是取之不尽用之不竭的。但是用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，物理随机数发生器的稳定性经常需进行检查和维修，费用昂贵。因此，大大降低了这种方法的使用价值3。1.1.2 数学方法在计算机上

3、用数学方法产生均匀随机数是指按照一定的计算方法而产生的数列，它们具有类似于均匀随机变量的独立抽样序列的性质，这些数既然是依照确定算法产生的，便不可能是真正的随机数，因此常把用数学方法产生的随机数称为伪随机数。用数学方法在计算机上产生随机数是目前使用最广、发展很快的一类方法。他的特地啊是占用内存少、速度快又便于复算。1.1.3 同余法目前应用最广泛的随机数产生的方法是同余法。它是由Lehmer在1951年提出的。此方法利用数论中的同于运算来产生随机数，故称为同余发生器。它包括混合同余发生器和乘同于发生器。它是目前使用最普遍、发展迅速的产生随机数的数学方法。基本方法如下： （1.1）其中为模数，为

4、乘子(乘数)，为增量（加数），且，均为非负整数。就是个0,1上的伪随机数。按（1.1）给出的为随机数生成方法称为线性同余法。当时称为混合同余法，时称为乘同余法，当且模是素数时，称为素数模乘同余法。线性同余法产生的随机数序列具有周期性，其最大周期为。著名的Coveyou与Macpherson混合同于发生器为：Kobayashi混合同余器是1.2 非均匀随机数生成非均匀随机数的生成是另一类随机数生成方法，是基于0,1上的均匀随机数产生给定概率密度函数的算法。如正态分布随机数发生器，指数分布随机数发生器等，可以利用分布的特殊性质，针对性极强。在一些实际问题或统计模拟时需要上万，几十万，甚至上百万个非

5、均匀随机数。这时，算法的速度就显得十分重要。目前产生非聚云随机数的一般方法有：直接抽样法（反函数法）、变换抽样法、值序抽样法和舍选抽样法。下面针对反函数变换法和舍选法进行介绍。1.2.1 反函数变换法设随机变量的分布函数是，且是连续的。那么，是随机变量，其分布是(0,1)上均匀分布。于是，有产生具有给定分布函数的非均匀随机数的算法: 这里是(0,1)上的均匀分布，是的反函数。是一个通用算法，对任何连续分布函数的反函数变换法，简称反函数变换法。该方法的缺点是:并不是任何分布函数都有快速算法。对常见分布，如正态分布，F分布等分布函数仅可以数值计算，没有解析表达式。1.2.2 舍选法舍选法是非常有用

6、的随机数生成方法。其应用最为广泛，研究成果也最多。设是随机变量的密度函数，如何生成有密度函数的随机数是随机数生成中最重要的和最基本的问题。最普遍应用的方法就是舍选法。方法如下:1.存在函数满足以下条件: （1.2）也是密度函数，而且密度为的随机变量容易生成。 2.生成随机变量和随机数; 3.若，则，否则转2; 4.输出, 的密度是 满足(1.2)式的函数叫做函数的优函数。不附加任何条件的优函数是容易找到的。1.3 随机数生成实例运用反函数法，用生成的100个随机数。因为指数分布的分布函数为：当时，反函数。利用R语言进行随机数生成。反函数法生成100个的指数分布：set.seed(100)FIn

7、verse<-function(y,lambda=1) -log(1-y)/lambdarandomSequence<-FInverse(runif(100,0,1)plot(randomSequence,type="l")为检验利用反函数法生成的随机数是否为指数分布，进行K-S检验，检验结果如下：从输出结果可以看出，故不能拒绝反函数法生成的序列来自于。2 MCMC算法马尔可夫链蒙特卡罗(MCMC)方法就是其中一种应用广泛的模拟方法。MCMC方法最早由蒙特卡罗(MC)方法发展而来，蒙特卡罗方法是一种根据大数定律重复模拟得到近似概率的过程。MCMC的发展最早是由统

8、计物理学家Metropolis等人将马氏链引入MC方法中，解开了普通Monte Carlo方法不改变抽样分布的瓶颈，使之能解决许多复杂的物理系统。在这之后随着计算机技术的蓬勃发展，极大的推动了MCMC方法的广泛应用，在统计学领域中，极大似然估计、非参数估计、贝叶斯统计学等方向上的很多问题得到了有效解决。2.1 MCMC的基本思想MCMC方法的基本思想是通过构建一个平稳分布为的Markov能得到服从的样本, 基于此样本对的分布特征进行各种统计推断。2.1.1 Markov链若随机变量序列满足则称此随机变量序列为Markov链。即在已知“现在”的条件下，“将来”与“过去”是相互独立的，“将来”仅与

9、“现在”有关，与“过去”无关。2.1.2 Markov链的转移核及平稳分布设为Markov链。其一步转移概率函数为：（离散）或对于任一可测集,有 (连续)称为该Markov链的转移核，通常假定与时间t无关(即时齐Markov链)。如果分布满足：则称为转移核的平稳分布。2.1.3 MCMC方法步骤MCMC方法可概括为如下三个基本步骤：(1)在样本空间上选择“合适”的Markov链，使为其相应的平稳分布；(2)由中某一点出发，利用(1)中的Markov链产生序列；(3)对给定的m和充分大的n，序列可认为是来自的遍历样本，可基此对或其数字特征进行估计。任一函数的期望估计为：MCMC方法的重点是构造转

10、移核, 使给定的概率分布为相应的Markov链的平稳分布。2.2 MCMC算法实例用MCMC算法的办法生成。设初始点为0.1，利用的随机数进行模拟变换。N <- 3000 #抽样个数 x <- c() x0 <- 0.1 # 初始值 x1 <- x0 k <- 0 #k表示拒绝转移的次数 u <- runif(N) #抽取均匀分布随机数 for (i in 2:N) y <- xi-1+runif(1,-0.2,0.2) if(0<y&y<1) if (ui <=(dbeta(y,2,7)/dbeta(xi-1,2,7) xi

11、 <- y else xi<- xi-1 k <-k + 1 else xi<-x0 得到的x是高度自相关的，要得到的分布，需要对x进行采样，采样的步长通过自相关系数可确定，从下图可看出自相关系数在延迟30阶后均在二倍标准差内，故确定步长为30。acf(x)进一步判断x是否达到了一个平稳分布，做出x的样本路径。plot(x,type="l")可看出，样本x已是平稳的，没有漂移的情况发生。可以认为这个马氏链的构造是成功的。开始采样，从1000开始，每隔30步取一个数，构成新的数列y。我们对y与x分别作K-S检验，有y<-xseq(1000,N,b

12、y=30)ks.test(y,"pbeta",2,7)从输出结果可以看出，故不能拒绝序列来自于的原假设。3 EM 算法在统计领域里，主要有两大类计算问题，一类是极大似然估计的计算，另一类是Bayes计算。从计算方法上讲，二者是可以合并讨论的，因极大似然估计的计算类似于Bayes的后验众数的计算。Bayes计算方法已经有很多，大体上可以分为两大类：一类是直接应用于后验分布以得到后验均值或后验众数的估计，以及这种估计的渐进方差或其近似。另一类算法总称为数据添加算法。这是近年发展很快且应用很广的一种算法，它不是直接对复杂的后验分布进行极大化或进行模拟，而是在观测数据的基础上加上一

13、些“潜在数据”，从而简化计算并完成一系列简单的极大化或模拟，该“潜在数据”可以使“缺损数据”或未知参数。3.1 EM算法基本原理设我们能观测的数据,关于的后验分布很复杂，难以直接进行各种统计计算，加入我们能假定一些没有能观测到的潜在数据为已知（譬如，为某变量的结尾观测值，而为该变量的真值），则可能得到一个关于的简单的添加后验分布，利用的简单性我们可进行各种统计计算，如极大化，抽样等，然而回过头来，我们又可以对的假定做检查和改进，如此进行，我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。3.2 EM算法基本思想与基本步骤EM算法是一种迭代方法，最初是由Dempster等提出，并

14、主要用来求后验分布的众数（即极大似然估计），它的每一次迭代由两步组成：E步（求期望）和M步（极大化）。一般地，以表示的基于观测数据的后验分布密度函数，成为观测后验分布，表示在给定和观测数据下潜在数据的条件分布密度函数。目的在于计算观测后验分布的众数，于是，EM算法如下进行。记为第次迭代开始时后验众数的估计值，则第次迭代的两步为：E步：将或关于的条件分布求期望，从而把积掉，即M步：将极大化，即找一个点，使如此形成了一次迭代。将上述E步和M步进行迭代直至或充分小时停止。3.3 EM算法实例考虑模型：在第组中的第个随机变量为，观测值为，且假设有如下观测值：x1=c(62,60,63,59)x2=c(

15、63,67,71,64,65,66)x3=c(68,66,71,67,68,68)x4=c(56,62,60,61,63,64,63,59)x=list(x1,x2,x3,x4)EM=function(u0,v0,w0,n,x) xx=unlist(x) ux=rep(0,n) vx=rep(0,n) wx=rep(0,n) ux1=u0 vx1=v0 wx1=w0 m=rep(0,length(x) x0=rep(0,length(x) u=matrix(0,n,length(x) v=matrix(0,n,length(x) j=1 while(j<n) for (i in 1:le

16、ngth(x) mi=length(xi) x0i=sum(xi)/length(xi) uj,i=(mi*x0i/vxj2+uxj/wxj2)/(mi/vxj2+1/wxj2) vj,i=(mi/vxj2+1/wxj2)(-1) uu=rep(uj,m) vv=rep(vj,m) u1=uj, v1=vj, uu vv j=j+1 uxj=sum(u1)/length(u1) vxj=(sum(uu-xx)2+vv)/length(xx)(1/2) wxj=(sum(u1-uxj)2+v1)/(length(x)-1)(1/2) return(list(ux=ux,vx=vx,wx=wx)E

17、M(64,2.29,3.56,8,x)EM(60,2,2,8,x)em=EM(23,5,1,90,x)empar(mfcol=c(1,2)plot(em1,type = "l",ylim=c(0,70),xlab = "迭代次数",ylab = "参数值",main = "EM算法在初始参数为23,5,1的迭代结果",col=6)lines(em2,col=4)lines(em3,col=5)参数 初值 1 2 3 4 5 6 7 8 64.00000 64.01042 64.01199 64.01225 64.0

18、1230 64.01230 64.01231 64.01231 64.01231 2.290000 2.350002 2.359641 2.361218 2.361477 2.361520 2.361528 2.361529 2.361529 3.560000 3.486979 3.473919 3.471498 3.471058 3.470979 3.470965 3.470963 3.470962 60.00000 63.42222 63.96079 64.00863 64.01208 64.01230 64.01231 64.01231 64.01231 2.000000 2.39892

19、7 2.369358 2.362940 2.361778 2.361572 2.361536 2.361530 2.361529 2.000000 3.185404 3.424427 3.463902 3.469882 3.470794 3.470936 3.470958 3.4709624 Bootstrap抽样Bootstrap方法是一种用于模拟频率论的推理方法.最初是由B.Efron作为推导任意估计值的标准误的一种方法，从B.Efron首次系统提出到现在的30多年间，学者已经将Bootstrap方法迅速的运用于各个领域，如医学、生物统计、金融等各领域的实际问题研究中。在当今这个计算机发展

20、相当迅速的时代，Bootstrap方法的应用领域也在不断的拓展，其发展前景是十分广阔的。4.1 Bootstrap方法简介Bootstrap方法也称“自助法”，是由B. Efron提出的一种高度依赖计算机的重抽样方法。它可以通过从观测样本中多次“有放回抽样”技术，近似模拟样本的真实统计分布函数和参数。根据统计推断理论，假设一个未知统计分布的有限样本容量，从观测数据中抽取B个独立的Bootstrap伪数据集，通过增大B可以使得模拟误差任意小。因此这是一种比传统方法更精确和有效的统计推断方法，可以用于模拟受到多种因素影响的样本。4.2 Bootstrap方法的基本思想Bootstrap方法的核心是通过再抽样来构造自助样本。Bootstrap方法经常用于下列情况: 1.标准假设无效(如样本容量n很小，样本不服从正态分布等); 2.需要解决复杂问题，且没有理论可依;除此之外，其他任何地方也可以用自助法。Bootstrap方法的基本思想是重抽样和Bootstrap分布。Bo