第3讲蒙特卡洛方法初探.ppt
随机数的产生及蒙特卡洛随机模拟方法,实验目的,实验内容,学习随机数的产生及蒙特卡洛随机模拟方法的基本过程与方法。,1、数学模拟的方法。,4、实验作业。,3、蒙特卡洛随机模拟实例。,2、产生随机数的计算机命令。,数学模拟的方法,在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。,在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。,计算机模拟可以反复进行，改变系统的结构和系数都比较容易。,蒙特卡洛（Monte Carlo）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的统计分析，求得所研究系统的某些参数,用蒙特卡洛方法进行计算机模拟的步骤：,1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象,一、产生模拟随机数的计算机命令,在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：,1产生m*n阶(a，b)均匀分布U（a，b）的随机数矩阵：unifrnd(a,b,m,n);产生一个a，b均匀分布的随机数：unifrnd(a,b),当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它。,2产生mm*nn阶离散均匀分布的随机数矩阵：R=unidrnd(N)R=unidrnd(N,mm,nn),当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。,若连续型随机变量X的概率密度函数为其中 0为常数，则称X服从参数为 的指数分布。,指数分布的期望值为,排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。,例 顾客到达某商店的间隔时间服从参数为10(分钟)的指数分布(指数分布的均值为10),-指两个顾客到达商店的平均间隔时间是10分钟.即平均10分钟到达1个顾客.顾客到达的间隔时间可用exprnd(10)模拟。,设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为其中 0为常数，则称X服从参数为 的泊松分布。,泊松分布在排队系统、产品检验、天文、物理等领域有广泛应用。,泊松分布的期望值为,6产生1个参数为n,p的二项分布的随机数binornd(n,p)，产生mn个参数为n,p的二项分布的随机数binornd(n,p,m,n)。,掷一枚均匀硬币，正面朝上的次数X服从参数为，p的二项分布,XB(1,p),总结：常见分布的随机数产生语句,补充：随机数的产生命令MATLAB可以直接产生满足各种分布的随机数具体命令如下:产生mn阶0,1上均匀分布的随机数矩阵 rand(m,n)产生一个0,1上均匀分布的随机数 rand 产生mn阶a,b上均匀分布的随机数矩阵 unifrnd(a,b,m,n)产生一个a,b上均匀分布的随机数 unifrnd(a,b)产生一个1:n的随机排列(元素均出现且不重复)p=randperm(n)注意:randperm(6)与unifrnd(1,6,1,6)的区别,产生mn阶均值为mu方差为sigma的正态分布的随机数矩阵 normrnd(mu,sigma,m,n)产生一个均值为mu方差为sigma的正态分布的随机数 normrnd(mu,sigma)产生mn阶期望值为mu(mu=1/)的指数分布的随机数矩阵 exprnd(mu,m,n)产生一个期望值为mu的指数分布的随机数 exprnd(mu)注意:产生一个参数为的指数分布的随机数应输入 exprnd(1/),产生mn阶参数为A1,A2,A3的指定分布name的随机数矩阵 random(name,A1,A2,A3,m,n)产生一个参数为为A1,A2,A3的指定分布name的随机数 random(name,A1,A2,A3)举例:产生24阶的均值为0方差为1的正态分布的随机数矩阵 random(Normal,0,1,2,4)name的取值可以是(详情参见help random)：norm or Normal/unif or Uniformpoiss or Poisson/beta or Betaexp or Exponential/gam or Gammageo or Geometric/unid or Discrete Uniform,二、频率的稳定性模拟,1.事件的频率 在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率 f=m/n,2.频率的稳定性,掷一枚均匀硬币，记录掷硬币试验中频率P*的波动情况。,function liti1(p,mm)pro=zeros(1,mm);randnum=binornd(1,p,1,mm)a=0;for i=1:mm a=a+randnum(1,i);pro(i)=a/i;end pro=pronum=1:mm;plot(num,pro),在Matlab中编辑.m文件输入以下命令：,在Matlab命令行中输入以下命令：,liti1(0.5,1000),在Matlab命令行中输入以下命令：,liti1(0.5,10000),练习掷一枚不均匀硬币，正面出现概率为0.3，记录前1000次掷硬币试验中正面频率的波动情况，并画图。,在Matlab命令行中输入以下命令：,liti1(0.3,1000),例2 掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况，并画图。,在Matlab中编辑.m文件输入以下命令：,function liti2(p,mm)pro=zeros(1,mm);randnum=binornd(1,p,2,mm)；a=0;for i=1:mm a=a+randnum(1,i)*randnum(2,i);pro(i)=a/i;end pro=pro，num=1:mm;plot(num,pro),liti2(0.4,100),liti2(0.4,10000),二、几何概率模拟,1.定义,向任一可度量区域G内投一点，如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比，而与g的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。,2.概率计算,P(A)=A的度量/S的度量,例5 两人约定于12点到1点到某地会面，先到者等20分钟后离去，试求两人能会面的概率？,解：设x,y分别为甲、乙到达时刻(分钟),令A=两人能会面=(x,y)|x-y|20，x60,y60,P(A)=A的面积/S的面积=（602-402）/602=5/9=0.5556,function proguji=liti5(mm)%mm 是随机实验次数frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;proguji=0;for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1;endendproguji=frq/mm,liti5(10000)proguji=0.5557,三、蒲丰投针实验：法国科学家蒲丰(Buffon)在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子。蒲丰投针实验的重要性并非是为了求得比其它方法更精确的值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导，由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花蒙特卡罗方法(MC)蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l(la)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求),蒲丰投针实验：如右图所示，以M表示针落下后的中点，以x表示M到最近一条平行线的距离，以表示针与此线的交角：针落地的所有可能结果满足：其样本空间视作矩形区域,面积是:针与平行线相交的条件：它是样本空间子集A，面积是：积分计算 syms l phi;int(l/2*sin(phi),phi,0,pi);%ans=l因此，针与平行线相交的概率为：从而有：特别当 时 p为统计频率,蒲丰投针实验的计算机模拟：format long;a=1;l=0.6;%显示精度,线宽和针长figure;axis(0,pi,0,a/2);%初始化绘图板set(gca,nextplot,add);%初始化绘图方式为叠加counter=0;n=2010;%初始化计数器和设定投针次数x=unifrnd(0,a/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.);frame(i)=getframe;%描点并取帧 title(Current Point,num2str(i),Total,num2str(n);counter=counter+1;%统计针与线相交的次数 endendfren=counter/n;pihat=2*l/(a*fren)%用频率近似计算%movie(frame,1)%播放帧动画1次,蒲丰投针实验计算圆周率,蒙特卡罗投点法是蒲丰投针实验的推广：在一个边长为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);,x,y,o,例 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的概率能毁伤敌人一门火炮，有1/6的概率能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来，利用频率稳定性，确定有效射击(毁伤一门炮或全部消灭)的概率.,五、复杂概率模拟（选学）,分析：这是一个复杂概率问题，可以通过理论计算得到相应的概率.,为了直观地显示我方射击的过程，现采用模拟的方式。,需要模拟出以下两件事：,1.问题分析,1 观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是1/2 因此，可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确，反之为不正确,2 当指示正确时，我方火力单位的射击结果情况,模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷骰子的方法来确定：如果出现的是、三个点：则认为没能击中敌人；如果出现的是、点：则认为毁伤敌人一门火炮；若出现的是点：则认为毁伤敌人两门火炮,2.符号假设,i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击(毁伤一门炮或两门炮)的概率,3.模拟框图,function liti6(p,mm)efreq=zeros(1,mm);randnum1=binornd(1,p,1,mm);randnum2=unidrnd(6,1,mm);k1=0;k2=0;k3=0;for i=1:mm if randnum1(i)=0 k1=k1+1;else if randnum2(i)=3 k1=k1+1;elseif randnum2(i)=6 k3=k3+1;else k2=k2+1;end end efreq(i)=(k2+k3)/i;end num=1:mm;plot(num,efreq),在Matlab中编辑.m文件输入以下命令：,在Matlab命令行中输入以下命令：,liti6(0.5,2000),在Matlab命令行中输入以下命令：,liti6(0.5,20000),5.理论计算,6.结果比较,模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程,作业：1.掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况，并画图。,2:两船欲停靠同一个码头,设两船到达码头的时间各不相干，而且到达码头的时间在一昼夜内是等可能的.如果两船到达码头后需在码头停留的时间分别是1 小时与2 小 时，试求在一昼夜内，任一船到达时，需要等待空出码头的概率.（频率估计概率）,function liti2(p,mm)pro=zeros(1,mm);randnum=binornd(1,p,2,mm)；a=0;for i=1:mm a=a+randnum(1,i)*randnum(2,i);pro(i)=a/i;end pro=pro，num=1:mm;plot(num,pro),在Matlab中编辑.m文件输入以下命令：,习题1解答：,liti2(0.4,100),liti2(0.4,10000),3.某厂生产的灯泡能用1000小时的概率为0.8,能用1500小时的概率为0.4,求已用1000小时的灯泡能用到1500小时的概率（频率估计概率）。,2.在一袋中有10 个相同的球，分别标有号码1,2,10。今有放回任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率（频率估计概率）,掷三枚不均匀硬币，每枚正面出现概率为0.3，记录前1000次掷硬币试验中至少两枚都为正面频率的波动情况，并画图。,选做作业:,4:两船欲停靠同一个码头,设两船到达码头的时间各不相干，而且到达码头的时间在一昼夜内是等可能的.如果两船到达码头后需在码头停留的时间分别是1 小时与2 小 时，试求在一昼夜内，任一船到达时，需要等待空出码头的概率.（频率估计概率）,5:在0,1,2,3,.,9中不重复地任取4个数，求它们能排成首位非零的四位偶数的概率.（频率估计概率）,6:从1,2,.,10十个数字中有放回地任取5个数字,求取出的5个数字中按由小到大排列,中间的那个数等于 4 的概率.（频率估计概率）,
