蒙特卡洛算法详讲

1、Monte Carlo 法概述Monte Carlo法不同于前面几章所介绍的确定性数值方法，它是用来解决数 学和物理问题的非确定性的（概率统计的或随机的）数值方法。Monte Carlo方法（MCM）,也称为统计试验方法，是理论物理学两大主要学科的合并：即随机 过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是 研究均匀介质的稳定状态1。它是用一系列随机数来近似解决问题的一种方法， 是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解 的处理数学问题的一种手段。运用该近似方法所获得的问题的解in spirit更接近于物理实验结果，而不是经典数值计算结果。普遍

2、认为我们当前所应用的 MC技术，其发展约可追溯至1944年，尽管在 早些时候仍有许多未解决的实例。MCM的发展归功于核武器早期工作期间 Los Alamos （美国国家实验室中子散射研究中心）的一批科学家。Los Alamos小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了 MCM在各种问题中的 应用2-4。Monte Carlo”的名称取自于Monaco （摩纳哥）内以赌博娱乐而闻名 的一座城市。Monte Carlo方法的应用有两种途径：仿真和取样。仿真是指提供实际随机 现象的数学上的模仿的方法。一个典型的例子就是对中子进入反应堆屏障的运动 进行仿真，用随机游动来模仿中子的锯齿形路径

3、。 取样是指通过研究少量的随机 的子集来演绎大量元素的特性的方法。例如，f（x）在a x b上的平均值可以通过间歇性随机选取的有限个数的点的平均值来进行估计。这就是数值积分的 Monte Carlo方法。MCM已被成功地用于求解微分方程和积分方程，求解本征 值，矩阵转置，以及尤其用于计算多重积分。任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数 的方法。这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以 及其它竞赛活动中都会出现。Monte Carlo计算方法需要有可得的、服从特定概 率分布的、随机选取的数值序列。机数和随机变量的产生 10全面的论述了产生随机数的各

4、类方法。其中较为普遍应用的产生随 机数的方法是选取一个函数g（x）,使其将整数变换为随机数。以某种方法选取孔，并按照xk1 g（xj产生下一个随机数。最一般的方程g（x）具有如下形式:g(x) (ax c) mod m其中x0初始值或种子（x0 0）a 乘法器（a 0）c 增值（c 0 ）m模数对于t数位的二进制整数，其模数通常为 2;例如，对于31位的计算机m即可取231 1。这里Xo,a和c都是整数，且具有相同的取值范围m a,m c, m X0。所需的随机数序Xn便可由下式得xn 1 (axn c) mod m该序列称为线性同余序列。例如，若xo a c 7且m 10,则该序列为7,6,

