7.3--随机模拟和统计分析

模块7MATLAB的应用随机模拟和统计分析7.3随机模拟(MonteCarlo算法)蒙特卡罗（MonteCarlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是为了适应当时原子能事业的发展而发展起来的。传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。该方法是20世纪40年代中期美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777年，法国数学家布丰（GeorgesLouisLecleredeBuffon，1707—1788）提出用投针实验的方法求圆周率π。这被认为是蒙特卡罗方法的起源。蒙特·卡罗方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。1.随机模拟概念2.随机模拟原理当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：（1）构造或描述概率过程；（2）实现从已知概率分布抽样；（3）建立各种估计量。3.随机模拟步骤（1）构造或描述概率过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。（2）实现从已知概率分布抽样构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。3.随机模拟步骤另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。（3）建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。4.随机模拟应用通常蒙特卡罗模拟通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗模拟是一种有效的求出数值解的方法。一般蒙特卡罗模拟在数学中最常见的应用就是蒙特卡罗积分。借助计算机技术，蒙特卡罗模拟实现了两大优点：一是简单，省却了繁复的数学报导和演算过程，使得一般人也能够理解和掌握；二是快速。简单和快速，是蒙特卡罗方法在现代项目管理中获得应用的技术基础。下面举三个例子例1分子模拟计算使用蒙特·卡罗方法进行分子模拟计算是按照以下步骤进行的：（1）使用随机数发生器产生一个随机的分子构型。（2）对此分子构型的其中粒子坐标做无规则的改变，产生一个新的分子构型。（3）计算新的分子构型的能量。（4）比较新的分子构型于改变前的分子构型的能量变化，判断是否接受该构型。若新的分子构型能量低于原分子构型的能量，则接受新的构型，使用这个构型重复再做下一次迭代。若新的分子构型能量高于原分子构型的能量，则计算玻尔兹曼因子，并产生一个随机数。若这个随机数大于所计算出的玻尔兹曼因子，则放弃这个构型，重新计算。若这个随机数小于所计算出的玻尔兹曼因子，则接受这个构型，使用这个构型重复再做下一次迭代。（5）如此进行迭代计算，直至最后搜索出低于所给能量条件的分子构型结束。例2项目管理项目管理中蒙特·卡罗模拟方法的一般步骤是：（1）对每一项活动，输入最小、最大和最可能估计数据，并为其选择一种合适的先验分布模型；（2）计算机根据上述输入，利用给定的某种规则，快速实施充分大量的随机抽样（3）对随机抽样的数据进行必要的数学计算，求出结果（4）对求出的结果进行统计学处理，求出最小值、最大值以及数学期望值和单位标准偏差（5）根据求出的统计学处理数据，让计算机自动生成概率分布曲线和累积概率曲线(通常是基于正态分布的概率累积S曲线)（6）依据累积概率曲线进行项目风险分析。在力学中，蒙特卡罗方法多被用来求解稀薄气体动力学问题，其中最为成功的是澳大利亚G.A.伯德等人发展的直接模拟统计试验法。此法通过在计算机上追踪几千个或更多的模拟分子的运动、碰撞及其与壁面的相互作用，以模拟真实气体的流动。它的基本假设与玻耳兹曼方程一致，但它是通过追踪有限个分子的空间位置和速度来代替计算真实气体中分布函数。模拟的相似条件是流动的克努曾数（Kn）相等，即数密度与碰撞截面之积保持常数。对每个分子分配以记录其位置和速度的单元。在模拟过程中分别考虑分子的运动和碰撞，在此平均碰撞时间间隔内，分别计算分子无碰撞的运动和典型碰撞。若空间网格取得足够小，其中任意两个分子都可以互相碰撞。具体决定哪两个刚体分子相撞，是随机取一对分子，计算它们的相对速度，根据此值与最大相对速度的比值和随机取样比较的结果，来决定该对分子是否入选。碰撞后分子的速度根据特定分子模型的碰撞力学和随机取样决定。分子与壁面碰撞后的速度，则根据特定的反射模型和随机取样决定。对于运动分子的位置和速度的追踪和求矩可以得出气体的密度、温度、速度等一些感兴趣的宏观参量。而对于分子与壁面间的动量和能量交换的记录则给出阻力、举力和热交换系数等的数学期望值。例3力学5.随机模拟(MonteCarlo)法应用实例例4蒙特卡洛法估计圆周率分析：设元的半径为1，则其外切正方形的面积为4，元的面积为pi,现在模拟产生在正方形ABCD中的均匀分布的点n个，如果这n个点中有m个在该园内，则圆的面积与正方形的面积之比可近似为m/n;即nmnm44用for循环编程：参考代码：n=input('请输入产生的点的个数：');m=0;fori=1:nx=-1+2*rand;y=-1+2*rand;ifx^2+y^2=1m=m+1;endendmy_pi=4*m/n请输入产生的点的个数：10000my_pi=3.1372请输入产生的点的个数：1000my_pi=3.0880请输入产生的点的个数：10000000000my_pi=3.1416用关系运算代替for循环：参考代码：n=input('请输入产生的点的个数：');x=unifrnd(-1,1,n,1);y=unifrnd(-1,1,n,1);m=sum((x.^2+y.^2)=1);my_pi=4*m/n请输入产生的点的个数：100000my_pi=3.1319请输入产生的点的个数：1000000my_pi=3.1426请输入产生的点的个数：10000000000错误使用rand内存不足。请键入HELPMEMORY查看选项。例5蒙特卡洛法计算积分根据几何意义，该积分是以f(x,y)为曲面顶，以A为底的柱体C的体积。用下列简单思路来求I的近似值：假设C被包在几何体D的内部，D的体积容易获得，若在D内部产生1个均匀分布的随机数，那么P(随机数落在C内)≈C的体积与D的体积之比现产生在D内的n个均匀分布的随机数，若其中m个在C内部，那么可以用下式近似计算定积分II≈D的体积*m/nAyxyxfdxdyyxfIA),(0),(,),(，若例5蒙特卡洛法计算积分12221yxdxdyx分析：显然该几何体C在立方体D:的内部，立方体D的体积为4.而C是参考代码：n=input('请输入产生的点的个数：');x=unifrnd(-1,1,n,1);y=unifrnd(-1,1,n,1);z=rand(n,1);%在积分区域内，z的最大值为1，最小值为0m=sum(((x.^2+y.^2)=1)&(x.^2+z.^2)1);I=4*m/n请输入产生的点的个数：10000I=2.661610,1,1zyxzzxyx0,1,12222symsxyint(int(sqrt(1-x^2),y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)%用符号计算的精确解来验证一下ans=8/3f=@(x,y)sqrt(1-x.^2);ymin=@(x)-sqrt(1-x.^2);ymax=@(x)sqrt(1-x.^2);I=integral2(f,-1,1,ymin,ymax)%用数值计算的近似解来验证一下I=2.6667例5蒙特卡洛法计算积分12221yxdxdyx练习：蒙特卡洛法计算积分其中区域D：围成.Ddxdyxy2xyxy22与练习参考答案：1.区域D的两条边界交点为(-1,1)和(4,2)；被积函数在D上的最大值为16，ezplot('y^2=x')holdonezplot('y=x-2')axistightgridonaxis([0,4,-1,2])该积分对应几何体位于下面的立方体(体积=192)之内160,21,40:,,zyxzyx参考代码1：k=input('请输入实验轮次=');%多轮试验，多次模拟n=input('请输入每轮实验次数n=');%每次模拟取n个数x=4*rand(n,k);%生成取值于（0,4）上的随机数y=unifrnd(-1,2,n,k);%生成取值于（-1,2）上的随机数z=16*rand(n,k);%生成取值于（0,16）上的随机数V=192;%大立方体体积V=4*3*16m=sum((x=y+2)&(x=y.^2)&(z=x.*y.^2));%落入要求的几何体内的随机满足的条件，这个需要自己分析I=mean(V*m/n)请输入实验轮次=10请输入每轮实验次数n=100000I=7.6212