第3讲蒙特卡洛方法基本思想课件.ppt

1、实验目的实验目的实验内容实验内容学习计算机模拟的基本过程与方法。学习计算机模拟的基本过程与方法。1 1、模拟的概念。、模拟的概念。4 4、实验作业。、实验作业。3 3、计算机模拟实例。、计算机模拟实例。2 2、产生随机数的计算机命令。、产生随机数的计算机命令。模拟的概念模拟的概念 模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。 模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息模拟的方法模拟的方法1、物理模拟物理模拟： 对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船

2、艇实验、沙盘作业等。 物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等。 在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。 在一定的假设条件下，运用数学运算模拟系统在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。计算机上进行的，称为计算机模拟。2、数学模拟数学模拟 计算机模拟可以反复进

3、行，改变系统的结构和系数都比较容易。 蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数 蒙特卡洛方法也称为蒙特卡洛方法也称为随机模拟方法,其起源最早可以追溯到18世纪下半叶的Buffon试验.例 在1777年,法国学者Buffon提出用试验方法求圆周率鸬闹.其原理如下：假设平面上有元数条距离为1的等矩平行线,现向该平面随机地投掷一根长度为KI1的针,则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解z针的中心点与最近的平行线间的距

4、离Z均匀地分布在区间0.1/2上,针与平行线的夹角以不管相交与否)均匀地分布在区间0,而上(见图6。.于是,针与线相交的充要条件是本寸,从而针线相交概率为1用蒙特卡洛方法进行计算机模拟的步骤用蒙特卡洛方法进行计算机模拟的步骤：1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象产生模拟随机数的计算机命令产生模拟随机数的计算机命令 在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：2产生m*n阶离散均匀分布的随机数矩阵：R = unidrnd(N)R = unidrnd(N,mm)R = u

5、nidrnd(N,mm,nn) 当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它。1产生m*n阶a，b均匀分布U（a，b）的随机数矩阵： unifrnd (a,b,m, n)产生一个a，b均匀分布的随机数：unifrnd (a,b)当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。机械加工得到的零件尺寸的偏差、射击命中点与目标的偏差、各种测量误差、人的身高、体重等，都可近似看成服从正态分布。若连续型随机变量X的概率密度函数为 其中 0为常数，则称X服从

6、参数为 的指数分指数分布布。0001)(/xxexft指数分布的期望值为 排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。注意：注意：Matlab中，产生参数为 的指数分布的命令为exprnd( )例例 顾客到达某商店的间隔时间服从参数为顾客到达某商店的间隔时间服从参数为1010的指的指数分布数分布 指数分布的均值为指数分布的均值为1010。 指两个顾客到达商店的平均间隔时间是指两个顾客到达商店的平均间隔时间是1010个单个单位时间位时间. .即平均即平均1010个单位时间到达个单位时间到达1 1个顾客个顾客. .

7、顾客顾客到达的间隔时间可用到达的间隔时间可用exprnd(10)exprnd(10)模拟。模拟。设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为其中 0为常数，则称X服从参数为 的帕松分布帕松分布。, 2 , 1 , 0,!)(kkekXPk帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。帕松分布的期望值为1 1 事件的频率事件的频率 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。 频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 掷一枚均匀硬币，记

8、录掷硬币试验中频率掷一枚均匀硬币，记录掷硬币试验中频率P P* *的波动情况。的波动情况。 R = binornd(N,P,mm,nn) 例例1 1 频率的稳定性频率的稳定性3 3 概率的频率定义概率的频率定义 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，记次试验，记m m是是n n次次试验中事件试验中事件A A发生的次数。当试验次数发生的次数。当试验次数n n很大时，很大时，如果频率如果频率m/nm/n稳定地在某数值稳定地在某数值p p附近摆动，而且一附近摆动，而且一般地说，随着试验次数的增加，这种摆动的幅度般地说，随着试验次数的增加，这种摆动的幅度越来越小，称数值越来

9、越小，称数值p p为事件为事件A A在这一组不变的条件在这一组不变的条件下发生的概率，记作下发生的概率，记作P(A)=p.P(A)=p.4 频率的基本性质频率的基本性质 （1） 对任意事件对任意事件A，有，有 1)(0AP（2）1)(SP0)(P（3）若）若A1，A2，An是互不相容的，则是互不相容的，则 )()(11nkknkkAPAP 频率定义的意义频率定义的意义:(1) 提供了估计概率的方法提供了估计概率的方法;(2)提供了一种检验理论正确与否的准则提供了一种检验理论正确与否的准则.理论依据：理论依据： 大数定律大数定律 大量的随机现象中平均结果的稳定性大量的随机现象中平均结果的稳定性

10、大数定律的客观背景大数定律的客观背景大量抛掷硬币大量抛掷硬币正面出现频率正面出现频率字母使用频率字母使用频率生产过程中的生产过程中的废品率废品率大数定律大数定律贝努里（Bernoulli） 大数定律设 nA 是 n 次独立重复试验中事件 A 发生的次数, p 是每次试验中 A 发生的概率，则0有0limpnnPAn或1limpnnPAn在概率的统计定义中，事件 A 发生的频率“ 稳定于”事件 A 在一次试验中发生的概率是指：nnAnnA频率与 p 有较大偏差pnnA是小概率事件, 因而在 n 足够大时, 可以用频率近似代替 p . 这种稳定称为依概率稳定.贝努里（贝努里（Bernoulli）

11、大数定律的意义大数定律的意义：定义a 是一常数，0limaYPnn(或)1limaYPnn则称随机变量序列,21nYYY依概率收敛于常数 a , 记作aYnPn故pnnnPA,21nYYY是一系列随机变量，设0有若在 Bernoulli 定理的证明过程中， Y n 是相互独立的服从 0-1分布的随机变量序列 Xk 的算术平均值, Y n 依概率收敛于其数学期望 p . 结果同样适用于服从其它分布的独立随机变量序列Chebyshev 大数定律,21nXXX相互独立，设随机变量序列(指任意给定 n 1, 相互独立)，且具有相同的数学期望和方差nXXX,21, 2 , 1,)(,)(2kXDXEkk

12、则0有01lim1nkknXnP或11lim1nkknXnP定理的意义定理的意义：当 n 足够大时，算术平均值几乎就是一个常数,可以用算术平均值近似地代替数学期望.具有相同数学期望和方差的独立随机变量序列的算术平均值依概率收敛于数学期望. 例如要估计某地区的平均亩产量，要例如要估计某地区的平均亩产量，要收割某些有代表性的地块，例如收割某些有代表性的地块，例如n n 块块. . 计算其平均亩产量，则当计算其平均亩产量，则当n n 较大时，可用较大时，可用它作为整个地区平均亩产量的一个估计它作为整个地区平均亩产量的一个估计. .辛钦大数定律辛钦大数定律 设,21nXXX相互独立，服从同一分布，且具

13、有数学期望 E(X k) = , k= 1,2,则对任意正数 001lim1nkknXnP,21nXXX相互独立，注3: 设随机变量序列, 2 , 1,)(iXEkki则0有01lim1knikinXnP具有相同的分布，且记knikiMXn1111nPA),(21kAAAgnP),(21kg则则22nPAknPkA),(21kxxxg连续，若 大数定律以严格的数学形式表达了随大数定律以严格的数学形式表达了随机现象最根本的性质之一：机现象最根本的性质之一：它是随机现象统计规律的具体表现它是随机现象统计规律的具体表现.大数定律在理论和实际中都有广泛的应用大数定律在理论和实际中都有广泛的应用.平均结

14、果的稳定性平均结果的稳定性例例1 1 频率的稳定性频率的稳定性1 1 事件的频率事件的频率 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。 频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 掷一枚均匀硬币，记录掷硬币试验中频率掷一枚均匀硬币，记录掷硬币试验中频率P P* *的波动情况。的波动情况。 R = binornd(N,P,mm,nn) function liti1(n,p,mm)pro=zeros(1,mm);randnum = binornd(n,p,1,mm)a=

15、0;for i=1:mm a=a+randnum(1,i); pro(i)=a/i;end pro=pronum=1:mm;plot(num,pro)在在MatlabMatlab中编辑中编辑.m.m文件输入以下命令：文件输入以下命令：在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.5,1000)在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.5,10000)练习练习 频率的稳定性频率的稳定性1 1 事件的频率事件的频率 R = binornd(N,P,mm,nn) 在一组不变的条件下，重复作在一组不变的条件