5、9,0,7,6,9,0(8可以证明，同余序列总会进入一个循环套； 也就是说，最终总会出现一个无休止 重复的数字的循环。(8.3)式中序列周期长度为4。当然,一个有用的序列必是 具有相对较长周期的序列。许多作者都用术语 乘同余法和混合同余法分别指代c 0和c 0时的线性同余法。选取x0,a,c和m的法则可参见6,10。这里我们只关心在区间(0,1)内服从均匀分布的随机数的产生。用字符 U来表示这些数字，则由式(8.2)可得这样U仅在数组0,1/m,2/m,(m 1)/m中取值。(对于区间(0,1)内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。其它检测方法可参见 3,6 0)产生区

6、间(a,b)内均匀分布的随机数X a (b a)U用计算机编码产生的随机数(利用式(8.2)和(8.4)并不是完全随机的； 事实上，给定序列种子，序列的所有数字U都是完全可预测的。一些作者为强调 这一点，将这种计算机产生的序列称为 伪随机数。但如果适当选取a,c和m,序列U的随机性便足以通过一系列的统计检测。它们相对于真随机数具有可快速产 生、需要时可再生的优点，尤其对于程序调试。Monte Carlo程序中通常需要产生服从给定概率分布F(x)的随机变量X。该步可用6, 1315中的几种方法加以实现，其中包括 直接法和舍去法。直接法(也称反演法或变换法)，需要转换与随机变量X相关的累积概率函数

7、F(x) prob(X x)(即：F(x)为X x的概率)。0 F(x) 1显然表明，通过产生(0,1)内均匀分布随机数U，经转换我们可得服从F(x)分布的随机样本X 为了得到这样的具有概率分布 F(x)的随机数X ,不妨设U F(x),即可得 _ 1 X F (U) 其中X具有分布函数F(x)0例如，若X是均值为呈指数分布的随机变量，且x/ F(x) 1 e ,0 x (8.7)在U F(x)中解出X可得X ln(1 U ) 由于(1 U)本身就是区间(0,1)内的随机数，故可简写为ln U(8有时(8.6)式所需的反函数F 1(x)不存在或很难获得。这种情况可用舍去法来 处理。令f(x)

8、尤凶为随机变量X的概率密度函数。令a x b时的f(x) 0, dx且f(x)上界为M (即：f(x) M ),如图8.1所示。我们产生区间(0,1)内的两个 随机数(U1,U2),则a (b a)U1f1 U 2M 分别为在(a,b)和(0,M)内均匀分布的随机数。若f1f(XJ则Xi为X的可选值，否则被舍去，然后再试新的一组(Ui,U2)。如此运用舍去法， 所有位于f(x)以上的点都被舍去，而位于f (x)上或以下的点都由X1 a (b a)U1 来产生 X1。f (x j A图8.1舍去法产生概率密度函数为f(x)的随机变量例8.1设计一子程序使之产生0,1之间呈均匀分布的随机数U 。用

9、该程序产生随机变，其概率分布由下式给定1T()万。cos ),0解：生成U的子程序如图8.2所示。该子程序中，m 221 1 2147483647,c 0,且a 75 16807。应用种子数(如1234),主程序中每调用一次子程序，就会生 成一个随机数U。种子数可取1到m间的任一整数。0001C*0002 C PROGRAM FOR GENERATING RANDOM VARIABLES WITH0003 C A GIVEN PROBABILITY DISTRIBUTION0004000500060005000600070008000900100011001200131000140015DOU

10、BLE PRECISION ISEEDISEED=1234.D0DO 10 I=1,100CALL RANSOM(ISEED,R)THETA=ACOSD(1.0-2.0*R)WRITE(6,*)I,THETACONTINUESTOPEND0001C*0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERSIN0003 C THE INTERVAL (0,1) 0004C* 0005 0006SUBROUTINE RANDOM (ISEED,R)0007DOUBLE PRECISION ISDDE,DEL,A0008DATA DEL,A/2147483647

11、.D0,16807.D0/0009 0010ISDDE=DMOD(A*ISDDE,DEL)0011R=ISDDE/DEL0012RETURN0013END图8.2 例8.1的随机数生成器 图8.2的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有 生成随机数的子程序。 为了生成随机变量，令 1 U T( ) -(1 cos )则有T 1(U) cos 1 (1 2U )据此，一系列具有给定分布的随机变量便可由图8.2所示主程序中生成。8.3 误差计算Monte Carlo程序给出的解按大量的检测统计都达到了平均值。因此，该解 中包含了平均值附近的浮动量，而且不可能达到100%的置信度。

12、要计算 MonteCarlo算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简 要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计13,16。设X是随机变量。则X的期望值或均值X定义为 x xf(x)dx这里f(x)是X的概率密度分布函数。如果从 f(x)中取些独立的随机样本Xi,X2，,Xn ,那么的x估计值就表现为N个样本值的均值。N n1Xnx是x的真正的平均值，而？只是x的有着准确期望值的无偏估计。虽然 ？ 的期望值等于x,但？ x。因此，我们还需要？的值在x附近的分布测度。为了估计X以及?在x附近的的值的分布，我们需要引入X的方差，其定义 为X与x差的平方的期望

13、值，即Var(x) 2 (x x)2 (x x)2 f (x)dx由(x x)2 x2 2xx x2 ,故有2 (x) x2 f (x)dx 2x xf (x)dx x2 f (x)dx或者2(x) x2 x2方差的平方根称为标准差，即(x) (x2 x2)1/2标准差给出了 x在均值x附近的分布测度，并由此给出了误差幅度的阶数。？的标准差与x的标准差的关系表示为(?)(x)(?)N该式表明，如果用根据(8.13)式由N个xn值构成的？来估计x,那么结果中0在x 附近的扩散范围便与 (x)成比例，且随着样本数N的增加而降低。为了估计？的扩散范围，我们定义样本方差为c1 NCS (xn 2N 1

14、 n 1(8.19)由此式还可看出(8.19)由此式还可看出S2的期望值等于2(x)。因此样本方差是2(x)的无偏估计。将平方项乘出来，便可得样本标准差为1/21 /21/21 N1xn父N n 1当N较大时，系数N /(N 1)可设为作为获得中心极限定理的一种方法,概率论的一个基本解可考虑二次项函数作为获得中心极限定理的一种方法,概率论的一个基本解可考虑二次项函数N!B(M N!B(M )M !(N M )!该式表明N次独立随机试验中有M次成功的概率。(8.21)中，p是一次试验中成功的概率，且q 1 p。当M和N M都较大时，便可用Stirling公式因此(8.21)式可近似为正态分布17

15、1( x)2B(M) f(X) exp 2,(X),22 (x)其中X Np且(X)丽。因此当N 时，中心极限定理表明，描述由N点 Monte Carlo算法获得的？的分布的概率密度函数是(8.23)式所示的正态分布函 数f (x)。也就是说，大量随机变量的集合趋于呈正态分布。将(8.18)式代入(8.23)式可得f(?楞f(?楞4rxpN(? x)22 2(x)正态(高斯)分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型 的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许 多不规则的、浮动的元素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变

16、量的算法。由于样本数N是有限的，所以Monte Carlo算法不可能达到绝对的确定性。 故而在x附近估计出某一范围或区间以确保我们估计的 ？落入该范围内。假设要 得到？位于x和x之间的概率。由定义xProb x 及 xx f(?)d?令，士 x 可得2/N xProb x 父 Prob x 父 x2 N /20erf N /2x(8.26句或Prob x z /2 )? x z /2 1,N. N(8.26b)其中erf (x)是误差函数，且z /2是标准正态差的上/2分位点。对(8.26)式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间x 的MonteCarlo程序以较大的N值反

17、复运行，则这些随机区间的erf J分位点将2 x近似等于农。随机区间x称为置信区间，erf jg称为置信度。大多数的Monte Carlo的Monte Carlo算法都使用误差x 行,这表明？是在实际均值x的标准0.6826,0.954,0.997,为已知。由于这种情况很少t-分布取代正态分布。我们知道当 N很大，比如N 30时，t分布近似趋于正态分布。(8.26)式等价于Prob xSt /2;N 1N其中t 0.6826,0.954,0.997,为已知。由于这种情况很少t-分布取代正态分布。我们知道当 N很大，比如N 30时，t分布近似趋于正态分布。(8.26)式等价于Prob xSt /

18、2;N 1N其中t /2;N1为自由度是N 1的学生氏t分布的上/ 2分位点；其值在任何标准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出St /2;N 1St /2;N 1.NMonte Carlo算法中关于误差估计的进步讨论，可参见18,19。例8.2利用中心极限定理产生一高斯(正态)分布的随机变量X 。根据中心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何,差范围内的。由(8.26)式可得，样本均值？位于区间？x/JN内的概率是0.6826或68.3%。若要求更高的置信度，可用两个或三个标准差的值Prob x M (x) x x M (x)N. N(

19、8.27)其中M是标准差的个数。在(8.26)和(8.27)式中均假设总体标准差 出现，故必须用由(8.20)式算得的样本标准差S来估计 的值，从而使学生氏总趋向于一种高斯分布。也就是说，对于任意随机数 丫1,2,.N ,均值为Y ,方差为Var(Y),NY NYZ i 1.NVar Y渐渐与N合并为零均值、单位标准差的正态分布。若 Yi是均匀分布的随机变量(即 Y UJ,则 Y 1/2,Var(Y) 1/02 ,故而NUi N /2且变量X Z近似等于均值为、方差为2的正态变量。N小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设 N 12,因为这样可消除(8.32)式中 的平

20、方根项。然而N的这种取值截掉了 6边界处的分布，且无法产生超过3 的值。对于分布曲线的两端比较重要的仿真， 就必须用其它方案来产生高斯分布(参见20 22)。这样，要产生一个均值为、标准差为的高斯随机变量X ,就要遵循以下步骤：(1)生成12个均匀分布的随机数Ui,U2,.,Ui212(2)求得 Z Ui 6 i 1(3)令 X Z 8.4数值积分对于一维积分，现已有一些求积公式，如3.10节中所述。而对多重积分的公式则相对较少。Monte Carlo技术对此类多重积分的重要性至少存在两方面的 原因。多重积分的求积公式非常复杂，而多重积分的MCM几乎保持不变。Monte Carlo积分的收敛性

21、与维数无关，但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维 积分的很有效的方法23。这里将讨论两类Monte Carlo积分方法，即简易MCM 和具有对立变量的MCM。其它类型，如成功失败法和控制变量法，可参见24 26。本节还将简要介绍MCM在不规范积分中的应用。 8.4.1 简易 Monte Carlo 积分 设要计算积分I fR其中R是n维空间。令X X1,X2,.Xn是在R内均匀分布的随机变量。则f(X)也是随机变量，具均值由下式给出27,28-f(X)也是随机变量，具均值由下式给出27,28-1f(X) RRfR RdX,

22、它们与R RdX,它们与X具有相同的分布，且构f Xi f X2Nf Xn1 Nf XiN i 1(8.35)方差力1Var f X 一R R其中如果取X的N个独立样本，即Xi,X2,.,Xn 成平均项我们便会期望该平均项接近于f(X)的均值。这样，由(8.35)和(8.38)式，即可得_R N_R NN i 1f XiMonte Carlo公式可用于有限区域R上的任何积分。为了举例说明，现将(8.39) 式应用于一维和二维积分中。对于一维积分，设bI f(x)dxabI f(x)dxa(8.40)由(8.39)式可得f Xi其中Xi是区间a,b内随机变量，即Xi a b a U ,0 U 1

23、对于二维积分，有I b d f X1,XgU dX1dX 2a cgU则相应的Monte Carlo公式为(b a) d cN(b a) d cN12f Xi,Xi其中X； a (b a)U1,0 U1 1 222Xi2 c (d c)U 2,0 U 2 1无偏估计值I收敛的很慢，这是由于估计值方差的级数是1/N 一种改良的方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛 性。 8.4.2具有对立变量的Monte Carlo积分方法术语对立变量29,30用来描述任一系列估计值，它们能互相抵消掉彼此的方 差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值1I

24、 0 g(U)dU1.一 .我们期望表达式1 g(U) g(1 U)与g(U)相比具有更小的方差。如果g(U)太小，那反过来g(1 U)就很可能太大。因此，我们定义估计值， 一 1、一,，、其中Ui是0和1之间的随机数。此估计值方差的级数是 二，这是在(8.39)式N4基础上的一个极大的进步。对于两维积分1 1I gU1,U2 dU1dU200、则相应的估计值为N 112121212IgUi,Ui gUi,1 Ui g1 3。 g 1 5,1 5N i 1 4根据相似性，可将该方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如(8.41)式到(8.45)式的转换。例如，b1f x d

25、x b a g U dUa0b a N 1g Ui g 1 UiN i 1 2其中g(U) f(X),且X a b aU。由(8.47)和(8.49)式可得，当积分维数增加时，每一维用来在简易Monte Carlo方法上提高效率的对立变量的最小 值也会增加。这样使得简易 Monte Carlo方法在多维运算中更具优越性。不规范积分积分式I 0 g(x)dx可用Monte Carlo仿真法进行估计31。对于概率密度为f(x)的随机变量X ,其中f(x)在区间(0,)上的积分等于1,则有g x-g x-dx0 f x0 g(x)dx因此，为计算(8.51)式中的I值，首先要得到N个服从同一在区间(

26、0,)上的 积分等于1的概率密度函数f(x)的独立随机变量。其样本均值1 N g x1 N g xiN1 f xi便是I的估计值 例8.3 用Monte Carlo方法计算积分2.I ej cos d d0 0解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得， 由此便可得Monte Carlo解的精确性。其闭合解为Ji其中J1是一阶Bessel函数。图8.3给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其中a 0,b 1,c 0以及d 2。该程序在Vax 11/780中调用子程序RAND

27、U来产生随机数UDU2。对于不同的N值，简易和对立变量Monte Carlo方法都可用于计算辐射积分。表8.1中对 5时的结果进行了精确的比较。在应用（8.49） 式时，用到了下面的对应式：U1 X1,U2 X2,1 U1 b X1 b a 1 U12221 U d X d c 1 U表8.1例8.3辐射积分的MC方法积分结果N简易MCM 对立变量MCM500-0.2892-j0.0-0.2887-j0.0585742-0.4982-j0.0080-0.4682-j0.0082-0.4216-j0.0323-0.4982-j0.0080-0.4682-j0.0082-0.4216-j0.032

28、3-0.3787-j0.0440-0.4139-j0.0241-0.4121-j0.02402000-0.4922-j0.00404000-0.3999-j0.03456000-0.3608-j0.02708000-0.4327-j0.037810,00-0.4229-j0.00237精确解：-0.4116+j0 8.5位势问题的解位势理论与布朗运动（或随机游动）的关系是由 Kakutani在1944年首次提 出的32。自此，所谓的概率位势理论便在诸如热传导33 38、静电学39 46以及电力工程等许多学科领域得到广泛应用。不同方程式的概率解或MC解的一个基本概念就是随机游动。不同类型的随机游

29、动应用不同的Monte Carlo方法。最常见的类型是固定随机游动和自动定位随机游动。其它不常见类型有迁离 法、缩减边界法、内切形法以及表面密度法。0001C*0002 C INTEGRATION USING CRUDE MONTE CARLO0003 C AND ANTITHETIC METHODS0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM0005 C FOR ANY MULTI DIMENSIONAL INTEGRATION 0006C*00070008DATAIS1,IS2,IS3,IS4/1234,5678,901

30、2,3456/0009DATAA,B,C/0.0,1.0,0.0/00100011 C SPECIFY THE INTEGRAND0012F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI)00130014J=(0.0,1.0)0015ALPHA=5.00016PIE=3=2.0*PIE0018DO 30 NRUN=500,10000,5000019SUM1=(0.0,0.0)0020SUM2=(0.0,0.0)0021DO 10 I=1,NRUN0022CALL RANDU(IS1,IS2,U1)0023CALL RANDU(IS3,IS

31、4,U2)0024X1=A+(B-A)*U10025X2=C+(D-C)*U20026X3=(B-A)*(1.0-U1)0027X4=(D-C)*(1.0-U2)0028SUM1=SUM1+F(X1,X2)0029SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4)003010 CONTINUE0031AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN)0032003300320033003400350036PRINT *,NRUN,AREA1,AREA2WRITE(6,*) NRUN,AREA1,AREA2WRITE(6,20) NRUN

32、,AREA1,AREA220003700380039004030200037003800390040301 CONTINUE STOP ENDF12.6,3X,F12.6,/) 8.5.1图8.3 例8.3中用Monte Carlo方法计算二维积分的程序 固定随机游动 8.5.1为具体起见，假设用固定随机游动的 MCM解拉普拉斯方程2V 0(在区域 R内)(8.54句满足Dirichlet边界条件V VP(在边界 B 上)(8.54b)首先将R分成网状结构，并用具有限差当量替代2。(8.54a)在二维R内的有限差表达式可由(3.26)式给出，即V x, yPxV x ,y Px V x ,y

33、Py V x, yPy V x, y(8.55句其中1PxPxPyPy4 (8.55b)，假设网络的一个方格尺寸是，如图8.4所示。下面给出该方程的概率解释。若随机游动粒子在某一瞬时处于x,y的概率解释。若随机游动粒子在某一瞬时处于x,y位置处，它从此点向(x ,y),(x,y)(x, y )及x, y移动的概率分别是Px , Px , Py和Py。决定粒子移动方式的方法是产生一个随机数的随机游动：定粒子移动方式的方法是产生一个随机数的随机游动：x,yx, yx,yx, yx,yx,yx,yx,yU,0 u 1,并按下面的方式指导粒子0 U 0.250.25 U0.50.5 U 0.750.75 U1如果不用方格而用矩形格，则有PxP如果不用方格而用矩形格，则有PxPx 且 PyP