Monte_Carlo方法第八讲

2020/3/141数学建模专题之蒙特卡罗方法主讲：张建侠时间：2012-08-09数学建模专题之MonteCarlo方法2020/3/142内容提纲1.引言2.MonteCarlo模拟基本思想3.随机数生成函数4.应用实例举例5.随机游走模拟6.机械零件的可靠度计算7.排队论模拟8.MonteCarlo方法求解规划问题9.MonteCarlo方法预测搜救区域数学建模专题之MonteCarlo方法2020/3/143引言(Introduction)MonteCarlo方法：蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来：•MonteCarlo是摩纳哥（monaco)最大的城市，该城以赌博闻名。Monte-Carlo,Monaco数学建模专题之MonteCarlo方法2020/3/144引言(Introduction)MonteCarlo方法的起源：二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用(中子的链锁反应)。但其基本思想并非新颖，可追溯到18世纪后半叶的蒲丰（Buffon）随机投针实验（1777年）。NicholasMetropolis(1915-1999)JohnVonNeumann(1903-1957)数学建模专题之MonteCarlo方法2020/3/145引言(Introduction)MonteCarlo方法的应用：物理：核物理，热力学与统计物理，粒子输运问题等数学：多重积分、解微分方程、非线性方程组求解等工程领域：真空技术，水力学，激光技术等经济学领域：期权定价、项目管理、投资风险决策等其他领域：化学、医学，生物，生产管理、系统科学、公用事业等方面。随着科学技术的发展，其应用范围将更加广泛。数学建模专题之MonteCarlo方法2020/3/146MonteCarlo方法的基本思想Buffon投针实验1768年，法国数学家ComtedeBuffon利用投针实验估计的值2,()LpLdddL数学建模专题之MonteCarlo方法2020/3/147SolutionThepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:[0,)[0,)AdTherandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.Theprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegralsin1002sinldlpPAldAdd数学建模专题之MonteCarlo方法2020/3/148例1.蒲丰投针问题利用关系式：求出π值其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。2lpd2212()lllNdpdpdn数学建模专题之MonteCarlo方法2020/3/149一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929数学建模专题之MonteCarlo方法2020/3/1410基本思想MS基本思想：为了求解数学、物理、工程技术或随机服务系统等方面的问题，首先构造一个模型（概率统计模型），使所求问题的解正好是该模型的参数或特征量或有关量，然后通过模拟（统计试验），给出模型参数或特征量的估计值，最后得出所求问题的近似解。MS特点：1.方法新颖、应用面广、实用性强；2.随机模拟方面的算法简单，但计算量大；3.模拟结果具有随机性，精度较低；4.模拟结果的收敛过程服从概率规律。数学建模专题之MonteCarlo方法2020/3/1411例1在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.为了能显示我方20次射击的过程，现采用模拟的方式。举例数学建模专题之MonteCarlo方法2020/3/1412需要模拟出以下两件事：[2]当指示正确时，我方火力单位的射击结果情况[1]观察所对目标的指示正确与否模拟试验有两种结果，每一种结果出现的概率都是1/2．因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确．模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)．这时可用投掷骰子的方法来确定：如果出现的是１、２、３三个点：则认为没能击中敌人；如果出现的是４、５点：则认为毁伤敌人一门火炮；若出现的是６点：则认为毁伤敌人两门火炮．问题分析数学建模专题之MonteCarlo方法2020/3/1413i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数；E：有效射击比率；E1：20次射击平均每次毁伤敌人的火炮数．符号说明数学建模专题之MonteCarlo方法2020/3/1414模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=(0*k1+1*k2+2*k3)/20停止硬币正面?YNNY1，2，34，56数学建模专题之MonteCarlo方法2020/3/1415模拟结果消灭敌人火炮数试验序号投硬币结　果指示正确指　示不正确掷骰子结　果０１２１正∨４∨２正∨４∨３反∨∨４正∨１∨５正∨２∨６反∨∨７正∨３∨８正∨６∨９反∨∨１０反∨∨数学建模专题之MonteCarlo方法2020/3/1416消灭敌人火炮数试验序号投硬币结果指示正确指示不正确掷骰子结果０１２１１正∨２∨１２反∨∨１３正∨３∨１４反∨∨１５正∨６∨１６正∨４∨１７正∨２∨１８正∨４∨１９反∨∨２０正∨６∨从以上模拟结果可计算出：170.35(0131423/200.520EE）数学建模专题之MonteCarlo方法2020/3/1417理论计算设：01j观察所对目标指示不正确观察所对目标指示正确A0：射中敌方火炮的事件；A1：射中敌方一门火炮的事件；A2：射中敌方两门火炮的事件．则由全概率公式：E=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E1=33.01212611数学建模专题之MonteCarlo方法2020/3/1418结果比较理论计算和模拟结果的比较分类项目无效射击有效射击平均值模拟0.650.350.5理论0.750.250.33虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程．要使结果接近理论计算值，必须加大模拟次数，这就要求使用计算机模拟了。用蒙特卡洛方法进行计算机模拟的步骤：[1]设计一个逻辑框图，即建立模拟模型．[2]根据流程图编写程序，模拟随机现象．可通过生成具有各种概率分布的随机数来模拟随机现象，进行模拟试验．[3]分析模拟结果，给出所求问题的近似解（求解）.数学建模专题之MonteCarlo方法2020/3/1419注：rand(n)=rand(n,n)Matlab中的随机数生成函数randperm(m)生成一个由1:m组成的随机排列randn(m,n)生成一个满足正态mn的随机矩阵rand(m,n)生成一个满足均匀分布的mn随机矩阵，矩阵的每个元素都在(0,1)之间。perms(1:n)生成由1:n组成的全排列，共n!个结果数学建模专题之MonteCarlo方法2020/3/1420name的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......random('name',A1,A2,A3,M,N)Matlab中的随机数生成函数数学建模专题之MonteCarlo方法2020/3/1421fix(x):截尾取整，直接将小数部分舍去floor(x):不超过x的最大整数ceil(x):不小于x的最小整数round(x):四舍五入取整Matlab中的取整函数数学建模专题之MonteCarlo方法2020/3/1422模拟随机投掷均匀硬币，验证国徽朝上与朝下的概率是否都是1/2n=10000;%给定试验次数m=0;%m表示试验成功（国徽朝上）的次数fori=1:nx=randperm(2)-1;%randperm(2）生成1和2的随机排列y=x(1);ify==0%0表示国徽朝上，1表示国徽朝下m=m+1;endendfprintf('国徽朝上的频率为：%f\n',m/n);小实例一：投掷硬币数学建模专题之MonteCarlo方法2020/3/1423随机投掷骰子，验证各点出现的概率是否为1/6n=10000;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;fori=1:nx=randperm(6);y=x(1);switchycase1,m1=m1+1;case2,m2=m2+1;case3,m3=m3+1;case4,m4=m4+1;case5,m5=m5+1;otherwise,m6=m6+1;endend...%输出结果小实例二：投掷骰子数学建模专题之MonteCarlo方法2020/3/1424用蒙特卡罗(MonteCarlo)投点法计算的值n=10000;a=2;m=0;%m表示落入圆内的次数fori=1:nx=rand(1)*a/2;y=rand(1)*a/2;if(x^2+y^2=(a/2)^2)m=m+1;endendfprintf('计算出来的pi为：%f\n',4*m/n);小实例三：蒙特卡罗投点法数学建模专题之MonteCarlo方法2020/3/1425小实例三：蒙特卡罗投点法ezplot('x^2+y^2-1',[-1.11.1]);holdonaxisequalplot([-1-111-1],[-111-1-1]);N=0;fork=1:100N_point=10000;xy=(rand(2,N_point)-0.5)*2;d=sqrt(xy(1,:).^2+xy(2,:).^2);N(k)=l