蒙特卡洛算法详讲

1、Monte Carlo 法§8.1 概述Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。Monte Carlo 方法（MCM），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态1.R.Hersch and R.J.Griego,“Brownian motion and potential theory,”Sci.Amer.,Mar.1969,pp.67-74.。它是用一系列随机数来近似解决问题的一种

2、方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解in spirit更接近于物理实验结果，而不是经典数值计算结果。普遍认为我们当前所应用的MC技术，其发展约可追溯至1944年，尽管在早些时候仍有许多未解决的实例。MCM的发展归功于核武器早期工作期间Los Alamos（美国国家实验室中子散射研究中心）的一批科学家。Los Alamos 小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MCM在各种问题中的应用。“Monte Carlo”的名称取自于Monaco（摩纳哥）内以赌博娱乐而闻名的一座城市。Monte

3、Carlo 方法的应用有两种途径：仿真和取样。仿真是指提供实际随机现象的数学上的模仿的方法。一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动来模仿中子的锯齿形路径。取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。例如，在上的平均值可以通过间歇性随机选取的有限个数的点的平均值来进行估计。这就是数值积分的Monte Carlo 方法。MCM已被成功地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以及其它竞赛活动

4、中都会出现。Monte Carlo 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。§8.2 随机数和随机变量的产生510全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数，使其将整数变换为随机数。以某种方法选取，并按照产生下一个随机数。最一般的方程具有如下形式： （8.1）其中 初始值或种子（） 乘法器（）增值（） 模数对于数位的二进制整数，其模数通常为。例如，对于31位的计算机即可取。这里和都是整数，且具有相同的取值范围。所需的随机数序便可由下式得 （8.2）该序列称为线性同余序列。例如，若且，则该序列为 7，6，9，0，7，6，9，

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

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

7、机数，不妨设，即可得 （8.6）其中具有分布函数。例如，若是均值为呈指数分布的随机变量，且 （8.7）在中解出可得 （8.8）由于本身就是区间(0,1)内的随机数，故可简写为 （8.9）有时(8.6)式所需的反函数不存在或很难获得。这种情况可用舍去法来处理。令为随机变量的概率密度函数。令时的，且上界为（即：），如图8.1所示。我们产生区间(0,1)内的两个随机数，则 （8.10）分别为在(a,b)和(0,M)内均匀分布的随机数。若 （8.11）则为的可选值，否则被舍去，然后再试新的一组。如此运用舍去法，所有位于以上的点都被舍去，而位于上或以下的点都由来产生。 图8.1 舍去法产生概率密度函数为

8、的随机变量 例8.1 设计一子程序使之产生0,1之间呈均匀分布的随机数。用该程序产生随机变，其概率分布由下式给定 解：生成的子程序如图8.2所示。该子程序中，且。应用种子数（如1234），主程序中每调用一次子程序，就会生成一个随机数。种子数可取1到间的任一整数。 0001 C*0002 C PROGRAM FOR GENERATING RANDOM VARIABLES WITH0003 C A GIVEN PROBABILITY DISTRIBUTION0004 C* 0005 0006 DOUBLE PRECISION ISEED 0007 0008 ISEED=1234.D0 0009 D

9、O 10 I=1,100 0010 CALL RANSOM(ISEED,R) 0011 THETA=ACOSD(1.0-2.0*R) 0012 WRITE(6,*)I,THETA 0013 10 CONTINUE 0014 STOP 0015 END 0001 C* 0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN 0003 C THE INTERVAL (0,1) 0004 C* 0005 0006 SUBROUTINE RANDOM (ISEED,R) 0007 DOUBLE PRECISION ISDDE,DEL,A 0008 DATA

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

11、e Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计13,16。设是随机变量。则的期望值或均值定义为 （8.12）这里是的概率密度分布函数。如果从中取些独立的随机样本，那么的估计值就表现为个样本值的均值。 （8.13）是的真正的平均值，而只是的有着准确期望值的无偏估计。虽然的期望值等于，但。因此，我们还需要的值在附近的分布测度。为了估计以及在附近的的值的分布，我们需要引入的方差，其定义为与差的平方的期望值，即 （8.14）由，故有 （8.15）或者 （8.16）方差的平方根称为标准差，即 （8.17）标准差给

