第5讲蒙特卡洛方法的应用

实验目的实验内容学习如何应用蒙特卡洛方法解决实际问题1、起源和发展2、原理3、计算机模拟应用实例4、实验作业蒙特卡洛方法的应用一、MC的起源和发展•随机模拟方法，也称为MonteCarlo方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。冯·诺伊曼是公理化方法和计算机体系的领袖人物，MonteCarlo方法也是他的功劳。•事实上，MonteCarlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。18世纪下半叶的法国学者Buffon提出用投针试验的方法来确定圆周率π的值。这个著名的Buffon试验是MonteCarlo方法的最早的尝试！•历史上曾有几位学者相继做过这样的试验。不过呢，他们的试验是费时费力的，同时精度不够高，实施起来也很困难。然而，随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验就可以了。•MonteCarlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的。Buffon试验•假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为的针()，•则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离均匀的分布在区间上，针与平行线的夹角（不管相交与否）均匀的分布在区间上。•因此，针与线相交的充要条件是l1l21,0x,021sinxBuffon试验•从而针线相交的概率为•根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以估计出针线相交的概率，从而得到的估计值。•针与线的位置关系：ldxdlXPpl22sin2ˆ0sin20pfunctionpiguji=buffon(llength,mm)%llength是针的长度%mm是随机实验次数frq=0;xrandnum=unifrnd(0,0.5,1,mm);phi=unifrnd(0,pi,1,mm);forii=1:mmif(xrandnum(1,ii)=(llength*sin(phi(1,ii))/2))frq=frq+1;endendpiguji=2*llength/(frq/mm)buffon(.6,1000)piguji=3.1662buffon(.6,10000)piguji=3.1072buffon(.6,100000)piguji=3.1522buffon(.6,1000000)piguji=3.1386buffon(.6,1000000)piguji=3.1451buffon(.6,1000000)piguji=3.1418buffon(.6,1000000)piguji=3.1448buffon(.6,1000000)piguji=3.1405buffon(.6,1000000)piguji=3.1394二、MC的原理•应用MonteCarlo方法求解工程技术问题可以分为两类：•确定性问题•随机性问题思路1、针对实际问题建立一个简单且便于实现的概率统计模型，使问题的解对应于该模型中随机变量的概率分布或其某些数字特征，比如，均值和方差等。所构造的模型在主要特征参量方面要与实际问题或系统相一致的。2、根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，再进行随机模拟试验。3、根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。4、按照所建立的模型进行仿真试验、计算，求出问题的随机解。5、统计分析模拟试验结果，给出问题的估计以及其精度估计。6、必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。•收敛性:由大数定律,Monte-Carlo模拟的收敛是以概率而言的．•误差:用频率估计概率时误差的估计，可由中心极限定理，给定置信水平的条件下，有：••模拟次数:由误差公式得NU2/1||))((XgVar22/1)(UN三、MC的应用举例1、定积分的MC计算随机投点法样本平均值法几种降低估计方差的MC方法2、系统的可靠性数值模拟计算问题1、定积分的MC计算•事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。•下面考虑一个简单的定积分•为了说明问题，我们首先介绍两种求的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。dxxfba•在计算积分上，MC的实用场合是计算重积分•其中是维空间的点，当较大时，用MC方法比一般的数值方法有优点，主要是它的误差与维数无关。dPPgIkPmmm随机投点法方法简述：设，有限，，，并设是在上均匀分布的二维随机变量，其联合密度函数为。则易见是中曲线下方的面积。假设我们向中进行随机投点，若点落在下方，(即称为中的，否则称为不中，则点中的概率为。若我们进行了次投点，其中次中的，则用频率来估计概率。即。abMxf0Mybxayx0,:,YX,MybxaIabM0,1dxxfbaxfyxfyxfyabMpn0npnnp0•那么我们可以得到的一个估计•具体试验步骤为nnabM01ˆ用上式来估计的个数统计，和，计算，，，随机数，个独立地产生0,...,11,02nyxfxfMvyabuaxnivuUniiiiiiiii求解定积分的算例例计算定积分事实上，其精确解为用随机投点法求解：注增加样本数目，可提高计算精度，但计算时间也会提高。sj(0,4,4,1000000)result=7.23367.2432=0.4sin0.8dxx400.2cosfunctionresult=sj(a,b,m,mm)fq=0;x=unifrnd(a,b,1,mm);y=unifrnd(0,m,1,mm);forii=1:mmif(cos(x(1,ii))+2=y(1,ii))fq=fq+1;endendresult=fq*m*(b-a)/mm•注1随机投点法的思想简单明了，且每次投点结果服从二项分布，故，其中•注2可证是的无偏估计。若用估计的标准差来衡量其精度，则估计的精度的阶为。•注3这里，定积分的解，就对应我们选定的随机变量的概率值。pnbn,~0abMp1ˆ1ˆ21n2019/12/21例的计算1.单位圆的面积等于2.3.0.20.40.60.810.20.40.60.8141102dxx4)1(1102dxx用随机投点法求的值sjtdf_pi1(1000)piguji=3.0520sjtdf_pi1(10000)piguji=3.1204sjtdf_pi1(100000)piguji=3.1296functionpiguji=sjtdf_pi1(mm)%mm是随机实验次数frq=0;xrandnum=unifrnd(0,1,1,mm);yrandnum=unifrnd(0,1,1,mm);forii=1:mmifxrandnum(1,ii)^2+yrandnum(1,ii)^2=1frq=frq+1;endendpiguji=4*frq/mmsjtdf_pi2(100)piguji=3.2000sjtdf_pi2(1000)piguji=3.2120sjtdf_pi2(10000)piguji=3.1260sjtdf_pi2(100000)piguji=3.1373functionpiguji=sjtdf_pi2(mm)%mm是随机实验次数frq=0;xrandnum=unifrnd(0,1,1,mm);yrandnum=unifrnd(0,1,1,mm);forii=1:mmif(sqrt(1-xrandnum(1,ii)^2)=yrandnum(1,ii))frq=frq+1;endEnd,piguji=4*frq/mm样本平均值法•基本原理：对积分，设是上的一个密度函数，改写•可见，任一积分均可以表示为某个随机变量（函数）的期望。由矩法，若有个来自的观测值，则可给出的一个矩估计。•最简单的，若，有限，可取。dxxfbaxgba,XgXfEdxxgxgxfbanxgababxg/1•设是来自的随机数，则的一个估计为•具体步骤为•注可证是的无偏估计。一般而言，样本均值法要比随机投点法更有效。nxx,...,1baU,niiniiixfnabxgxfn1121ˆ用上式来估计，，，和计算，随机数个独立地产生nixfuabaxuuUniiin,...,1,...,1,012ˆ求解定积分的算例例计算定积分事实上，其精确解为样本平均值法求解：注增加样本数目，可提高计算精度，但计算时间也会提高。ybjzf(0,4,4,1000)result=7.3036ybjzf(0,4,4,10000)result=7.2970ybjzf(0,4,4,100000)result=7.2578dxx400.2cos7.2432=0.4sin0.8functionresult=ybjzf(a,b,m,mm)%a是积分的下限%b是积分的上限%积分函数cos(x)+2%mm是随机实验次数sum=0;xrandnum=unifrnd(a,b,1,mm);forii=1:mmsum=sum+cos(xrandnum(1,ii))+2;endresult=sum*(b-a)/mm2019/12/21例的计算1.单位圆的面积等于2.3.0.20.40.60.810.20.40.60.8141102dxx4)1(1102dxx用样本平均值法求的值functionresult=ybjzf1(a,b,m,mm)%a是积分的下限%b是积分的上限%积分函数%mm是随机实验次数xrandnum=unifrnd(a,b,1,mm);sum=sum(sqrt(1-xrandnum.^2));result=sum*(b-a)/mm;result=result*4ybjzf1(0,1,1,100)result=3.08745746887753ybjzf1(0,1,1,1000)result=3.15500646827616ybjzf1(0,1,1,10000)result=3.11911714237706ybjzf1(0,1,1,100000)result=3.14014654862983ybjzf1(0,1,1,1000000)result=3.14093979612119functionresult=ybjzf1(a,b,m,mm)%a是积分的下限%b是积分的上限%积分函数%mm是随机实验次数sum=0;xrandnum=unifrnd(a,b,1,mm);forii=1:mmsum=sum+1/(1+xrandnum(1,ii)^2));endresult=sum*(b-a)/mm;result=result*4ybjzf2(0,1,1,100)result=3.04500162146030ybjzf2(0,1,1,1000)result=3.14857120401090ybjzf2(0,1,1,10000)res