蒙特卡罗方法与计算机模拟资料

1、蒙特卡罗与计算机模拟,内容：计算机模拟(或称仿真)是一种广义数值计算 方法，适合解决一些规模大、难以解析化 以及受随机因素影响的不确定数学模型 目的：了解蒙特卡罗方法的基本思想，掌握利用 Matlab对离散/连续系统进行模拟的方法 要求：掌握Matlab随机数函数，处理应用问题 了解蒙特卡罗方法的起源和基本思想 掌握离散系统和连续系统计算机模拟实例 掌握随机函数 rand unifrnd normrnd exprnd 了解Matlab仿真模块 Simulink,蒙特卡罗方法的起源和基本思想,蒙特卡罗方法(Monte Carlo method)，或称计算机随机模拟方法，是一种基于“随机数”的计算

2、方法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用蒲丰投针的方法来计算圆周率，上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,蒲丰投针实验近似计算圆周率,蒲丰投针实验： 法国科学家蒲丰(Buffon)在1777年提 出的蒲丰投针实验是早期几何概率一个 非常著名

3、的例子。蒲丰投针实验的重要 性并非是为了求得比其它方法更精确的 值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导，由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花蒙特卡罗方法(MC) 蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l(la)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求),蒲丰投针实验近似计算圆周率,蒲丰投针实验： 如右图所示，以M 表示针落下后的中点， 以x表示M到最近一条 平行线的距离，以表示针与此线的交角： 针落地的所有可能结果满足： 其样本空间视作矩形区域,