16、下，重复作n n次试验，次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。 频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 练习练习掷一枚不均匀硬币，正面出现概率为掷一枚不均匀硬币，正面出现概率为0.30.3，记录前记录前10001000次掷硬币试验中正面频率的波动次掷硬币试验中正面频率的波动情况，并画图。情况，并画图。 在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.3,1000)例例2 2 掷两枚不均匀硬币，每枚正面出现概率掷两枚不均匀硬币，每枚正面出现概率为为0.40.4，记录前，记录前100

17、01000次掷硬币试验中两枚都为次掷硬币试验中两枚都为正面频率的波动情况，并画图。正面频率的波动情况，并画图。 在在Matlab中编辑中编辑.m文件输入以下命令：文件输入以下命令：function liti2(n,p,mm)pro=zeros(1,mm);randnum = binornd(n,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) 熊宇乐 y=zeros(1,1000); a=binornd(1,0.4,1,1000);b=bino

18、rnd(1,0.4,1,1000); c=0;d=0; for i=1:1000 c=c+a(1,i).*b(1,i); y(i)=c/i; end y=y; num=1:1000; plot(num,y)孟亚function bino (n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);for i=1:m if x(i)=1&y(i)=1 s(i)=1; else s(i)=0; endend for i=1:m y(i)=sum(s(1,1:i)/i;endplot(y)liti2(1,0.4,100)liti2(1,0.4,10000)在一袋中有在一

19、袋中有10 10 个相同的球，分别标有号码个相同的球，分别标有号码1,2,1,2,10,10。每次任取一个球，记录其号码。每次任取一个球，记录其号码后不放回袋中，再任取下一个。这种取法叫后不放回袋中，再任取下一个。这种取法叫做做“不放回抽取不放回抽取”。今不放回抽取。今不放回抽取3 3个球，个球，求这求这3 3个球的号码均为偶数的概率。个球的号码均为偶数的概率。（用频率（用频率估计概率）估计概率） 例例3:解：令解：令A=不放回抽取不放回抽取3个球，求这个球，求这3个球的号码个球的号码均为偶数均为偶数=(2,4,6)，(2,4,8)，.，(6,8,10)121)(31035PPNkAPfunc

20、tion proguji=liti3(nn,num,mm)%nn 是每盒中的火柴数 %num 是剩余的火柴数%mm 是随机实验次数frq=0; randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a120)&(a2=5 frq=frq+1; end% a1=a1,a2=a2,frq % pauseend proguji=frq/mmfunction proguji=liti3(nn,num,mm)%nn 是每盒中的火柴数 %num 是剩余的火柴数%mm 是随机实验次数frq=0; randnum=b

21、inornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a1nn)&(a2=num frq=frq+1; end% a1=a1,a2=a2,frq % pauseend proguji=frq/mm例例4 4 两盒火柴，每盒两盒火柴，每盒2020根。每次随机在任一根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有完而另一盒中至少还有5 5根火柴的概率有多大？根火柴的概率有多大？（用频率估计概率）（用频率估计概率） liti4(20,5,100)proguji

22、 = 0.4800 liti4(20,5,1000)proguji = 0.4970 liti4(20,5,10000)proguji = 0.4910 liti4(20,5,100000)proguji = 0.4984function proguji=liti4(mm)%mm 是随机实验次数frq=0; randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a120)&(a2=5 frq=frq+1; end% a1=a1,a2=a2,frq % pauseend proguji=frq/mm二二

23、. 几何概率几何概率1.1.定义定义 向任一可度量区域向任一可度量区域G G内投一点，如果所投内投一点，如果所投的点落在的点落在G G中任意可度量区域中任意可度量区域g g内的可能内的可能性与性与g g的度量成正比，而与的度量成正比，而与g g的位置和形的位置和形状无关，则称这个随机试验为几何型随状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。机试验。或简称为几何概型。2. 概率计算概率计算 1. P(A)=A的度量的度量/S的度量的度量两人约定于两人约定于12点到点到1点到某地会面，先点到某地会面，先到者等到者等20分钟后离去，试求两人能会面分钟后离去，试求两人能会面的概率？的概

24、率？ 例例5:解：设解：设x,y分别为甲、乙到达时刻分别为甲、乙到达时刻(分钟分钟)令令A=两人能会面两人能会面=(x,y)|x-y|20，x60,60,y y6060P(A)=A的面积的面积/S的面积的面积=（602-402）/602=5/9=0.5556function 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

25、frq=frq+1; endendproguji=frq/mmliti5(10000)proguji = 0.5557例例2 2在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现，我方指挥所对敌方目标的指示有30是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的1打击结果显现出来，利用频率稳定性，确定有效射击的概率分析分析： 这是一个概率问题，可以通过理论计算得到相应的概率和期望值. 为了能显

26、示我方射击的过程，现采用模拟的方式。 需要模拟出以下两件事： 1. 1. 问题分析问题分析1 1 观察所对目标的指示正确与否观察所对目标的指示正确与否模拟试验有两种结果，每一种结果出现的概率都是1/2 因此，可用投掷一枚硬币的方式予以确可用投掷一枚硬币的方式予以确定定，当硬币出现正面时为指示正确，反之为不正确2 2 当指示正确时，我方火力单位的射击结当指示正确时，我方火力单位的射击结果情况果情况 模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6) 这时可用投掷骰子的方法来确定可用投掷骰子的方法来确定：如果出现的是、

