蒙特卡罗方法介绍及其建模应用Part II 2012-07-07

Monte-Carlo方法介绍及其建模应用朱连华Tel：13675122648南京信息工程大学数学与统计学院E-mail:[e-mail]2020/4/1423:39南京信息工程大学课程说明公用邮箱：[e-mail]：ahualian2008参考书目:黄燕、吴平.SAS统计分析及应用，机械工业出版社.陈杰.Matlab宝典，电子工业出版社.张文彤等.SPSS11.0统计分析教程，北京希望电子出版社.薛益、陈立萍.统计建模与R软件，清华大学出版社.2020/4/1423:39南京信息工程大学主要内容蒙特卡洛方法应用实例2排队论模拟介绍3蒙特卡洛方法介绍12009-B眼科病床安排应用42020/4/1423:39南京信息工程大学蒙特卡洛方法应用实例概率计算模拟分析1定积分的MC计算2系统可靠性模拟计算32020/4/1423:39南京信息工程大学概率计算模拟分析12020/4/1423:39南京信息工程大学频率的稳定性模拟频率:–在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数，频率f=m/n频率的稳定性：f----P–例1：掷一枚均匀硬币，记录掷硬币试验中频率P的波动情况–functionliti21(p,mm)–pro=zeros(1,mm);–randnum=binornd(1,p,1,mm)–a=0;–fori=1:mm–a=a+randnum(1,i);–pro(i)=a/i;–end–pro=pro–num=1:mm;–plot(num,pro)2020/4/1423:39南京信息工程大学liti21(0.4,100)liti21(0.4,10000）liti21(0.5,1000)liti21(0.5,10000)2020/4/1423:39南京信息工程大学例1'：掷一枚不均匀硬币，正面出现概率为0.3，记录前1000次掷硬币试验中正面频率的波动情况liti21(0.3,1000)2020/4/1423:39南京信息工程大学例2:掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况functionliti22(p,mm)pro=zeros(1,mm);randnum=binornd(1,p,2,mm)；a=0;fori=1:mma=a+randnum(1,i)*randnum(2,i);pro(i)=a/i;endpro=pro，num=1:mm;plot(num,pro)2020/4/1423:39南京信息工程大学古典概率模拟例3:在一袋中有10个相同的球，分别标有号码1,2,…,10。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取3个球，求这3个球的号码均为偶数的概率。（用频率估计概率）解：有放回取3个球,所有取法有103种;有放回取3个偶数号码的球,所有取法有53种.所以81105)(33APfunctionproguji=liti23(n,mm)frq=0;randnum=unidrnd(n,mm,3);proguji=0;fori=1:mma=(randnum(i,1)+1)*(randnum(i,2)+1)*(randnum(i,3)+1)；ifmod(a,2)==1frq=frq+1endend;proguji=frq/mm2020/4/1423:39南京信息工程大学古典概率模拟例4:两盒火柴，每盒20根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有5根火柴的概率有多大？（用频率估计概率）functionproguji=liti24_0(mm)%mm是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*20);proguji=0;fori=1:mma1=0;a2=0;j=1;while(a120)&(a220)ifrandnum(i,j)==1a1=a1+1;elsea2=a2+1;endj=j+1;endifabs(a1-a2)=5frq=frq+1;endend;proguji=frq/mmliti24_0(100)proguji=0.4800liti24_0(1000)proguji=0.4970liti24_0(10000)proguji=0.4910liti24_0(100000)proguji=0.49842020/4/1423:39南京信息工程大学古典概率模拟例4':两盒火柴，每盒n根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有k根火柴的概率有多大？（用频率估计概率）functionproguji=liti24_1(n,k,mm)%n是每盒中的火柴数%k是剩余的火柴数%mm是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*n);proguji=0;fori=1:mma1=0;a2=0;j=1;while(a1n)&(a2n)ifrandnum(i,j)==1a1=a1+1;elsea2=a2+1;endj=j+1;endifabs(a1-a2)=k,frq=frq+1;end%a1=a1,a2=a2,frq%pauseend;proguji=frq/mmliti24_1(20,5,100)proguji=0.4800liti24_1(20,5,1000)proguji=0.4970liti24_2(20,5,10000)proguji=0.4910liti4(20,5,100000)proguji=0.49842020/4/1423:39南京信息工程大学几何概率模拟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的面积=（602-402）/602=5/9=0.55562020/4/1423:39南京信息工程大学functionproguji=liti25(mm)%mm是随机实验次数frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;proguji=0;forii=1:mmifabs(randnum(ii,1))=20frq=frq+1;endendproguji=frq/mm几何概率模拟liti25(10000)proguji=0.55572020/4/1423:39南京信息工程大学复杂概率模拟例6:在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的概率能毁伤敌人一门火炮，有1/6的概率能全部消灭敌人．现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来，利用频率稳定性，确定有效射击(毁伤一门炮或全部消灭)的概率.2020/4/1423:39南京信息工程大学复杂概率模拟分析：–这是一个复杂概率问题，可以通过理论计算得到相应的概率.–为了直观地显示我方射击的过程，现采用模拟的方式。1.问题分析需要模拟出以下两件事：[1]观察所对目标的指示正确与否–模拟试验有两种结果，每一种结果出现的概率都是1/2．–因此，可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确，反之为不正确．2020/4/1423:39南京信息工程大学复杂概率模拟[2]当指示正确时，我方火力单位的射击结果情况–模拟试验有三种结果：•毁伤一门火炮的可能性为1/3(即2/6)•毁伤两门的可能性为1/6•没能毁伤敌火炮的可能性为1/2(即3/6)–这时可用投掷骰子的方法来确定：•如果出现的是１、２、３三个点：则认为没能击中敌人；•如果出现的是４、５点：则认为毁伤敌人一门火炮；•若出现的是６点：则认为毁伤敌人两门火炮．2020/4/1423:39南京信息工程大学复杂概率模拟2.符号假设–i：要模拟的打击次数；–k1：没击中敌人火炮的射击总数；–k2：击中敌人一门火炮的射击总数；–k3：击中敌人两门火炮的射击总数．–E：有效射击(毁伤一门炮或两门炮)的概率2020/4/1423:39南京信息工程大学复杂概率模拟3.模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜mm?E=(k2+k3)/mm停止硬币正面?YNNY1，2，34，562020/4/1423:39南京信息工程大学复杂概率模拟functionliti26(p,mm)efreq=zeros(1,mm);randnum1=binornd(1,p,1,mm);randnum2=unidrnd(6,1,mm);k1=0;k2=0;k3=0;fori=1:mmifrandnum1(i)==0k1=k1+1;elseifrandnum2(i)=3k1=k1+1;elseifrandnum2(i)==6k3=k3+1;elsek2=k2+1;endendefreq(i)=(k2+k3)/i;endnum=1:mm;plot(num,efreq)2020/4/1423:39南京信息工程大学复杂概率模拟liti26(0.5,2000)liti26(0.5,20000)2020/4/1423:39南京信息工程大学复杂概率模拟5.理论计算模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程．设：观察所对目标指示正确确观察所对目标指示不正10jA0：射中敌方火炮的事件；A1：射中敌方一门火炮的事件；A2：射中敌方两门火炮的事件．则由全概率公式：P(A0)=P(j=0)P(A0∣j=0)+P(j=1)P(A0∣j=1)=25.02121021P(A1)=P(j=0)P(A1∣j=0)+P(j=1)P(A1∣j=1)=613121021P(A2)=P(j=0)P(A2∣j=0)+P(j=1)P(A2∣j=1)=1216121021E=25.0121612020/4/1423:39南京信息工程大学定积分的MC计算22020/4/1423:39南京信息工程大学定积分的MC计算事实上，不少的统计问题最后都归结为定积分的近似计算问题!相对于其它方法，用MC方法比一般的数值方法有优点，主要体现在它的误差与维数m无关！下面考虑一个简单的定积分为了说明问题，我们首先介绍两种求的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。dxxfba2020/4/1423:39南京信息工程大学方法简述：设a，b有限，0f(x)M，={(x,y):axb，0yM},并设(X,Y)是在上均匀分布的二维随机变量，其联合密度函数为MybxaIabM0,1则易见是中y=f(x)曲线下方的面积假设我们向中进行随机投点，则点落在y=f(x)下方的概率p,abMdxxfabMdxdyabMdxdyIabMXfYPpbabaxfMybxaxfy)(1]1[1)}({)(00,)(dxxfba随机投点法2020/4/1423:39南京信息工程大学若我们进行了n次投点，其中n0次点落入y=f(x)曲线下方，则用频率n0/n来估计概率p。即