《蒙特卡罗方法完整教程.docx》由会员分享，可在线阅读，更多相关《蒙特卡罗方法完整教程.docx（14页珍藏版）》请在课桌文档上搜索。
1、MonteCarlo方法法1概述MonteCarlo法不同于确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。MonteCarlo方法（MCM）,也称为统计试验方法,是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动试验）和位势理论，主要是争论匀称介质的稳定状态。它是用一系列随机数来近似解决问题的一种方法,是通过查找一个概率统计的相像体并用试验取样过程来获得该相像体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解inspirit更接近于物理试验结果，而不是经典数值计算结果。普遍认为我们当前所应用的MC技术，其
2、进展约可追溯至1944年，尽管在早些时候仍有很多未解决的实例。MCM的进展归功于核武器早期工作期间LOSAlamOS（美国我国试验室中子散射争论中心）的一批科学家。LosAlamoS小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓舞了MCM在各种问题中的应用RMjMonteCarlo”的名称取自于Monaco（摩纳哥）内以赌博消遣而著名的一座城市。MonteCarlo方法的应用有两种途径：仿真和取样。仿真是指供应实际随机现象的数学上的仿照的方法。一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动来仿照中子的锯齿形路径。取样是指通过争论少量的随机的子集来演绎大量元素的特性的方法。
3、例如，/（x）在vx。上的平均值可以通过间歇性随机选取的有限个数的点的平均值来进行估量。这就是数值积分的MonleCark）方法。MCM已被胜利地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略嬉戏，以及其它竞赛活动中都会消失。MonteCarlo计算方法需要有可得的、听从特定概率分布的、随机选取的数值序列。2随机数和随机变量的产生5一10全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数g（x）,使其将整数变
4、换为随机数。以某种方法选取与,并依据Z+1=g（xk）产生下一个随机数。最一般的方程g（x）具有如下形式：g（x）=（0r+c）mod机（1）其中=初始值或种子（X。0）乘法器（0）C=增值（c0）m二模数对于/数位的二进制整数，其模数通常为2、例如，对于31位的计算机机即可取231。这里x0,和C都是整数，且具有相同的取值范围相a,mc,mx0o所需的随机数序卜便可由下式得(2)+=(4x,+c)md7该序列称为固性闻余序闻例如，若XO=c=7且机=10,则该序列为7,6,9,0,7,6,9,0(3)可以证明，同余序列总会进入一个循环套；也就是说，最终总会消失一个无休止重复的数字的循环。(3
5、)式中序列周期长度为4。当然，一个有用的序列必是具有相对较长周期的序列。很多作者都用术语乘何M和据今闻M分别指代C=O和CHO时的线性同余法。选取%,c和加的法则可参见6,10。这里我们只关怀在区间(0,1)内听从匀称分布的随机数的产生。用字符U来表示这些数字，则由式(2)可得XU=q(4)m这样U仅在数组,lm,2町，(加一1)/帆中取值。(对于区间(0)内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。其它检测方法可参见3,6。)产生区间(,加内匀称分布的随机数X,可用下式X=a-(b-a)U(5)用计算机编码产生的随机数(采用式(2)和(4)并不是完全随机的；事实上，给定序列
6、种子，序列的全部数字U都是完全可猜测的。一些作者为强调这一点，将这种计算机产生的序列称为伪婚抗教但假如适当选取凡。和加，序列U的随机性便足以通过一系列的统计检测。它们相对于真随机数具有可快速产生、需要时可再生的优点，尤其对于程序调试。MonteCarlo程序中通常需要产生听从给定概率分布尸(X)的随机变量X。该步可用6,中的几种方法加以实现，其中包括苴雌和舍去滋；直接法(也称反演法或变换法)，需要转换与随机变量X相关的累积概率函数F(X)=PrOb(XQ(即：/(%)为Xx的概率)。0/(x)1明显表明，通过产生(0)内匀称分布随机数U，经转换我们可得听从尸(幻分布的随机样本X。为了得到这样的
7、具有概率分布户(%)的随机数X,不妨设U=JF(X),即可得X=K(U)(6)其中X具有分布函数/)。例如，若X是均值为呈指数分布的随机变量，且F(x)=l-,0x人时的/(X)=O,且F(X)上界dx为M(即：f(x)M),如图1所示。我们产生区间(0,1)内的两个随机数(U,U?),则X1=a+(b-a)Ui1=U2M(10)分别为在(a,b)和(0,M)内匀称分布的随机数。若1(X1)(11)则K为X的可选值，否则被舍去，然后再试新的一组(。1,口2)。如此运用舍去法，全部位于f(x)以上的点都被舍去，而位于F(X)上或以下的点都由X1=a+(b叫来产生XO图1 舍去法产生概率密度函数为
8、f(x)的随机变量例1设计一子程序使之产生0,1之间呈匀称分布的随机数Uo用该程序产生随机变。，其概率分布由下式给定T(9)=-(1-cos9),09J-SJ-X或者2(x)=x2-X2(16)方差的平方根称为斯准塞，即-(x)=(-x2)lz2(17)标准差给出了X在均值元四周的分布测度，并由此给出了误差幅度的阶数。戈的标准差与X的标准差的关系表示为。()=噌(18)该式表明，假如用依据(13)式由N个乙值构成的来估量工，那么结果中S在工四周的集中范围便与b(x)成比例，且随着样本数N的增加而降低。为了估量攵的集中范围，我们定义律初建为1NS2=y-x)2(19)由此式还可看出S?的期望值等
9、于2()。因此样本方差是(X)的无偏估量。将(19)式中平方项乘出来，便可得样本标准差为(20)当N较大时，系数N(N-1)可设为1。作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数B(M) =-PMqN-MM!(N M)!(21)该式表明N次独立随机试验中有“次胜利的概率。(21)中，P是一次试验中胜利的概率,且g=l-p0当M和N-M都较大时，便可用SHrling公式n,eny2n(22)因此(21)式可近似为正态分布17B(M)f(x)=1exp-,二?(23)(x)y2L2x)J其中元=NP且(x)=NNPq。因此当N时，中心极限定理表明，描述由N点MonteCarl
10、o算法获得的的分布的概率密度函数是(23)式所示的正态分布函数/(x)。也就是说,大量随机变量的集合趋于呈正态分布。将(18)式代入(23)式可得IV 1/V(X-X)2I/一、=lk-TTexP-G2/、(24)V 2(x)_2(X)正态(高斯)分布在工程，物理以及统计学的各类问题中都特别有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型常常用于其影响程度取决于由很多不规章的、浮动的元素构成的集合的状况。例2中我们给出了依据中心极限定理产生高斯随机变量的算法。由于样本数N是有限的，所以MomeCar1。算法不行能达到肯定的确定性。故而在工四周估量出某一范围或区间以确保我们估量的片落入该
11、范围内。假设要得到位于了-b和5+3之间的概率。由定义Pr0。元-xx+y=j*f(x)dx(25)令人=7(工二?、可得2V(x)Probx-xx+=-=，历北一万(26a)xx+z(26b)其中e炉(X)是误差函数，且/2是标准正态差的上二2分位点。对(26)式可做如下解释:假如用来呈现独立随机观测值并构建相关随机区间元5的MonteCarlo程序以较大的N值反复运行，则这些随机区间的分位点将近似等于戈。随机区间MS称为置信区间，Problx-M。(戈)v X 30时，一分布近似趋于正态分布。(26)式等价于nR-Sfa/2;AM_StaN-I.z-oProbX7=-xx+J=-=-a(2
12、8)W7V其中G/2：NT为自由度是N-I的同学氏/一分布的上2分位点；其值在任何标准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出上限建U(29)yN下限=5一当竺L(30)yjNMonteCarlo算法中关于误差估量的进一步争论，可参见例2采用中心极限定理产生一高斯(正态)分布的随机变量X。依据中心极限定理，大量均值四周的独立随机变量的总和,无论其个体变量的分布如何,总趋向于一种高斯分布。也就是说，对于任意随机数匕,i=l,2,.N,均值为方差为Var(丫)，NEYLNy(31)Z=I4NVariY慢慢与N合并为零均值、单位标准差的正态分布。若匕是匀称分布的随机变量(即Y1
13、=Ui),则P=I2,%Ny)=i,故而S-N/2Z=-,L(32)JN/12且变量X=Z+(33)近似等于均值为、方差为O?的正态变量。N小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设N=12,由于这样可消退(32)式中的平方根项。然而N的这种取值截掉了6b边界处的分布，且无法产生超过3。的值。对于分布曲线的两端比较重要的仿真，就必需用其它方案来产生高斯分布(参见2022)。这样，要产生一个均值为M、标准差为b的高斯随机变量X,就要遵循以下步骤：(1)生成12个匀称分布的随机数q,U2,u212(2)求得Z=Zq-6/=1(3)令X=Z+4数值积分对于一维积分，现已有一些求
14、积公式，如3.10节中所述。而对多重积分的公式则相对较少。MOnIeCar1。技术对此类多重积分的重要性至少存在两方面的缘由。多重积分的求积公式特别简单，而多重积分的MCM几乎保持不变。MonteCarlo积分的收敛性与维数无关,但对求积公式并非如此。人们已经发觉，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法23。这里将争论两类MonteCar1。积分方法，即简易MCM和具有对立变量的MCMo其它类型，如胜利一失败法和掌握变量法,可参见2426。本节还将简要介绍MCM在不法律规范积分中的应用。 4.1 简易MonteCarlo积分设要计算积分其中R
15、是维空间。令*=(*1*2,.*)是在/?内匀称分布的随机变量。则/(X)也是随机变量，其均值由下式给出27,28/(X)=/=(35)JIRrR方差为/2)=向JJ-扁川(36)其中IM=依(37)假如取X的N个独立样本，即X,X2,，Xn，它们与X具有相同的分布，且构成平均项/(Xj+f(X,+i+/(XM)=lg(j(38)我们便会期望该平均项接近于/(X)的均值。这样，由(35)和(38)式，即可得(39)Y(X,)MonteCarlo公式可用于有限区域R上的任何积分。为了举例说明，现将(39)式应用于一维和二维积分中。对于一维积分，设fb/=ff(x)dx(40)Ja由(39)式可得
16、(41)其中X,是区间(。/)内随机变量，即Xi =a + (b-a)U,0U 对于二维积分，有 = Cl(x2 WdX2 则相应的MonIeCarlo公式为/ =一一 W,x;) N /=I其中X =a + (h-a)Ul,0Ui 1Xr c + (d-c)U2U2 1(42)(43)(44)(45)(39)式中无偏估量值/收敛的很慢，这是由于估量值方差的级数是1/N。一种改良的方法，即掌握变量法，可通过减小估量值方差来提高其精确性和收敛性。 4.2 有对立变量的MonteCarlo积分方法术语就立交第29,30用来描述任一系列估量值，它们能相互抵消掉彼此的方差。便利起见，我们假设积分范围为
17、区间(0,1)。设要得到下面单重积分的估量值=fg(U)dU(46)我们期望表达式gg(U)+g(l-U)与g(U)相比具有更小的方差。假如g(U)太小，那反过来g(l-U)就很可能太大。因此，我们定义估量值二(ig(三)+g(l-Uj(47)NiZ其中Uj是0和1之间的随机数。此估量值方差的级数是,这是在(39)式基础上的一N4个极大的进步。对于两维积分I=s(fjU2)dUxdU2(48)则相应的估量值为1N1=g(5)+g(uu-U：)+g(lU；,，2)+g(l-U-哨(49)IV1=14依据相像性，可将该方法延拓至更高阶的积分。对于不是(U)的区间，可应用诸如(41)式到(45)式的
18、转换。例如，门(沙戈=3-g(uMu年(a)+g(iM)(5。)NJ=IZ其中g(U)=/(X),且X=a+0aQ。由(47)和(49)式可得，当积分维数增加时，每一维用来在简易MonteCarlo方法上提高效率的对立变量的最小值也会增加。这样使得简易MomeCar1。方法在多维运算中更具优越性。4.3不法律规范积分积分式/=g(x)clx(51)可用MonteCarlo仿真法进行估量31。对于概率密度为f(x)的随机变量X,其中/(x)在区间(0,8)上的积分等于1,则有瑞EZ(XMX(52)因此，为计算(51)式中的/值，首先要得到N个听从同一在区间(0,8)上的积分等于1的概率密度函数/
19、(X)的独立随机变量。其样本均值(53)便是/的估量值。例3用MonleCark)方法计算积分l=eiapdpd解：该积分式表示振幅和相位分别听从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是由于它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得MonteCarlo解的精确性。其闭合解为/().2)a其中，()是一阶BeSSel函数。图3给出由(44)和(45)式计算该积分的一个简洁的程序，其中=0,6=Lc=O以及d=2乃。该程序在VaX11/780中调用子程序RANDU来产生随机数Ul和U?。对于不同的N值，简易和对立变量MonteCar1。方法都可用于计算辐射积分。表1
20、中对a=5时的结果进行了精确的比较。在应用(49)式时，用到了下面的对应式：UiXU2X2,-Uyb-Xi=(b-a-Ux-U2d-X2=(d-c-U2)表1例3辐射积分的MC方法积分结果N简易MCM对立变量MCM5001000-0.5737+j0.0808200040006000800010,000精确解：-0.4116+j05位势问题的解位势理论与布朗运动(或随机游动)的关系是由KakUtani在1944年首次提出的32。自此，所谓的概率位势理论便在诸如热传导3引一38、静电学3946以及电力工程等很多学科领域得到广泛应用。不同方程式的概率解或MC解的一个基本概念就是随机游动。不同类型的随
21、机游动应用不同的MonteCarlo方法。最常见的类型是固定随机游动种自动定OOOl000200030004000500060007000800090010001100120013001400150016001700180019002000210022002300240025002600270028002900300031003200330034003500360037003800390040位随机游动其它不常见类型有迂离法、缩减边界法、内切形法以及表面密度法。C*CINTEGRATIONUSINGCRUDEMONTECARLOCANDANTITHETICMETHODSCONLYFEWLINE
22、SNEEDBECHANGEDTOUSETHISPROGRAMCFORANYMULTI-DIMENSIONALINTEGRATIONC*DATAISIJS2,IS3JS41234,5678,9012,3456/DATAA,B,C0.0,i.0,0.0CSPECIFYTHEINTEGRANDF(RHO,PHI)=RHO*CEXP(J35cALPHA*RHO*COS(PHI)J=(O,OJ-O)ALPHA=5.0PIE=3.I415927D=2.0*PIEDO30NRUN=500,10000,500SUMl=(OAO-O)SUM2=(0.0,0.0)DO10I=l,NRUNCALLRANDU(ISI
23、JS2,U1)CALLRANDU(IS3JS4,U2)X1=A(B-A)*UIX2=C+(D-C)*U2X3=(B-A)*(1.0-Ul)X4=(D-C)*(1.0-U2)SUMl=SUMl+F(X1,X2)SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4)10CONTINUEAREA1=(B-A)*(D-C)*SUM1FLOAT(NRUN)AREA2=(B-A)*(D-C)*SUM2(4.0*FLOAT(NRUN)PRINT*,NRUN,AREA1,AREA2WRITE(6,*)NRUN,AREA1,AREA2WRITE(6,20)NRUN,AREAl
24、,AREA220FoRMAT(2X,NRUN=I5,3X,AREA1=F12.6,3X,Fl2.6,AREA2=,1F12.6,3X,F12.6)30CONTINUESTOPEND图3例3中用MonteCarlo方法计算二维积分的程序5.1 固定随机游动为详细起见，假设用固定随机游动的MCM解拉普拉斯方程V2V=O(在区域R内)(54a)满意DiriChIet边界条件V=Vp(在边界B上)(54b)首先将R分成网状结构，并用其有限差当量替代(54a)在二维R内的有限差表达式可由(3.26)式给出，即V(x,y)=P,v+V(x+,y)+Pa-V(x-,y)+pv+V(x,y+)+pv_V(x,
25、y-)(55a)其中px+=Px-=Py+=Py-=-(55b)(55)式中，假设网络的一个方格尺寸是,如图4所示。下面给出该方程的概率解释。若随机游动粒子在某一瞬时处于(x,y)位置处，它从今点向(x+Ay),(x-A,y)(x,y+A)及(x,y-A)移动的概率分别是P,P,Ph和py_。打算粒子移动方式的方法是产生一个随机数U,0Uvl,并按下面的方式指导粒子的随机游动:0U0.250.25 U 0.50.5 (/ 0.750.75U )的MonteCarlo解为1N2ng-i.v(x)=-vp(/)+-g(j,yj)/V/=I-/Vi=1LJT对刚才所描述MCM的一个好玩的分析就是游走醉鬼问题115,35。我们将随机游动粒子当作一个“醉鬼”，将网格方块当作“城市街区”，结点当作“十字路口。边界8当作“城市边界”，而且把边界B上的结点当作“警察”。尽管醉鬼想步行走回家，但他却醉的一塌糊涂以至于在整个城市中随机游走。警察的工作是在城市边界上一观察醉鬼便将其抓获，并勒令醉鬼交付罚金Vp。那么醉鬼要支付的预期罚金将会是多少呢？其答案便是(57)式。在电介质边界上，指定边界条件为=。2”。