27、三个点：则认为没能击中敌人；如果出现的是、点：则认为毁伤敌人一门火炮；若出现的是点：则认为毁伤敌人两门火炮2. 2. 符号假设符号假设i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击比率； 3. 3. 模拟框图模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1imm?E=(k2+k3)/mm 停止硬币正面?YNNY1，2，34，56function binomoni(p,mm)efreq=zeros(1,mm);randnum1

28、= 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)在在MatlabMatlab中编辑中编辑.m.m文件输入以下命令：文件输入以下命令：在在MatlabMatlab命令行中输入以下命令：命令行

29、中输入以下命令：liti2_moni(0.5,2000)在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti2_moni(0.5,20000)5. 5. 理论计算理论计算6. 6. 结果比较结果比较 模拟结果与理论计算近似一致，能更模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程加真实地表达实际战斗动态过程 3. 某厂生产的灯泡能用1000小时的概率为0.8, 能用1500小时的概率为0.4 , 求已用1000小时的灯泡能用到1500小时的概率2. 在一袋中有在一袋中有10 个相同的球，分别标有号码个相同的球，分别标有号码1,2,10。今任取两个球，求取

30、得的第一个球号码为奇数，。今任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率。第二个球的号码为偶数的概率。1. 掷三枚不均匀硬币，每枚正面出现概率为掷三枚不均匀硬币，每枚正面出现概率为0.3，记，记录前录前1000次掷硬币试验中至少两枚都为正面频率的次掷硬币试验中至少两枚都为正面频率的波动情况，并画图。波动情况，并画图。 作业作业: 某厂生产的灯泡能用1000小时的概率为0.8, 能用1500小时的概率为0.4 , 求已用1000小时的灯泡能用到1500小时的概率解解 令 A 灯泡能用到1000小时 B 灯泡能用到1500小时所求概率为)()(APABPABPAB218 . 04 . 0)()(APBP例例 精品课件精品课件！精品课件精品课件！例例:在一袋中有在一袋中有10 个相同的球，分别标有号个相同的球，分别标有号码码1,2,10。今任取两个球，求取得的第。今任取两个球，求取得的第一个球号码为奇数，第二个球的号码为一个球号码为奇数，第二个球的号码为偶数的概率。偶数的概率。解：令解：令A=A=抽取抽取2 2个球，第一个球号码为奇数，第二个球个球，第一个球号码为奇数，第二个球的号码为偶数的号码为偶数=(1,2)=(1,2)，(1,4)(1,4)，. .，(9,10)(9,10)1859025)(2101515PPPNkAP