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

统计计算文献报告摘要:本文主要围绕随机数的产生算法、MCMC算法、EM算法、Bootstrap算法四种统计计算方法，介绍各方法的算法原理；并针对每个算法给出了实例。关键词：随机数，MCMC方法，EM算法，bootstrap方法1．随机数的产生统计模拟的最基本问题是随机数生成，是按统计模型生成数据的基础。随机数生成可分成两类:［0,1］区间上均匀随机数生成和非均匀随机数生成，是随机数生成问题的两个基本研究领域［1］。1.1［0,1］区间上均匀随机数生成［0,1］上均匀随机数生成是随机数生成之基础，非均匀随机数是由［0,1］上均匀随机数经相应运算而产生的。［0,1］上均匀随机数的质量决定了非均匀随机数的质量。如何快速地生成［0,1］上高质量均匀随机数是普遍关注的问题。生成方法分为两类:物理方法和数学方法［2］。1.1.1物理方法基本原理为：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。把具有随机性质的物理过程变换为随机数，使用物理随机随机数发生器，在计算机上可以得到真正的随机数，其随机性和均匀性都是很好的、而且是取之不尽用之不竭的。但是用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，物理随机数发生器的稳定性经常需进行检查和维修，费用昂贵。因此，大大降低了这种方法的使用价值［3］。1.1.2数学方法在计算机上用数学方法产生均匀随机数是指按照一定的计算方法而产生的数列，它们具有类似于均匀随机变量的独立抽样序列的性质，这些数既然是依照确定算法产生的，便不可能是真正的随机数，因此常把用数学方法产生的随机数称为伪随机数。用数学方法在计算机上产生随机数是目前使用最广、发展很快的一类方法。他的特地啊是占用内存少、速度快又便于复算。

1.1.3同余法目前应用最广泛的随机数产生的方法是同余法。它是由Lehmer在1951年提出的。此方法利用数论中的同于运算来产生随机数，故称为同余发生器。它包括混合同余发生器和乘同于发生器。它是目前使用最普遍、发展迅速的产生随机数的数学方法。基本方法如下：x=ax+x=ax+b1xn n—1(modm)=(ax+c) (modm)1.1)u=u=x/mnnn=1,2,…其中m为模数，a为乘子（乘数），c为增量（加数），且xn，m，a，c均为非负整数。u,u，…,u就是k个［0,1］上的伪随机数。按（1.1）给出的为随机数生1 2 k成方法称为线性同余法。当b>0时称为混合同余法，b=0时称为乘同余法，当b=0且模是素数时，称为素数模乘同余法。线性同余法产生的随机数序列具有周期性，其最大周期为m。著名的Coveyou与Macpherson混合同于发生器为：x=515x+1mod（235）TOC\o"1-5"\h\z<n n—1u=x/235 n=1,2,••-n nKobayashi混合同余器是”x=314159269x+453806245mod（231）<n n—1u=x/231 n=1,2,••-n n1.2非均匀随机数生成非均匀随机数的生成是另一类随机数生成方法，是基于［0,1］上的均匀随机数产生给定概率密度函数的算法。如正态分布随机数发生器，指数分布随机数发生器等，可以利用分布的特殊性质，针对性极强。在一些实际问题或统计模拟时需要上万，几十万，甚至上百万个非均匀随机数。这时，算法的速度就显得十分重要。目前产生非聚云随机数的一般方法有：直接抽样法（反函数法）、变换抽样法、值序抽样法和舍选抽样法。下面针对反函数变换法和舍选法进行介绍。1.2.1反函数变换法设随机变量X的分布函数是F(。)，且F是连续的。那么，U=F(X)是随机变量，其分布是(0,1)上均匀分布。于是，有产生具有给定分布函数F(。)的非均匀随机数X的算法:X=F-i(U),U〜U(0,1)这里U(0,1)是(0,1)上的均匀分布，F-1(。)是F的反函数。是一个通用算法，对任何连续分布函数的反函数变换法，简称反函数变换法。该方法的缺点是:并不是任何分布函数都有快速算法。对常见分布，如正态分布，F分布等分布函数仅可以数值计算，没有解析表达式。1.2.2舍选法舍选法是非常有用的随机数生成方法。其应用最为广泛，研究成果也最多。设f(。)是随机变量X的密度函数，如何生成有密度函数f(。)的随机数是随机数生成中最重要的和最基本的问题。最普遍应用的方法就是舍选法。方法如下:1.存在函数M(。)满足以下条件：f(x)<M(x),JM(x)dx=C<g,m(x)=M"X) (1.2)C—gm(。)也是密度函数，而且密度为m(。)的随机变量Y容易生成。生成随机变量y〜m(。)和随机数V〜U(0,1)；若M(Y)V<f(Y)，则X=Y，否则转2;输出X,X的密度是f(。)满足(1.2)式的函数M(。)叫做函数f的优函数。不附加任何条件的优函数是容易找到的。

