单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,Monte Carlo,模拟,内容提纲,1.,引言,2.Monte Carlo,模拟基本思想,3.,随机数生成函数,4.,应用实例举例,5.,排队论模拟,6.Monte Carlo,模拟求解规划问题,引言,(Introduction),Monte Carlo,方法:,蒙特卡罗方法,又称随机模拟方法,属于计算数学的一个分支,它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的亦称统计模拟方法,,statistical simulation method,利用随机数进行数值模拟的方法,Monte Carlo,名字的由来:,Monte Carlo,是摩纳哥(,monaco,),的首都,该城以赌博闻名,Nicholas Metropolis(1915-1999),Monte-Carlo,Monaco,Monte Carlo,方法的基本思想,蒙特卡罗方法,或称计算机随机模拟方法,是一种基于“随机数”的计算方法源于美国在第,二,次世界大战研制原子弹的“曼哈顿计划”,,,该计划的主持人之一数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法,为它蒙上了一层神秘色彩。
蒙特卡罗,方法的基本思想很早以前就被人们所发现和利用早在17世纪,人们就知道用事件发生的“频率”来决定事件的“概率”1,9,世纪人们用,蒲丰,投针的方法来,计算,圆周率,,上,世纪40年代电子计算机的出现,特别是近年来高速电子计算机的出现,使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能蒲丰投针实验,:,法国科学家蒲丰,(Buffon),在,1777,年提出的蒲丰,投针实验是早期几何概率一个非常著名的例子蒲丰,投针实验,的重要性并非是为了求得比其它方法更精确,的,值,而是它开创了使用随机数处理确定性数学问,题的先河,是用偶然性方法去解决确定性计算的前导由此,可以,领略到从“概率土壤”上开出的一朵瑰丽的鲜花,-,蒙特卡罗,方法,(,MC),蒲丰投针实验可归结为下面的数学问题:平面上画有距离为,a,的一些平行线,向平面上任意投一根长为,l(la),的针,假设针落在任意位置的可能性相同,试求针与平行线相交的概率,P(,从而求,),蒲丰投针实验,:,如右图所示,以,M,表示针落下,后的中点,以,x,表示,M,到最近一条平行,线的距离,以,表示针与此线的交角:,针落地的所有可能结果满足:,其样本空间视作矩形区域,面积是,:,针与平行线相交的条件:,它是样本空间,子集,A,,面积是:,syms,l,phi;int(l/2*sin(phi),phi,0,pi),%ans=,l,因此,针与平行线相交的概率为:,从而有:,蒲丰投针实验,的计算机模拟:,format long;%,设置,15,位显示精度,a=1;,l=0.6;,%两平行线间的宽度,和,针长,figure;axis(0,pi,0,a/2);%,初始化绘图板,set(gca,nextplot,a,dd,);,%,初始化绘图方式为叠加,counter=0;,n=,2010,;,%,初始化计数器和设定,投针次数,x=unifrnd(0,a/2,1,n);,p,hi,=unifrnd(0,pi,1,n);,%,样本空间,for i=1:n,if x(i)l*sin(,p,hi,(i)/2,%满足此条件表示针与线的相交,plot(phi(i),x(i),r.),;,counter=counter+1;,%,统计针与线相交的次数,frame(,counter,)=getframe;%,描点并取帧,end,end,fren=counter/n;,pihat=2*l/(a*fren),%,用,频率,近似计算,figure(2),movie(frame,1)%,播放帧动画,1,次,一些人进行了实验,其结果列于下表,:,实验者,年份,投计次数,的实验值,沃尔弗,(Wolf),1850,5000,3.1596,斯密思,(Smith),1855,3204,3.1553,福克斯,(Fox),1894,1120,3.1419,拉查里尼,(Lazzarini),1901,3408,3.1415929,蒙特卡罗投点法是,蒲丰投针实验,的推广:,在一个边长为,a,的正方形内随机投点,,该点落在此正方形的内切圆中的概率,应为该内切圆与正方形的面积比值,,即,n=10000;a=2;m=0;,for i=1:n,x=rand(1)*a;y=rand(1)*a;,if(x-a/2)2+(y-a/2)2=(a/2)2),m=m+1;,end,end,disp(,投点法近似计算的,为,:,num2str(4*m/n);,x,y,o,(a/2,a/2),基本思想,由上面的例子可以看出,当所求问题的解是某个事件的概率,或者是某个随机变量的数,学期望,或者是与之有关的量时,通过某种试验的方法,得出该事件发生的频率,再通过它得到问题的解。
这就是蒙特卡罗方法的基本思想蒙特卡罗方法的关键步骤在于随机数的产生,计算机产生的随机数都不是真正的随机数,(,由算法确定的缘故,),,如果伪随机数能够通过一系列统计检验,我们也可以将其当作真正的随机数使用rand(seed,0.1);,rand(1),每次运行程序产生的值是相同的,rand(state,sum(100*clock)*rand);,rand(1)%,每次重新启动,matlab,时,输出的随机数不一样,注意,:,产生一个参数为,的指数分布的随机数应输入,exprnd(1/,),产生mn阶,参数为,A1,A2,A3,的指定分布,name,的,随机数矩阵,random(name,A1,A2,A3,m,n),举例,:,产生,2,4阶的均值为,0,方差为,1,的正态分布的随机数矩阵,random(Normal,0,1,2,4),name,的取值可以是,(,详情参见,help random),:,norm or Normal/unif or Uniform,poiss or Poisson/beta or Beta,exp or Exponential/gam or Gamma,geo or Geometric/unid or Discrete Uniform,非常见分布的随机数的产生,逆变换方法,Acceptance-Rejection,方法,最早由,Von Neumann,提出,现在已经广泛应用于各种随机数的生成。
基本思路:,通过一个容易生成的概率分布,g,和一个取舍,准则生成另一个与,g,相近的概率分布,f,为要模拟服从给定分布的随机变量,用生成一个易于实现的不可约遍历链 作为随机样本,使其平稳分布为给定分布的方法,称为,马氏链蒙特卡罗方法,.,马氏链蒙特卡罗方法,1,连续型随机变量,(,以指数分布为例,),:,syms t x lambda;,Fx=int(,lambda,*exp(-,lambda,*t),t,0,x)%,分布函数,syms r;,Fxinv=finverse(Fx,x);%,求反函数,Fxinv=subs(Fxinv,x,r)%,替换反函数变量,x,为,r,Fxinv=inline(Fxinv),x=Fxinv(3,rand)%,产生参数,lambda=3,指数分布的随机数,%,指数分布随机数产生函数已经提供,exprnd(1/3,1,1),2,离散型随机变量,(,以离散分布为例,),:,x=2,4,6,8;px=0.1,0.4,0.3,0.2;%,以下为程序片段,Fx=0;for n=1:length(px),Fx=Fx,sum(px(1:n);end,r=rand;index=find(rMAXK,或,PMAXP,时停止迭代,框 图,初始化,:,给定,MAXK,MAXP;k=0,p=0,Q:,大整数,x,j,=a,j,+R(b,j,-a,j,)j=1,2,n,j=0,j=j+1,p=p+1,PMAXP?,Y,N,x,j,=a,j,+R(b,j,-a,j,),g,i,(X)0?,i=1,2n,Y,N,jMAXK?,Y,N,输出,X,Q,停止,Y,N,在,Matlab,软件包中编程,共需三个文件,:randlp.m,mylp.m,lpconst.m.,主程序为,randlp.m.,%mylp.m,function,z=mylp(x),%,目标函数,z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2);,%,转化为求最小值问题,%randlp.m,function,sol,r1,r2=randlp(a,b,n),%,随机模拟解非线性规划,debug=1;,a=0;,%,试验点下界,b=10;,%,试验点上界,n=1000;,%,试验点个数,r1=unifrnd(a,b,n,1);,%,a,b,均匀分布随机数矩阵,r2=unifrnd(a,b,n,1);,sol=r1(1)r2(1);,z0=inf;,for,i=1:n,x1=r1(i);,x2=r2(i);,lpc=lpconst(x1 x2),;,if,lpc=1,z=mylp(x1 x2),;,if,zz0,z0=z,;,sol=x1 x2,;,end,end,end,与,Monte Carlo,方法相似,但理论基础不同的方法,“,拟蒙特卡罗方法”,(Quasi-Monte Carlo,方法,),近年来也获得迅速发展。
这种方法的基本思想是“用确定性的超均匀分布序列,(,数学上称为,Low Discrepancy Sequences),代替,Monte Carlo,方法中的随机数序列对某些问题该方法的实际速度一般可比,Monte Carlo,方法提高数百倍。