MonteCarlo模拟教程.pptx

Monte Carlo Monte Carlo 模拟模拟内容提纲1.引言2.Monte Carlo模拟基本思想3.随机数生成函数4.应用实例举例5.排队论模拟6.Monte Carlo模拟求解规划问题引言引言(Introduction)(Introduction)Monte Carlo方法：方法：蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte Carlo名字的由来：名字的由来：Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博闻名的首都，该城以赌博闻名Nicholas Metropolis(1915-1999)Monte-Carlo,MonacoMonte CarloMonte Carlo方法的基本思想方法的基本思想 蒙特卡罗方法，或称计算机随机模拟方法，是一种基蒙特卡罗方法，或称计算机随机模拟方法，是一种基于于“随机数随机数”的计算方法。源于美国在第的计算方法。源于美国在第二二次世界大战研次世界大战研制原子弹的制原子弹的“曼哈顿计划曼哈顿计划”，该计划的主持人之一数学家该计划的主持人之一数学家冯冯诺伊曼用驰名世界的赌城诺伊曼用驰名世界的赌城摩纳哥的摩纳哥的Monte CarloMonte Carlo来命名这种方法，为它蒙上了一层神秘色彩。来命名这种方法，为它蒙上了一层神秘色彩。蒙特卡罗蒙特卡罗方法的基本思想很早以前就被人们所发现和方法的基本思想很早以前就被人们所发现和利用。早在利用。早在1717世纪，人们就知道用事件发生的世纪，人们就知道用事件发生的“频率频率”来来决定事件的决定事件的“概率概率”。1 19 9世纪人们用世纪人们用蒲丰蒲丰投针的方法来投针的方法来计计算算圆周率圆周率，上，上世纪世纪4040年代电子计算机的出现，特别是近年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。大量、快速地模拟这样的试验成为可能。蒲丰投针实验蒲丰投针实验：法国科学家蒲丰法国科学家蒲丰(Buffon)(Buffon)在在17771777年提出的蒲丰年提出的蒲丰投针实验是早期几何概率一个非常著名的例子。投针实验是早期几何概率一个非常著名的例子。蒲丰蒲丰投针实验投针实验的重要性并非是为了求得比其它方法更精确的重要性并非是为了求得比其它方法更精确的的值，而是它开创了使用随机数处理确定性数学问值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导。题的先河，是用偶然性方法去解决确定性计算的前导。由此由此可以可以领略到从领略到从“概率土壤概率土壤”上开出的一朵瑰丽的鲜上开出的一朵瑰丽的鲜花花-蒙特卡罗蒙特卡罗方法方法(MC)MC)蒲丰投针实验可归结为下面的数学问题：平面上画蒲丰投针实验可归结为下面的数学问题：平面上画有距离为有距离为a a的一些平行线，向平面上任意投一根长为的一些平行线，向平面上任意投一根长为l l(la)(la)的针，假设针落在任意位置的可能性相同，试求针的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率与平行线相交的概率P(P(从而求从而求)蒲丰投针实验蒲丰投针实验：如右图所示，以如右图所示，以M M表示针落下表示针落下后的中点，以后的中点，以x x表示表示M M到最近一条平行到最近一条平行线的距离，以线的距离，以表示针与此线的交角：表示针与此线的交角：针落地的所有可能结果满足：针落地的所有可能结果满足：其样本空间视作矩形区域其样本空间视作矩形区域,面积是面积是:针与平行线相交的条件：针与平行线相交的条件：它是样本空间它是样本空间子集子集A A，面积是：，面积是：syms l phi;int(l/2*sin(phi),phi,0,pi)%ans=l因此，针与平行线相交的概率为：因此，针与平行线相交的概率为：从而有：从而有：蒲丰投针实验蒲丰投针实验的计算机模拟：的计算机模拟：format long;%format long;%设置设置1515位显示精度位显示精度a=1;a=1;l=0.6;l=0.6;%两平行线间的宽度两平行线间的宽度和和针长针长figure;axis(0,pi,0,a/2);%figure;axis(0,pi,0,a/2);%初始化绘图板初始化绘图板set(gca,nextplot,set(gca,nextplot,a adddd););%初始化绘图方式为叠加初始化绘图方式为叠加counter=0;counter=0;n=n=20102010;%初始化计数器和设定初始化计数器和设定投针次数投针次数x=unifrnd(0,a/2,1,n);x=unifrnd(0,a/2,1,n);p phihi=unifrnd(0,pi,1,n);=unifrnd(0,pi,1,n);%样本空间样本空间for i=1:nfor i=1:n if x(i)l*sin(if x(i)l*sin(p phihi(i)/2(i)/2%满足此条件表示针与线的相交满足此条件表示针与线的相交 plot(phi(i),x(i),plot(phi(i),x(i),r.r.);counter=counter+1;counter=counter+1;%统计针与线相交的次数统计针与线相交的次数frame(frame(countercounter)=getframe;%)=getframe;%描点并取帧描点并取帧 end endendendfren=counter/n;fren=counter/n;pihat=2*l/(a*fren)pihat=2*l/(a*fren)%用用频率频率近似计算近似计算figure(2)figure(2)movie(frame,1)%movie(frame,1)%播放帧动画播放帧动画1 1次次一些人进行了实验，其结果列于下表一些人进行了实验，其结果列于下表 ：实验者年份投计次数的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.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;endenddisp(投点法近似计算的投点法近似计算的为:,num2str(4*m/n);xyo(a/2,a/2)基本思想基本思想 由上面的例子可以看出，当所求问题的解由上面的例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数是某个事件的概率，或者是某个随机变量的数学期望，或者是与之有关的量时，通过某种试学期望，或者是与之有关的量时，通过某种试验的方法，得出该事件发生的频率，再通过它验的方法，得出该事件发生的频率，再通过它得到问题的解。这就是蒙特卡罗方法的基本思得到问题的解。这就是蒙特卡罗方法的基本思想。想。蒙特卡罗方法的关键步骤在于随机数的产蒙特卡罗方法的关键步骤在于随机数的产生，计算机产生的随机数都不是真正的随机数生，计算机产生的随机数都不是真正的随机数(由算法确定的缘故由算法确定的缘故)，如果伪随机数能够通过，如果伪随机数能够通过一系列统计检验，我们也可以将其当作真正的一系列统计检验，我们也可以将其当作真正的随机数使用。随机数使用。rand(seed,0.1);rand(1)每次运行程序产生的值是相同的每次运行程序产生的值是相同的rand(state,sum(100*clock)*rand);rand(1)%每次重新启动每次重新启动matlab时，输出的随机数不一样时，输出的随机数不一样注意注意:产生一个参数为产生一个参数为的指数分布的随机数应输入的指数分布的随机数应输入 exprnd(1/exprnd(1/)产生产生mnmn阶阶参数为参数为A1,A2,A3A1,A2,A3的指定分布的指定分布namename的的随机数矩阵随机数矩阵 random(name,A1,A2,A3,m,n)random(name,A1,A2,A3,m,n)举例举例:产生产生2 24 4阶的均值为阶的均值为0 0方差为方差为1 1的正态分布的随机数矩阵的正态分布的随机数矩阵 random(Normal,0,1,2,4)random(Normal,0,1,2,4)namename 的取值可以是的取值可以是(详情参见详情参见help random)help random)：norm or Normal/unif or Uniformnorm or Normal/unif or Uniformpoiss or Poisson/beta or Betapoiss or Poisson/beta or Betaexp or Exponential/gam or Gammaexp or Exponential/gam or Gammageo or Geometric/unid or Discrete Uniformgeo 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为为rFxinv=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);endr=rand;index=find(rMAXK或PMAXP时停止迭代框框 图图初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数xj=aj+R(bj-aj)j=1,2,nj=0j=j+1,p=p+1PMAXP?YNxj=aj+R(bj-aj)gi(X)0?i=1,2nYNjMAXK?YN输出X,Q,停止YN 在Matlab软件包中编程，共需三个文件:randlp.m,mylp.m,lpconst.m.主程序为randlp.m.%mylp.m%mylp.mfunction z=mylp(x)%目标函数z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2);%转化为求最小值问题%randlp.m%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 endend 与Monte Carlo方法相似，但理论基础不同的方法“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)近年来也获得迅速发展。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提高数百倍。