1.3随机数生成实例运用反函数法，用U(0,1)生成exp(l)的100个随机数。因为指数分布的分布函数为：/、 [1—e—儿x>0F(x)斗[0 x<0当x>0时，反函数F-1(x)=—log(1—y)。九利用R语言进行随机数生成。反函数法生成100个九=1的指数分布:set.seed(100)FInverse<-function(y,lambda=1){-log(1-y)/lambda}randomSequence<-FInverse(runif(100,0,1))plot(randomSequence,type="l")Q 2D 鈕 "IMinam为检验利用反函数法生成的随机数是否为指数分布，进行K-S检验，检验结果如下：> k：3.匸mm匸[imnilciniSeque二eefrir?expr,fL)OnE-BampleKoIrr.ogoiov-Sm_irnovtest-data: randoruSeqaenceD=0.10026fp-valae=0.26^2altemat-^VEhypothESis:two-sided从输出结果可以看出p=0.2672》0.05，故不能拒绝反函数法生成的序列来自于exp(1)。2．MCMC算法马尔可夫链蒙特卡罗(MCMC)方法就是其中一种应用广泛的模拟方法。MCMC方法最早由蒙特卡罗(MC)方法发展而来，蒙特卡罗方法是一种根据大数定律重复模拟得到近似概率的过程。MCMC的发展最早是由统计物理学家Metropolis等人将马氏链引入MC方法中，解开了普通MonteCarlo方法不改变抽样分布的瓶颈，使之能解决许多复杂的物理系统。在这之后随着计算机技术的蓬勃发展，极大的推动了MCMC方法的广泛应用，在统计学领域中，极大似然估计、非参数估计、贝叶斯统计学等方向上的很多问题得到了有效解决。MCMC的基本思想MCMC方法的基本思想是通过构建一个平稳分布为兀(x)的Markov能得到服从兀(x)的样本，基于此样本对兀(x)的分布特征进行各种统计推断。Markov链若随机变量序列｛X(0),X⑴,…,X(n),・・・｝满足P(X(n+1)<xIX(0)=x(0),X(1)=X⑴…，X(n)=x(n))=P(X(n+1)<x|X(n)=x(n))贝席此随机变量序列为Markov链。即在已知“现在”的条件下，“将来”与“过去”是相互独立的，“将来”仅与“现在”有关，与“过去”无关。2.1.2Markov链的转移核及平稳分布设｛X(0),X(1),X⑵…｝为Markov链。其一步转移概率函数为：p(x,xr)=P(xTx')=P(X(t+1)=x'IX(t)=x)(离散)或对于任一可测集B ,有JVP(x—B)=Jp(x,x')dxx(连续)BP(-,•)称为该Markov链的转移核，通常假定P(-,•)与时间t无关(即时齐Markov链)。如果分布兀(x)满足：Jp(x,x')兀(x)dx=兀(x')则称兀(x)为转移核P(・,•)的平稳分布。2.1.3MCMC方法步骤MCMC方法可概括为如下三个基本步骤：x(1)在样本空间 上选择“合适”的Markov链，使兀(x)为其相应的平稳分布；x⑵由 中某一点X(。)=x(o)出发，利用⑴中的Markov链产生序列X(1),X(2)，…,X(n)；(3)对给定的m和充分大的n序列X(m+i),X(m+2)，…,XG)可认为是来自兀(x)的遍历样本，可基此对兀(x)或其数字特征进行估计。任一函数f(x)的期望估计为：Ef=—^£f(X(J兀 n一mt=m+1MCMC方法的重点是构造转移核，使给定的概率分布兀(x)为相应的Markov链的平稳分布。2.2MCMC算法实例用MCMC算法的办法生成Beta(2,7)。设初始点为0・1,利用U(-0.2,0・2)的随机数进行模拟变换。N<-3000#抽样个数x<-c()x0<-0.1#初始值x[1]<-x0k<-0#*表示拒绝转移的次数u<-runif(N)#抽取均匀分布随机数for(iin2:N){y<-x[i-1]+runif(1,-0.2,0.2)if(0<y&y<1){if(u[i]<=(dbeta(y,2,7)/dbeta(x[i-1],2,7)))x[i]<-yelse{x[i]<-x[i-1]k<-k+1}}elsex[i]<-x0}得到的x是高度自相关的，要得到Beta(2,7)的分布，需要对x进行采样，采样的步长通过自相关系数可确定，从下图可看出自相关系数在延迟30阶后均在二倍标准差内，故确定步长为30。acf(x)进一步判断x是否达到了一个平稳分布，做出x的样本路径。plot(x,type="l")构造是成功的。开始米样，从1000开始，每隔30步取一个数，构成新的数列y。我们对y与x分别作K-S检验，有y<-x[seq(1000,N,by=30)]ks.test(y,"pbeta",2,7)> ；3eq(L000fN^y=30)]、k3.te3t(yfOne-sai^leXolir.ogora-1?-Sn-imovtestddtd: yD=0.0S1761fp-vdlae=0.7616alternativehypothesis:fwo-sided从输出结果可以看出p0.7616》0.05，故不能拒绝序列来自于Beta(2,7)的原假设。3．EM算法在统计领域里，主要有两大类计算问题，一类是极大似然估计的计算，另一类是Bayes计算。从计算方法上讲，二者是可以合并讨论的，因极大似然估计的计算类似于Bayes的后验众数的计算。Bayes计算方法已经有很多，大体上可以分为两大类：一类是直接应用于后验分布以得到后验均值或后验众数的估计，以及这种估计的渐进方差或其近似。另一类算法总称为数据添加算法。这是近年发展很快且应用很广的一种算法，它不是直接对复杂的后验分布进行极大化或进行模拟，而是在观测数据的基础上加上一些“潜在数据”，从而简化计算并完成一系列简单的极大化或模拟，该“潜在数据”可以使“缺损数据”或未知参数。3.1EM算法基本原理设我们能观测的数据Y,9关于Y的后验分布p(9IY)很复杂，难以直接进行各种统计计算，加入我们能假定一些没有能观测到的潜在数据Z为已知(譬如，Y为某变量的结尾观测值，而Z为该变量的真值)，则可能得到一个关于9的简单的添加后验分布p(9IY,Z)，利用p(9IY,Z)的简单性我们可进行各种统计计算,如极大化，抽样等，然而回过头来，我们又可以对Z的假定做检查和改进，如此进行，我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样3.2EM算法基本思想与基本步骤EM算法是一种迭代方法，最初是由Dempster等提出，并主要用来求后验分布的众数(即极大似然估计)，它的每一次迭代由两步组成：E步(求期望)和M步(极大化)。一般地，以p(9IY)表示9的基于观测数据的后验分布密度函数，成为观测后验分布，p(Z19,Y)表示在给定9和观测数据Y下潜在数据Z的条件分布密度函数。目的在于计算观测后验分布p(9IY)的众数，于是，EM算法如下进行。记9⑴为第i+1次迭代开始时后验众数的估计值，则第i+1次迭代的两步为：E步：将p(9IY,Z)或logp(9IY,Z)关于Z的条件分布求期望，从而把Z积掉，即Q(919(i),Y)仝E[logp(9IY,Z)19⑴,Y]Z二门og[p(9IY,Z)]p(Z19(i),Y)dZM步：将Q(919(i),Y)极大化，即找一个点9(i+i)，使Q(9(i+1)19(i),Y)二maxQ(919⑴,Y)9如此形成了一次迭代9(i)t9(i+1)。将上述E步和M步进行迭代直至卩(i+1)-0(i)||或IQ(0(i+i)10(i),Y)-Q(0(i)10(i),Y)|充分小时停止。3.3EM算法实例考虑模型：在第i组中的第j个随机变量为X..，观测值为，且ij ij（XI卩Q）~N（卩Q2）,j=1,2,…,miji i假设有如下观测值： 表L观测数据 组别观测值如n3 Xi.16260635946126367716465666663686671676868668456626()6163646359861x1=c(62,60,63,59)x2=c(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)ux[1]=u0vx[1]=v0wx[1]=w0m=rep(0,length(x))x0=rep(0,length(x))u=matrix(0,n,length(x))v=matrix(0,n,length(x))j=1while(j<n){for(iin1:length(x)){m[i]=length(x[[i]])x0[i]=sum(x[[i]])/length(x[[i]])u[j,i]=(m[i]*x0[i]/vx[j]A2+ux[j]/wx[j]A2)/(m[i]/vx[j]A2+1/wx[j]A2)v[j,i]=(m[i]/vx[j]A2+1/wx[j]A2)A(-1)uu=rep(u[j,],m)vv=rep(v[j,],m)}u1=u[j,]v1=v[j,]uuvvj=j+1ux[j]=sum(u1)/length(u1)vx[j]=(sum((uu-xx)人2+vv)/length(xx))人(1/2)wx[j]=(sum((ul-ux[j])人2+vl)/(length(x)-l))人(1/2)}return(list(ux=ux,vx=vx,wx=wx))}EM(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(em[[1]],type="l",ylim=c(0,70),xlab="迭代次数”,ylab="参数值",main="EM算法在初始参数为23,5,1的迭代结果”,col=6)lines(em[[2]],col=4)lines(em[[3]],col=5)参数初值 1 2 3 4 5 6 7 8卩 64.00000 64.01042 64.01199 64.01225 64.01230 64.01230 64.01231 64.01231 64.012312.290000 2.350002 2.359641 2.361218 2.361477 2.361520 2.361528 2.361529 2.361529T 3.560000 3.486979 3.473919 3.471498 3.471058 3.470979 3.470965 3.470963 3.47096260.0000063.42222 63.96079 64.00863 64.01208 64.01230 64.01231 64.01231 64.012312.0000002.398927 2.369358 2.362940 2.361778 2.361572 2.361536 2.361530 2.3615292.0000003.185404 3.424427 3.463902 3.469882 3.470794 3.470936 3.470958 3.4709624．Bootstrap抽样Bootstrap方法是一种用于模拟频率论的推理方法.最初是由B.Efron作为推导任意估计值的标准误的一种方法，从B.Efron首次系统提出到现在的30多年间，学者已经将Bootstrap方法迅速的运用于各个领域，如医学、生物统计、金融等各领域的实际问题研究中。在当今这个计算机发展相当迅速的时代，Bootstrap方法的应用领域也在不断的拓展，其发展前景是十分广阔的。4.1Bootstrap方法简介Bootstrap方法也称“自助法”是由B.Efron提出的一种高度依赖计算机的重抽样方法。它可以通过从观测样本中多次“有放回抽样”技术，近似模拟样本的真实统计分布函数和参数。根据统计推断理论，假设一个未知统计分布的有限样本容量，从观测数据中抽取B个独立的Bootstrap伪数据集，通过增大B可以使得模拟误差任意小。因此这是一种比传统方法更精确和有效的统计推断方法，可以用于模拟受到多种因素影响的样本。样本总体样本估计量e估计量©怙计量0样本总体样本估计量e估计量©怙计量0图1重抽样方法的思想示意图〔左；传统的总体一样本*右：重抽样方法的样.本一子样本）4.2Bootstrap方法的基本思想Bootstrap方法的核心是通过再抽样来构造自助样本。Bootstrap方法经常用于下列情况:标准假设无效(如样本容量n很小，样本不服从正态分布等);需要解决复杂问题，且没有理论可依;除此之外，其他任何地方也可以用自助法。Bootstrap方法的基本思想是重抽样和Bootstrap分布。Bootstrap分布既包含了总体数据的信息，也包含了统计量抽样分布的信息。设总体X〜F(x;0),0e©为感兴趣的未知参数，F(x;0)=F未知。X=(X,X，…,X)为来自总体F的简单随机样本。在统计推断(如区间估计和1 2 n 0假设检验)中，一般都会涉及某个随机变量R(X;F)的分布0J(x;F)=P(R(X;F)<x)n0 F0 0经典统计方法中，求J(x;F)的方法主要包括：n0求J(x;F)的精确分布；n0求J(x;F)的极限