(精品)蒙特卡洛方法初探

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,蒙特卡洛方法初探,实验目的,实验内容,学习计算机模拟的基本过程与方法。,1,、模拟的概念。,4,、实验作业,。,3,、计算机模拟实例。,2,、产生随机数的计算机命令。,模拟的概念,模拟,就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。,模拟的基本思想,是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息,模拟的方法,1,、,物理模拟,：,对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等。,物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等。,在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。,这时，计算机模拟几乎成为唯一的选择。,在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。,2,、,数学模拟,计算机模拟可以反复进行，改变系统的结构和系数都比较容易。,蒙特卡洛（,Monte Carlo,）,方法,是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的统计分析，求得所研究系统的某些参数,用蒙特卡洛方法进行计算机模拟的步骤,：,1,设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。,2,模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象,产生模拟随机数的计算机命令,在,Matlab,软件中，可以直接产生满足各种分布的随机数，命令如下：,1,产生,m,*,n,阶,(a,，,b),均匀分布,U,（,a,，,b,）的随机数矩阵：,unifrnd,(,a,b,m,n);,产生一个,a,，,b,均匀分布的随机数：,u,nifrnd,(,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,),，,产生,m,n,个参数为,n,p,的二项分布的随机数,binornd(n,p,m,n,),。,掷一枚均匀硬币，正面朝上的次数,X,服从参数为，,p,的二项分布,XB(1,p),频率的稳定性模拟,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=,pro,num=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),例,3:,在一袋中有,10,个相同的球，分别标有号码,1,2,10,。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取,3,个球，求这,3,个球的号码均为偶数的概率。,（用频率估计概率）,解：有放回取,3,个球,所有取法有,10,3,种,;,有放回取,3,个偶数号码的球,所有取法有,5,3,种,.,所以,古典概率模拟,function,proguji,=liti3(n,mm),frq,=0;,randnum,=unidrnd(n,mm,3);proguji=0;,for i=1:mm a=(randnum(i,1)+1)*(randnum(i,2)+1)*(randnum(i,3)+1),；,if mod(a,2)=1,frq,=frq+1,end,end;,proguji,=,frq,/mm,例,4,两盒火柴，每盒,20,根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有,5,根火柴的概率有多大？（用频率估计概率）,function,proguji,=liti40(mm),%mm,是随机实验次数,frq,=0;,randnum,=binornd(1,0.5,mm,2*20);proguji=0;,for i=1:mm,a1=0;a2=0;j=1;,while(a120)&(a2=5,frq,=frq+1;end,end;,proguji,=,frq,/mm,liti40(100),proguji=0.4800,liti40(1000),proguji=0.4970,liti40(10000),proguji=0.4910,liti40(100000),proguji=0.4984,例,4,两盒火柴，每盒,n,根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有,k,根火柴的概率有多大？（用频率估计概率）,function,proguji,=liti4(n,k,mm),%n,是每盒中的火柴数,%k,是剩余的火柴数,%mm,是随机实验次数,frq,=0;,randnum,=binornd(1,0.5,mm,2*,n);proguji,=0;,for i=1:mm,a1=0;a2=0;j=1;,while(a1n)&(a2=k ,frq,=frq+1;end,%a1=a1,a2=a2,frq%pause,end;,proguji,=,frq,/mm,liti4(20,5,100),proguji=0.4800,liti4(20,5,1000),proguji=0.4970,liti4(20,5,10000),proguji=0.4910,liti4(20,5,100000),proguji=0.4984,几何概率模拟,1.,定义,向任一可度量区域,G,内投一点，如果所投的点落在,G,中任意可度量区域,g,内的可能性与,g,的度量成正比，而与,g,的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。,2.,概率计算,P,(,A,)=,A,的度量,/,S,的度量,例,5,两人约定于,12,点到,1,点到某地会面，先到者等,20,分钟后离去，试求两人能会面的概率？,解：设,x,y,分别为甲、乙到达时刻,(,分钟,),令,A,=,两人能会面,=(x,y)|x-y|,20,，,x,60,y,60,P(A)=A,的面积,/S,的面积,=,（,60,2,-40,2,）,/60,2,=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;,end,end,proguji,=,frq,/mm,liti5(10000),proguji,=0.5557,例,6,在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有,50,是准确的，而我方火力单位，在指示正确时，有,1/3,的概率能毁伤敌人一门火炮，有,1/6,的概率能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的,1,次打击结果显现出来，利用频率稳定性，确定有效射击,(,毁伤一门炮或全部消灭,),的概率,.,复杂概率模拟,分析,：,这是一个复杂概率问题，可以通过理论计算得到相应的概率,.,为了直观地显示我方射击的过程，现采用模拟的方式。,需要模拟出以下两件事：,1.,问题分析,1,观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是,1/2,因此，,可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确，反之为不正确,2,当指示正确时，我方火力单位的射击结果情况,模拟试验有三种结果：毁伤一门火炮的可能性为,1/3(,即,2/6),，毁伤两门的可能性为,1/6,，没能毁伤敌火炮的可能性为,1/2(,即,3/6),这时,可用投掷骰子的方法来确定,：,如果出现的是、三个点：则认为没能击中敌人；,如果出现的是、点：则认为毁伤敌人一门火炮；,若出现的是点：则认为毁伤敌人两门火炮,2.,符号假设,i,：,要模拟的打击次数；,k,1,：,没击中敌人火炮的射击总数；,k,2,：,击中敌人一门火炮的射击总数；,k,3,：,击中敌人两门火炮的射击总数,E,：有效射击,(,毁伤一门炮或两门炮,),的概率,3.,模拟框图,初始化,:i=0,k,1,=0,k,2,=0,k,3,=0,i=i+1,骰子,点数,?,k,1,=k,1,+1,k,2,=k,2,+1,k,3,=k,3,+1,k,1,=k,1,+1,i,mm?,E=(k,2,+k,3,)/mm,停止,硬币正面,?,Y,N,N,Y,1,，,2,，,3,4,，,5,6,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 ran