蒙特卡罗算法

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,Monte Carlo Simulation Methods(,蒙特卡罗模拟方法,),主要内容,:,一,. M-C,方法概述,.,二,.,随机数的生成,.,三,.,模拟训练,.,四,.,实验题目,.,成信院数学与信息科学系,1,蒙特卡洛,(Monte Carlo),方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯,诺伊曼用驰名世界的赌城,摩纳哥的,Monte Carlo,来命名这种方法，为它蒙上了一层神秘色彩。,一,.M-C,方法概述,2,基本思想很早以前就被人们所发现和利用。,17,世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。,19,世纪人们用投针试验的方法来决定,。高速计算机的出现，使得用数学方法在计算机上大量模拟这样的试验成为可能。,3,从Buffon(蒲丰)投针问题谈起,4,5,试验者,时间,(,年,),针长,投针次数,相交次数,的估计值,Wolf,1850,0.80,5000,2532,3.15956,Smith,1855,0.60,3204,1218,3.15665,Fox,1884,0.75,1030,489,3.15951,Lazzarini,1925,0.83,3408,1808,3.14159292,6,数值积分问题,选取,(0, 1),中随机数序列,x,1,，,x,2,，,x,3,，,x,n,。则,误差约 ，它并不能和一些高级的数值积分算法比拟，,但对多维情况，,MC,方法却很有吸引力。,7,我们可产生一系列随机数,可简单取,3,个随机数构成一个随机点，即,相应地，,一般地，,8,Monte Carlo数值积分的优点,与一般的数值积分方法比较，,Monte Carlo,方法,具有以下优点：,9,M-C的基本思路,1.,针对实际问题建立一个简单且便于实现的概率统计模,型，使所求的量（或解）恰好是该模型某个指标的概率分布或者数字特征。,2.,对模型中的随机变量建立抽样方法，在计算机上进行,模拟测试，抽取足够多的随机数，对有关事件进行统计,3.,对模拟试验结果加以分析，给出所求解的估计及其精,度,(,方差,),的估计,4.,必要时，还应改进模型以降低估计方差和减少试验费,用，提高模拟计算的效率,10,回顾几种连续型分布,1.,均匀分布,U(a,b),其概率密度函数为,有,11,均匀性特点,:,均匀分布随机变量,X,落在,(a,b),内,任意子区间的概率只与子区间的长度有关,而与,子区间的位置无关,.,可假设有这种特性的随机变量服从均匀分布,.,均匀分布随机变量,X,的取值具有,”,均匀性”,12,2.,正态分布,正态分布随机变量,X,的概率密度函数是,正态分布由两个参数 和 唯一确定,.,其中 是,X,的均值,(,数学期望,): =E(X),它确定了概率曲线的,中心位置,而 是,X,的标准差,: ,它确定了,概率曲线的,”,宽窄”程度,.,13,在许多实际问题中,有一类随机变量可以表,示成为许多相互独立的随机变量之和,而其中每,个随机变量对总和只起微小的影响,这类随机变,量往往服从或近似服从正态分布,.,在实际应用中,如果我们分析到一个随机变量受到较多独立的,微小因素的叠加影响,就可以用正态分布来模拟,这个变量,.,如,:,工厂产品的测量尺寸,农作物的收,获量,某地区成年人的身高,体重等可看成服从正,态分布的随机变量,.,14,3.,指数分布,指数分布随机变量,X,的概率密度为,指数分布常用来描述寿命问题,.,15,二.随机数的生成,1.,蒙特卡罗模拟的关键是生成优良的随机数。,2.,在计算机实现中,我们是通过确定性的算法生成,随机数,所以这样生成的序列在本质上不是随机,的,只是很好的模仿了随机数的性质,(,如可以通过,统计检验,),。我们通常称之为伪随机数,(pseudo-random numbers),。,3.,在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布,U(0,1),的随机数。,16,U(0,1)随机数的生成,乘同余法：,称为种子,a,是乘因子,m,是模数,17,一个简单的例子,18,一个简单的例子(续),上面的例子中，第一个随机数生成器的周期长度是,10,，而后两个生成器的周期长度只有它的一半。我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布。,19,线性同余生成器,(,混合同余法,)(Linear Congruential Generator ),c,是非负整数,.,通过适当选取参数,c,可以改善,随机数的统计性质,(,独立性,均匀性,).,20,常用的线性同余生成器,21,复杂一些的生成器,Multiple recursive generator,22,算法实现,许多程序语言中都自带生成随机数的方法，如,c,中的,random(),函数，,Matlab,中的,rand(),函数等。,但这些生成器生成的随机数效果很不一样，比如,c,中的函数生成的随机数性质就比较差，如果用,c,，最好自己再编一个程序。,Matlab,中的,rand(),函数，经过了很多优化。可以产生性质很好的随,机数，可以直接利用。,23,从,U(0,1),到其它概率分布的随机数,1.,离散型随机数的模拟,2.,连续型随机数的模拟,3.,正态随机数的模拟,24,1.,离散型随机数的模拟,设随机变量 的分布律为,令,将 作为区间,(0,1),的分点,.,若随机变量,有,25,令,则有据此,可得产生 的随机数的具体过程为,:,每产生一个,(0,1),区间上均匀分布随机数,若 则令 取值,.,26,例,1:,离散型随机变量,X,有如下分布律,:,X 0 1 2,P(x) 0.3 0.3 0.4,设 是,(0,1),上均匀分布的随机数,令,则 是具有,X,分布律的随机数,.,27,2.连续型随机数的模拟,a.,逆变换方法,(,常用,),(Inverse Transform Method),b.,舍取方法,(Acceptance-Rejection Method),定理,:,设随机变量,Y,的分布函数,F(y),是连续函数,而,U,是在,(0,1),上均匀分布的随机变量,令,则,Y,与,X,有相同的分布,.,28,29,30,31,舍取方法,舍取方法,(Acceptance-Rejection ),方法最早由,Von Neumann,提出，现在已经广泛应用于各,种随机数的生成。,基本思路：,通过一个容易生成的概率分布,g,和一个取舍,准则生成另一个与,g,相近的概率分布,f,。,32,具体步骤：,33,下面我们验证由上述步骤生成的随机数,Y,确实,具有密度函数,f(x),34,所以为了提高舍取法的效率，我们应该使,c,的取值尽,可能的小，也就是使,f,和,g,的分布更为相近。,35,3.标准正态分布随机数的生成,正态分布是概率统计中最重要的分布，在此,我们着重讨论如何生成标准正态分布随机数。,引理：,36,Box-Muller 算法,37,利用中心极限定理,设 是,n,个相互独立的在,(0,1),上均匀分布,的随机变量,由中心极限定理知,渐近服从正态分布,N(0,1).,一般取,n=10,即可,若取,n=12,则上式简化为 再由公式 即可,得到正态分布 的随机数,.,38,Matlab程序,Function r=rnd-u(a,b),%,产生在,a,b,间均匀分布的随机数,r=a+(b-a)*rand;,return,39,Matlab程序,Function r=rnd-beta(lamda),%,模拟指数分布,%lamda,表示指数分布的参数,r=-log(rand)/lamda;,return,40,Matlab程序,Function y= rnd(mean, segema),%,模拟均值为,mean,方差为,segma,的正态分布,r=rand(1,12);,x=sum(r)-6;,y=segma*x+mean;,return,41,三. 模拟训练,例,1.,模拟求近似圆周率,在边长为,1,的正方形内有一半径为,0.5,的内切圆,.,现在模拟产生在正方形内均匀分布的点,n,个,.,如,果有,m,个在圆内,则圆面积与正方形的面积比可,近似为,m/n.,即,/4,m/n,4m/n,42,n=10000,m=0;,For i=1:n,if rand2+rand2 =0,则可以通过模拟,估算,.,构造一个矩形包含曲边梯形,dmax f(x).,产生,n(,足够大,),个在矩形区域内的点,如果落在由函数,f(x),构成的曲边梯形内的点为,m,个,则所求定积分为,.,44,n=106;,a=0;,b=1;,d=1;,m=0;,for i=1:n,x=a+rand*(b-a);,y=d*rand;,if y=x2,m=m+1;,end,end,s=m/n*(b-a)*d,45,n=106;,x=rand(1,n);,y=x.2;,s=sum(y)/n,采用前面讲的方法,:,46,例,3.,渡口模型,问题描述,:,一个渡口的渡船营运者拥有一只甲板长,32,米,可以并排,停放两列车辆的渡船,.,他在考虑怎样在甲板上安排过河,车辆的位置,才能安全地运过最多数量的车辆,.,分析,:,怎样安排过河车辆,关心一次可以运多少辆各类,车,.,准备工作,:,观察数日,发现每次情况不尽相同,得到下列,数据和情况,:,(1),车辆随机到达,形成一个等待上船的车列,;,(2),来到车辆,轿车约占,40%,卡车约占,55%,摩托车约占,5%;,(3),轿车车身长,3.55.5,米,卡车车身长为,810,米,.,47,问题分析,:,这是一个机理较复杂的随机问题,是遵,”,先到先,服务”的随机排队问题,.,解决方法,:,采用模拟模型方法,.,因此需考虑以下问题,: (1),应该怎样安排摩托车,?,(2),下一辆到达的车是什么类型,?,(3),怎样描述一辆车的车身长度,?,(4),如何安排到达车辆加入甲板上两列车队中的哪一列中去,?,本实验主要模拟装载车辆的情况,暂时不考虑渡船的安全,.,48,模型建立,:,设到达的卡车,轿车长度分别为随机变量, .,结合实际,这里不妨设卡车,轿车的车身长度,均服从正态分布,.,由于卡车车身长,810,米,所以卡车车长 的均,值为,(8+10)/2=9,米,由概率知识中的,”3,”,原则,其标准差为,(9-8)/3=1/3,所以得到,.,同理可得,.,49,模拟程序设计,:,由以上的分析,程序设计时应划分的主要模块,(,函数,),如下,:,(1),确定下一辆到达车辆的类型,:,(2),根据车的类型确定到达车辆的长度,;,(3),根据一定的停放规则,确定放在哪一列,.,50,function r =makeid,%,模拟下一辆到达车的类型,t=rand;,if t=0.55,r=1;%,到达卡车,elseif tL(2),if L(2)+newlen32,pos=2;,else,full=1;,end,else,if L(1)+newlen32,pos=1;,else,full=1;,end,end,53,function sim_dukou,%,渡口模型的模拟,n=input(,输入模拟次数,:);,if isempty(n)|(n500),n=500;,end,N=zeros(1,3);,for i=1:n,isfull=0;,L=0,0;,while isfull,id=makeid;,N(id)=N(id)+1;,newlen=getlength(id);,isfull,pos=getiffull(L,newlen);,if isfull,L(pos)=L(pos)+newlen;,end,end,end,disp(,平均每次渡船上的车数,),mean_n=N/n,54,四.实验题目,赶上火车的概率,问题描述,:,如图,一列火车从,A,站开往,B,站,某人每,天赶往,B,站上这趟火车,.,他已了解到,:,(1),火车从,A,站到,B,站的运行时间是均值为,30,分钟,标准差为,2,分钟的随机变量,;,(2),火车在下午大约,1,点离开,A,站,离开时刻的频率分布如下,:,55,出发时刻,午后,1:00,午后,1:05,午后,1:10,频率,0.7,0.2,0.1,此人到达,B,站的时刻频率分布为,:,时刻,午后,1:28,午后,1:30,午后,1:32,午后,1:34,频率,0.3,0.4,0.2,0.1,问他能赶上火车的概率是多少,?,56,变量说明,:,火车从,A,站出发的时刻,;,:,火车从,A,站到,B,站的运行时间,单位,:,分钟,;,:,他到达,B,站的时刻,问题分析与假设,:,此题包含多个随机因数,.,这里假设, ,都是随机变量,其中 服从正态分布,.,57,模型建立,很显然,他能及时赶上火车的条件是,: + .,为了简化计算,将下午,1,点记为,0,时刻,.,和 的分布,律如下,:,/min,0,5,10,P(t),0.7,0.2,0.1,/min,28,30,32,34,P(t),0.3,0.4,0.2,0.1,58,如果,r,为在,(0,1),均匀分布的随机数,为了模拟随,机变量, ,可以通过如下方法,:,其中,和 分别用来模拟随机变量 和,.,59,function pro_huoche,%,模拟赶上火车的概率,n=input(,输入模拟次数,:),m=0;,for i=1:n,r1=rand;,if r1=0.7,t1=0;,elseif r1=0.9,t1=5;,else,t1=10;,end,r3=rand;,if r3=0.3,t3=28;,elseif r3=0.7,t3=30;,elseif r3=0.9,t3=32;,else,t3=34;,end,t2=30+randn*2;,if t3(t1+t2),m=m+1;,end,end,p=m/n,end,60,