4、面积是: 针与平行线相交的条件： 它是样本空间子集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,add); %初始化绘图方式为叠加 counter=0; n=2010; %初始化计数器和设定投针次数 x=unifrnd(0,a/

5、2,1,n); phi=unifrnd(0,pi,1,n); %样本空间 for i=1:n if x(i)l*sin(phi(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次,蒲丰投针实验近似计算圆周率,蒲丰投针实验的计算机模拟：,意大利数学家拉泽里尼得到了准确到6

6、位小数的值，不过他的实验因为太准确而受到了质疑,蒲丰投针实验计算圆周率,蒙特卡罗投点法是蒲丰投针实验的推广： 在一个边长为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);,常见分布的随机数产生语句,蒙特卡罗方法的关键步骤在于随机数的产生，计算机产生的随机数都不是真正的随机数(由算法确定

7、的缘故)，如果伪随机数能够通过一系列统计检验，我们也可以将其当作真正的随机数使用：,MATLAB随机数的“重置”问题,Matlab的随机数是伪随机数，但在一定的信度之下可以看作真正的随机数。问题是rand函数产生的随机数从一个随机数序列中取出来，而每次启动Matlab时，rand的状态都会被重置(相当于把序列的指针移到了随机数序列的开始)，换言之第一次启动Matlab调用的第n次rand函数与下一次启动调用的第n个rand函数产生相同的数值。 如果想打乱这种状态，可以为rand指定一个与当前时间相关的初始状态，而不用默认状态： rand(state,sum(100*clock); 或者 ran

8、d(state,sum(100*clock)*rand); 若使每次运行程序产生的值是相同的 rand(seed,0.1); rand(1),非常见分布的随机数的产生,1.逆变换方法,非常见分布的随机数的产生,举例说明： 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) %产生参

9、数 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(rFx); x(index(1)-1) %已编写通用离散分布随机数产生程序 scatrnd(x,px,n),2.Acceptance-Rejection 方法 最早由 Von Neumann提出，现在已经 广泛应用于

10、各种随机数的生成。,离散系统的计算机模拟实例,范例 报童售报问题:一报童每天清晨从邮局订购报纸后零售，每份报纸进价0.35元，售价0.5元，邮局要求最低订购数量为60份，根据过去经验一个报童一天平均售出报纸120份(且满足泊松分布)，未售出的报纸只要没有破损可退给邮局，试求报童每天清晨订购多少份报纸可获最大利润？ 1 数学建模,Pin=0.35元 Pout=0.5元 Nin60,200区间 Nout=poissrnd(120)随机数,离散系统的计算机模拟实例,2 计算机模拟 当订购 122 份报纸时，可获得最大利润 17.3584 元,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题,离

11、散系统的计算机模拟实例,范例 海港系统的卸载货物问题1 程序片段(船只到港时间均匀分布,船只卸货时间均匀分布) ShipBetweenTime(1)=unifrnd(15,145,1,1); %船只到港间隔时间随机化(均匀分布) ShipUnloadTime(1)=unifrnd(45,90,1,1); %船只卸货时间随机化(均匀分布) 通用程序haibor.m可实现多次模拟，并将结果保存到H.txt delete H.txt %清除历史数据 harbor(100,15,145,45,90) load H.txt;Hmean=mean(H); %导入H并按列取平均值,离散系统的计算机模拟实例,

12、范例 海港系统的卸载货物问题2 程序片段(船只到港时间指数分布,船只卸货时间均匀分布) ShipBetweenTime(1)=exprnd(60,1,1); %船只到港间隔时间随机化(指数分布) ShipUnloadTime(1)=unifrnd(45,90,1,1); %船只卸货时间随机化(均匀分布) 通用程序haibor2.m可实现多次模拟，结果保存到H2.txt delete H2.txt %清除历史数据 harbor2(100,60,45,90) load H2.txt;Hmean2=mean(H2); %导入H2并按列取平均值,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3

13、 程序片段(船只到港时间离散分布,船只卸货时间离散分布) 1 编写船只到港间隔离散累积分布函数并作阶梯图： xs=15:10:145; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;end px=0.009,0.029,0.035,0.051,0.090,0.161,0.200,0.172,0.125,0.071,0.037,0.017,0.003; Fx=0; for i=1:length(px),Fx=Fx,sum(px(1:i);end plot(10,x,Fx,-rs); hold on; stairs(0,x-5,145,Fx,1); set(

14、gca,xtick,0:5:145); set(gca,xgrid,on); axis tight;,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3 程序片段(船只到港时间离散分布,船只卸货时间离散分布) 2 编写船只到港间隔离散累积分布反函数并作线性插值： Fxi=0:0.001:1-eps; xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n); rnd=; for i=1:n index=find(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); en

15、d end %以上程序已编写通用M函数文件 harborrnd(xs,px,n) %即给出n个满足离散分布(x,px)的船只到港间隔随机数,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3 程序片段(船只到港时间离散分布,船只卸货时间离散分布) 3 编写船只卸货时间离散累积分布函数并作阶梯图： xs=45:5:90; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;end px=0.017,0.045,0.095,0.086,0.130,0.185,0.208,0.143,0.091; Fx=0; for i=1:length(px),Fx=Fx,

16、sum(px(1:i);end plot(40,x,Fx,-rs); hold on; stairs(40,x-2.5,90,Fx,1); set(gca,xtick,40:2.5:90); set(gca,xgrid,on); axis tight;,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3 程序片段(船只到港时间离散分布,船只卸货时间离散分布) 4 编写船只卸货时间离散累积分布反函数并作线性插值： Fxi=0:0.001:1-eps; xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n); rnd=; for i=1:n index=find

17、(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); end end %以上程序已编写通用M函数文件 harborrnd(xs,px,n) %即给出n个满足离散分布(x,px)的船只卸货时间随机数,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3 程序片段(船只到港时间指数分布,船只卸货时间均匀分布) 5 模拟船只到港间隔 / 卸货时间均为离散分布的海港系统 ShipBetweenTime(1)=harborrnd(sbtxs,sbtpx,1); %船只到港间隔时间随机化(离散分布) Ship

18、UnloadTime(1)=harborrnd(sutxs,sutpx,1); %船只卸货时间随机化(离散分布) 通用程序haibor3.m可实现多次模拟，结果保存到H3.txt delete H3.txt %清除历史数据 load harbor.mat %载入数据 harbor3(100,sbtxs,sbtpx,sutxs,sutpx) load H3.txt;Hmean3=mean(H3); %导入H3并按列取平均值,连续系统的计算机模拟实例,范例 某军导弹基地发现正北方向120km处有一艘敌舰以90km/h的速度向正东方向行驶，该基地即刻发射导弹进行拦击，导弹速率450km/h，制导系统

19、确保在任一时刻导弹都能对准敌舰 (问题1) 试问导弹何时何处击中敌舰,syms x y t dydt dxdt; solve(dydt/dxdt=(120-y)/(90*t-x),dydt2+dxdt2=4502,dydt,dxdt); ans=ans.dxdt,ans.dydt; dxdt=ans(2,1); dydt=ans(2,2); pretty(dxdt); pretty(dydt);,连续系统的计算机模拟实例,1 将随等距时间连续变化的状态变量轨迹x(t),y(t)用欧拉法离散化：,连续系统的计算机模拟实例,2 编写程序模拟导弹拦击敌舰过程 x1=0; y1=0; x2=0; y2

20、=120; t=0.001; v1=450; v2=90; dis=120; axis(0,40,0,140); grid on; set(gca,nextplot,add); for k=1:1000 x1=x1+v1*t*(v2*k*t-x1)/sqrt(v2*k*t-x1)2+(dis-y1)2); y1=y1+v1*t*(dis-y1)/sqrt(v2*k*t-x1)2+(dis-y1)2); x2=x2+v2*t; y2=y2; plot(x1,y1,ro,x2,y2,bs); frame(k)=getframe; if sqrt(x1-x2)2+(y1-y2)2)=0.1,brea