12、出了在均值附近的分布测度，并由此给出了误差幅度的阶数。的标准差与的标准差的关系表示为 （8.18）该式表明，如果用根据(8.13)式由个值构成的来估计，那么结果中在附近的扩散范围便与成比例，且随着样本数的增加而降低。为了估计的扩散范围，我们定义样本方差为 （8.19）由此式还可看出的期望值等于。因此样本方差是的无偏估计。将（8.19）式中平方项乘出来，便可得样本标准差为 （8.20）当较大时，系数可设为1。 作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数 （8.21）该式表明次独立随机试验中有次成功的概率。(8.21)中，是一次试验中成功的概率，且。当和都较大时，便可用St

13、irling 公式 （8.22）因此(8.21)式可近似为正态分布17 (8.23)其中且。因此当时，中心极限定理表明，描述由点Monte Carlo算法获得的的分布的概率密度函数是（8.23）式所示的正态分布函数。也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得 （8.24）正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。由于样本数是有限的，所以Monte Carlo

14、算法不可能达到绝对的确定性。故而在附近估计出某一范围或区间以确保我们估计的落入该范围内。假设要得到位于和之间的概率。由定义 （8.25）令可得 （8.26a）或 （8.26b）其中是误差函数，且是标准正态差的上分位点。对（8.26）式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间的Monte Carlo 程序以较大的值反复运行，则这些随机区间的分位点将近似等于。随机区间称为置信区间，称为置信度。大多数的Monte Carlo 算法都使用误差，这表明是在实际均值的标准差范围内的。由（8.26）式可得，样本均值位于区间内的概率是0.6826或68.3。若要求更高的置信度，可用两个或三个

15、标准差的值。如 （8.27）其中是标准差的个数。 在（8.26）和（8.27）式中均假设总体标准差为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差来估计的值，从而使学生氏分布取代正态分布。我们知道当很大，比如时，分布近似趋于正态分布。（8.26）式等价于 （8.28）其中为自由度是的学生氏分布的上分位点；其值在任何标准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出 上限 （8.29） 下限 （8.30） Monte Carlo 算法中关于误差估计的进一步讨论，可参见18,19。例8.2 利用中心极限定理产生一高斯（正态）分布的随机变量。根据中心极限定理，大

16、量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数，均值为，方差为， （8.31）渐渐与合并为零均值、单位标准差的正态分布。若是均匀分布的随机变量（即），则，故而 （8.32）且变量 （8.33）近似等于均值为、方差为的正态变量。小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设，因为这样可消除（8.32）式中的平方根项。然而的这种取值截掉了边界处的分布，且无法产生超过的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布（参见2022）。 这样，要产生一个均值为、标准差为的高斯随机变量，就要遵循以下步骤： （

17、1） 生成12个均匀分布的随机数 （2） 求得 （3） 令§8.4 数值积分 对于一维积分，现已有一些求积公式，如3.10节中所述。而对多重积分的公式则相对较少。Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常复杂，而多重积分的MCM几乎保持不变。Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法23。这里将讨论两类Monte Carlo 积分方法，即简易MCM 和具有对立变量的MCM。其它类型，如成功失败法和控制变

18、量法，可参见2426。本节还将简要介绍MCM在不规范积分中的应用。§8.4.1 简易Monte Carlo 积分设要计算积分 （8.34）其中是维空间。令是在内均匀分布的随机变量。则也是随机变量，其均值由下式给出27,28 (8.35)方差为 （8.36）其中 （8.37）如果取的个独立样本，即，它们与具有相同的分布，且构成平均项 （8.38）我们便会期望该平均项接近于的均值。这样，由（8.35）和（8.38）式，即可得 （8.39）Monte Carlo 公式可用于有限区域上的任何积分。为了举例说明，现将（8.39）式应用于一维和二维积分中。 对于一维积分，设 （8.40）由（8.

19、39）式可得 （8.41）其中是区间内随机变量，即 （8.42） 对于二维积分，有 （8.43）则相应的Monte Carlo 公式为 （8.44）其中 （8.45） （8.39）式中无偏估计值收敛的很慢，这是由于估计值方差的级数是。一种改良的方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛性。§8.4.2 具有对立变量的Monte Carlo 积分方法 术语对立变量29,30用来描述任一系列估计值，它们能互相抵消掉彼此的方差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值 （8.46）我们期望表达式与相比具有更小的方差。如果太小，那反过来就很可

