蒙特卡洛算法.ppt

《蒙特卡洛算法.ppt》由会员分享，可在线阅读，更多相关《蒙特卡洛算法.ppt（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、蒙特卡洛算法一.蒙特卡洛算法的介绍二.随机点的产生1.伪随机数的产生2.准随机数的产生3.随机点产生程序分析蒙特卡洛算法的介绍算法简介蒙特卡洛算法，也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要数值计算方法，蒙特卡罗方法在金融工程学，宏观经济学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。蒙特卡洛算法的介绍算法描述以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。例如pi的计算：在一个面积为一的正方形里面画

2、一个圆，半径0.5，与正方形四边相切，在正方形中生成一系列随机点，统计单位圆内的点数与总点数，（圆面积和正方形面积之比为pi/4=m/n），当随机点取得越多（但即使取10的9次方个随机点时，其结果也仅在前4位与圆周率吻合）时，其结果越准确。随机点的产生伪随机算法随机数生成算法是一类重要的算法，广泛应用于仿真技术等场合。伪随机数是由确定的算法生成的，其分布函数与相关性均能通过统计测试。与真实随机数的差别在于，它们是由算法产生的，而不是一个真实的随机过程。一般地，伪随机数的生成方法很多，有线性同余法，直接法，逆转法等随机点的产生C/C+语言中伪随机数生成算法实际上是采用了“线性同余法“。具体的计算

3、如下：Xi=(Xi-1*A+C)modM其中A,C,M都是常数（一般会取质数）。当C=0时，叫做乘同余法。引出一个概念叫seed，它会被作为X0被代入上式中，然后每次调用rand()函数都会用上一次产生的随机值来生成新的随机值。可以看出实际上用rand()函数生成的是一个递推的序列，一切值都来源于最初的seed。所以当初始的seed取一样的时候，得到的序列都相同。C语言里面有RAND_MAX这样一个宏，定义了rand()所能得到的随机值的范围。在C里可以看到RAND_MAX被定义成0 x7fff，也就是32767。rand()函数里递推式中M的值就是32767。线性同余法生成的是伪随机数，粗略

4、符合均匀分布。随机点的产生准随机算法伪随机算法都存在差异性，不均匀性。因此，不要求新的发生器模拟真实的均匀分布，而力求任意大小的样本（尤其是小样本）都能满足低差异性。换言之，以牺牲随机性为代价，换来均匀性的提高，称其为准随机模拟器。目前有3种准随机序列可用来辅助生成均匀分布随机数，分别是Halton序列、Sobol序列、Latin超立方体序列。随机点的产生Halton序列以质数为基底产生的序列假设是以质数2为基底的Halton序列，范围在01之间。则从产生的列为1/21/43/41/85/83/87/81/169/16.怎么产生的呢？3=1+25=1+47=3+49=1+8.假设以质数3为基底

5、，Halton序列如下1/32/31/94/97/92/95/98/91/2710/2719/27.随机点的产生在二维中，01之间产生的点的序列就是（1/2，1/3）（1/4，2/3）（3/4，1/9）（1/8，4/9）（5/8，7/9）（3/8，2/9）（7/8，5/9）（1/16，8/9）（9/16，1/27）.核心代码分析vprivatestaticclassHaltonSequence/Halton序列的产生vstaticfinalintP=2,3;/以2、3为基底产生序列vstaticfinalintK=63,40;/使2和3产生的序列中元素的个数相同vprivatelonginde

6、x;/索引项vprivatedoublex;/x数组定义vprivatedoubleq;/二维数组q定义vprivateintd;/二维数组d定义vHaltonSequence(longstartindex)/输入参数值，得到序列中的元素vindex=startindex;vx=newdoubleK.length;/K.length=2vq=newdoubleK.length;vd=newintK.length;vfor(inti=0;iK.length;i+)/给q，d赋值vqi=newdoubleKi;/q0=63,q1=40vdi=newintKi;/d0=63,d1=40vvfor(i

7、nti=0;iK.length;i+)vlongk=index;/把参数传给kvxi=0;vvfor(intj=0;jKi;j+)/这部分通过运算来演示vqij=(j=0?1.0:qij-1)/Pi;vdij=(int)(k%Pi);vk=(k-dij)/Pi;vxi+=dij*qij;vvvvdoublenextPoint()/根据上面求得的点来产生第二个点vindex+;vfor(inti=0;iK.length;i+)vfor(intj=0;jKi;j+)vdij+;vxi+=qij;vif(dijPi)vbreak;vvdij=0;vxi-=(j=0?1.0:qij-1);vvvret

8、urnx;vvv改写后的程序（C+版）v#includev#includevusingnamespacestd;vintP2=2,3;vintK2=63,40;vlongindex;vdoublex2;vdoubleq2100;vintd2100;vvoidHaltonSequence()vvindex=1;vfor(inti=0;i2;i+)vvlongk=index;vxi=0;vvfor(intj=0;jKi;j+)vvqij=(j=0?1.0:qij-1)/Pi;vdij=(int)(k%Pi);vk=(k-dij)/Pi;vxi+=dij*qij;vvvvvoidnextPoint()vvindex+;vfor(inti=0;i2;i+)vvfor(intj=0;jKi;j+)vvdij+;vxi+=qij;vif(dijPi)vvbreak;vvdij=0;vxi-=(j=0?1.0:qij-1);vvvvvoidmain()vvinti;vHaltonSequence();vcout(x0,x1)endl;vfor(i=0;i500;i+)vvnextPoint();vcout(x0,x1)endl;vvv500个点的产生