21、k;end end T=k*t,x1,y1 %微分方程求解和计算机模拟过程已整合进daodan1.m,连续系统的计算机模拟实例,3 导弹拦击敌舰过程模拟动画(AVI视频) warning off movie2avi(frame,daodan1,compression,Indeo5,quality,100,fps,6); %将模拟过程编码成AVI视频 导弹于0.2770小时在(25.0018千米, 119.9306千米)击中敌舰,连续系统的计算机模拟实例,(问题2) 如果敌舰即刻发现导弹，并以垂直导弹方向135km/h速度逃逸，试问导弹何时何处击中敌舰 1 建立模型,连续系统的计算机模拟实例,2

22、 求解并模拟模型 (下为片段,完整见daodan2.m) syms x1 y1 x2 y2 t dy1dt dx1dt dy2dt dx2dt; solve(dy1dt/dx1dt=(y2-y1)/(x2-x1),dy2dt/dx2dt=(x2-x1)/(y2-y1),dy1dt2+dx1dt2=4502,dy2dt2+dx2dt2=1352,dy1dt,dx1dt,dy2dt,dx2dt); for k=1:1000 x1=x1+dx1dt(x1,x2,y1,y2)*t; y1=y1+dy1dt(x1,x2,y1,y2)*t; x2=x2+dx2dt(x1,x2,y1,y2)*t; y2=y

23、2-dy2dt(x1,x2,y1,y2)*t; plot(x1,y1,ro,x2,y2,bs); frame(k)=getframe; if sqrt(x1-x2)2+(y1-y2)2)=0.3,break;end end,连续系统的计算机模拟实例,3 过程动画 (AVI视频) warning off movie2avi(frame,daodan2,compression,Indeo5,quality,100,fps,6); %将模拟过程编码成AVI视频 导弹于0.2660小时在(32.8916千米, 109.8963千米)击中敌舰,排队问题随机模拟,排队论主要研究随机服务系统的工作过程。,在

24、排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对随机服务现象的统计研究，找出反映这些随机现象平均特性的规律指标，如排队的长度、等待的时间及服务利用率。,1 系统的假设： （1） 顾客源是无穷的； （2） 排队的长度没有限制； （ 3）到达系统的顾客按先后顺序依次进入服务，“先到先服务”。,在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店设： 1顾客到来间隔时间服从参数为0.1的指数分布 对顾客的服务时间服从，上的均匀分布 排队按先到先服务规则，队长无限制,假定一个工作日为8小时，时间以分钟为单位。 1模拟一个工作日内完成服务的顾客个数及顾客平均等待时间t 2模拟100个工作日，求出平均每日完成服务的顾客个数及每日顾客的平均等待时间。,单服务员的排队模型模拟,w：总等待时间； ci：第i个顾客的到达时刻； bi：第i个顾客开始服务时刻； ei：第i个顾客服务结束时刻； xi：第i-1个顾客与第i个顾客之间到达的间隔时间； yi：对第i个顾客的服务时间。,符号说明,c1,b1,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,e4,b5,ci=ci-1+ xi ei=bi+yi bi=max(ci,ei-1),