20、能太大。因此，我们定义估计值 （8.47）其中是0和1之间的随机数。此估计值方差的级数是，这是在（8.39）式基础上的一个极大的进步。对于两维积分 （8.48）则相应的估计值为 （8.49）根据相似性，可将该方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如（8.41）式到（8.45）式的转换。例如， （8.50）其中，且。由（8.47）和（8.49）式可得，当积分维数增加时，每一维用来在简易Monte Carlo 方法上提高效率的对立变量的最小值也会增加。这样使得简易Monte Carlo 方法在多维运算中更具优越性。§8.4.3 不规范积分 积分式 （8.51）可用Mo

21、nte Carlo 仿真法进行估计31。对于概率密度为的随机变量，其中在区间(0,)上的积分等于1，则有 （8.52）因此，为计算（8.51）式中的值，首先要得到个服从同一在区间(0,)上的积分等于1的概率密度函数的独立随机变量。其样本均值 （8.53）便是的估计值。例8.3 用Monte Carlo 方法计算积分 解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得Monte Carlo 解的精确性。其闭合解为 其中是一阶Bessel 函数。 图8.3给出由（8.44）和（8.45）式计算该积

22、分的一个简单的程序，其中以及。该程序在Vax 11/780 中调用子程序RANDU来产生随机数和。对于不同的值，简易和对立变量Monte Carlo 方法都可用于计算辐射积分。表8.1中对时的结果进行了精确的比较。在应用（8.49）式时，用到了下面的对应式： 表8.1 例8.3辐射积分的MC方法积分结果N简易MCM对立变量MCM500-0.2892-j0.0742-0.2887-j0.05851000 -0.5737+j0.0808 -0.4982-j0.00802000 -0.4922-j0.0040 -0.4682-j0.00824000 -0.3999-j0.0345-0.4216-j0

23、.03236000-0.3608-j0.0270 -0.3787-j0.04408000-0.4327-j0.0378 -0.4139-j0.024110,000 -0.4229-j0.0237 -0.4121-j0.0240精确解：-0.4116+j0§8.5 位势问题的解位势理论与布朗运动（或随机游动）的关系是由Kakutani在1944年首次提出的32。自此，所谓的概率位势理论便在诸如热传导3338、静电学3946以及电力工程等许多学科领域得到广泛应用。不同方程式的概率解或MC解的一个基本概念就是随机游动。不同类型的随机游动应用不同的Monte Carlo方法。最常见的类型是固

24、定随机游动和自动定位随机游动。其它不常见类型有迁离法、缩减边界法、内切形法以及表面密度法。 0001 C* 0002 C INTEGRATION USING CRUDE MONTE CARLO 0003 C AND ANTITHETIC METHODS 0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM 0005 C FOR ANY MULTIDIMENSIONAL INTEGRATION 0006 C* 0007 0008 DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/ 0009 DATA A

25、,B,C/0.0,1.0,0.0/ 0010 0011 C SPECIFY THE INTEGRAND 0012 F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI) 0013 0014 J=(0.0,1.0) 0015 ALPHA=5.0 0016 PIE=3.1415927 0017 D=2.0*PIE 0018 DO 30 NRUN=500,10000,500 0019 SUM1=(0.0,0.0) 0020 SUM2=(0.0,0.0) 0021 DO 10 I=1,NRUN 0022 CALL RANDU(IS1,IS2,U1) 0023 CALL RAN

26、DU(IS3,IS4,U2) 0024 X1=A+(B-A)*U1 0025 X2=C+(D-C)*U2 0026 X3=(B-A)*(1.0-U1) 0027 X4=(D-C)*(1.0-U2) 0028 SUM1=SUM1+F(X1,X2) 0029 SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4) 0030 10 CONTINUE 0031 AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN) 0032 AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN) 0033 PRINT *,NRUN,AREA1,

27、AREA2 0034 WRITE(6,*) NRUN,AREA1,AREA2 0035 WRITE(6,20) NRUN,AREA1,AREA2 0036 20 FORMAT(2X,NRUN=,I5,3X,AREA1=,F12.6,3X,F12.6,AREA2=, 0037 1 F12.6,3X,F12.6,/) 0038 30 CONTINUE 0039 STOP 0040 END 图8.3 例8.3中用Monte Carlo方法计算二维积分的程序§8.5.1 固定随机游动 为具体起见，假设用固定随机游动的MCM解拉普拉斯方程 （在区域内） （8.54a）满足Dirichlet边界条件 （在边界上） （8.54b）首先将分成网状结构，并用其有限差当量替