1、第八章Monte Carlo 法8.1 概述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、随机游动来模仿中子的锯齿形路径。取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。例如， 在 上的平均值)(xfba可以通过间歇性随机选取的有限个数的点的平均值来进行估计。这就是数值积分的 Monte Carlo 方法。 MCM 已被成功地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以及其它竞赛活动中都会出现。Monte Carlo 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。8.2 随机数和随机变量的产
4、生510全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数 ，使其将整数变换为随机数。以某种方法选)(xg取 ，并按照 产生下一个随机数。最一般的方程 具有如下形式：0x)(1kkx )(xgmcaxgod)()（8.1）其中 初始值或种子（ ）0x0x乘法器（ ）a0a增值（ ）cc模数m对于 数位的二进制整数，其模数通常为 。例如，对于 31 位的计算机 即可t t2m取 。这里 和 都是整数，且具有相同的取值范围 。132ax,0 0,xca所需的随机数序 便可由下式得 nmcaxnnod)(1（8.2）该序列称为 线性同余序列 。例如，若 且 ，则该序列
5、为70107，6，9，0，7，6，9，0 （8.3）可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。 （8.3）式中序列周期长度为 4。当然，一个有用的序列必是具有相对较长周期的序列。许多作者都用术语 乘同余法 和 混合同余法 分别指代 和 时的线性同余法。选取 和 的法则可参见6,10。0ccax,0m这里我们只关心在区间 内服从均匀分布的随机数的产生。用字符 来)1,( U表示这些数字，则由式（8.2）可得mxUn1（8.4）这样 仅在数组 中取值。 （对于区间（0,1）内的随机Um/)(,./21,0数，一种快速检测其随机性的方法是看其均值是否为
6、0.5。其它检测方法可参见3,6。 ）产生区间 内均匀分布的随机数 ，可用下式),(baXUab)(（8.5）用计算机编码产生的随机数（利用式（8.2）和（8.4） ）并不是完全随机的；事实上，给定序列种子，序列的所有数字 都是完全可预测的。一些作者为强调这一点，将这种计算机产生的序列称为 伪随机数 。但如果适当选取 和 ，ca,m序列 的随机性便足以通过一系列的统计检测。它们相对于真随机数具有可快U速产生、需要时可再生的优点，尤其对于程序调试。Monte Carlo 程序中通常需要产生服从给定概率分布 的随机变量 。)(xFX该步可用6， 1315中的几种方法加以实现，其中包括 直接法 和
7、舍去法 。直接法（也称反演法或变换法） ，需要转换与随机变量 相关的累积概率X函数 （即： 为 的概率） 。 显然表明，)()(xXprobxF)(xFX1)(0xF通过产生(0,1) 内均匀分布随机数 ，经转换我们可得服从 分布的随机样本U。为了得到这样的具有概率分布 的随机数 ，不妨设 ，即可X)(x)(xU得)(1FX（8.6）其中 具有分布函数 。例如，若 是均值为 呈指数分布的随机变量，且X)(xFxex0,1/（8.7）在 中解出 可得)(xU)1ln(UX（8.8）由于 本身就是区间(0,1)内的随机数，故可简写为)1(ln（8.9）有时(8.6) 式所需的反函数 不存在或很难获
8、得。这种情况可用舍去法)(1xF来处理。令 为随机变量 的概率密度函数。令 时的dxf)()(Xbxa，且 上界为 （即： ） ，如图 8.1 所示。我们产生区间0)(xf Mxf)(0,1)内的两个随机数 ，则),(21U1abX MUf21（8.10）分别为在(a,b)和(0,M)内均匀分布的随机数。若)(11Xf（8.11）则 为 的可选值，否则被舍去，然后再试新的一组 。如此运用舍去1X ),(21U法，所有位于 以上的点都被舍去，而位于 上或以下的点都由)(xf )xf来产生 。 11)(UabX1X图 8.1 舍去法产生概率密度函数为 的随机变量 )(xf例 8.1 设计一子程序使
9、之产生 0,1 之间呈均匀分布的随机数 。用该程序产U生随机变 ，其概率分布由下式给定0),cos1(2)T解：生成 的子程序如图 8.2 所示。该子程序中，U且 。应用种子数（如 1234） ，主程,021478361cm6875a序中每调用一次子程序，就会生成一个随机数 。种子数可取 1 到 间的任一Um整数。0001 C*0002 C PROGRAM FOR GENERATING RANDOM VARIABLES WITH0003 C A GIVEN PROBABILITY DISTRIBUTION0004 C*0005 0006 DOUBLE PRECISION ISEED0007 0
10、008 ISEED=1234.D00009 DO 10 I=1,1000010 CALL RANSOM(ISEED,R)0011 THETA=ACOSD(1.0-2.0*R)0012 WRITE(6,*)I,THETA0013 10 CONTINUE0014 STOP0015 END0001 C*0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN0003 C THE INTERVAL (0,1)0004 C*00050006 SUBROUTINE RANDOM (ISEED,R)0007 DOUBLE PRECISION ISDDE,DEL,A
11、0008 DATA DEL,A/2147483647.D0,16807.D0/0009 0010 ISDDE=DMOD(A*ISDDE,DEL)0011 R=ISDDE/DEL0012 RETURN0013 END图 8.2 例 8.1 的随机数生成器 图 8.2 的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有生成随机数的子程序。为了生成随机变量 ，令)cos1(2)TU则有 1U据此，一系列具有给定分布的随机变量 便可由图 8.2 所示主程序中生成。8.3 误差计算Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平均值附近的浮动量，而且不可能
12、达到 100的置信度。要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计13,16。设 是随机变量。则 的 期望值 或 均值 定义为XXxdfx)(（8.12）这里 是 的概率密度分布函数。如果从 中取些独立的随机样本)(xf )(xf，那么的 估计值就表现为 个样本值的均值。N,.21xNnx1（8.13）是 的真正的平均值，而 只是 的有着准确期望值的无偏估计。虽然xX的期望值等于 ，但 。因此，我们还需要 的值在 附近的分布测度。xx x为了估计 以及 在 附近的的值的分布，我们需要引入
13、 的方差，其定X X义为 与 差的平方的期望值，即dxfxVar )()()( 222（8.14）由 ，故有22)(xx dxfdxfdf )()()()( 2（8.15）或者 22)(x（8.16）方差的平方根称为 标准差 ，即 2/12)()x（8.17）标准差给出了 在均值 附近的分布测度，并由此给出了误差幅度的阶数。 的x x标准差与 的标准差的关系表示为Nx)(（8.18）该式表明，如果用根据(8.13)式由 个 值构成的 来估计 ，那么结果中 在nxxx附近的扩散范围便与 成比例，且随着样本数 的增加而降低。x)(x为了估计 的扩散范围，我们定义 样本方差 为NnxS122)(（8
14、.19）由此式还可看出 的期望值等于 。因此样本方差是 的无偏估计。2)(2)(2x将（8.19）式中平方项乘出来，便可得样本标准差为2/1122/1NnxS（8.20）当 较大时，系数 可设为 1。N)/(N作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数MNqpMB)!(!)(（8.21）该式表明 次独立随机试验中有 次成功的概率。(8.21) 中， 是一次试验中N p成功的概率，且 。当 和 都较大时，便可用 Stirling 公式pq1nen2!（8.22）因此(8.21)式可近似为正态分布17)(2exp)(1)(xfMB(8.23)其中 且 。因此当 时，中心极限
15、定理表明，描述由pNxpqx)(N点 Monte Carlo 算法获得的 的分布的概率密度函数是（8.23）式所示的正态x分布函数 。也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）)(f式代入（8.23）式可得)(2exp)(12)( xNxf （8.24）正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例 8.2 中我们给出了根据中心极限定理产生高斯随机变量的算法。由于样本数 是有限的，所以 Monte Carlo 算法不可能达到绝对的确定性。N故而
16、在 附近估计出某一范围或区间以确保我们估计的 落入该范围内。假设要x x得到 位于 和 之间的概率。由定义xxdfxob)(Pr（8.25）令 可得xN/2/202Pr Ndexobxrf/（8.26a）或 1Pr 2/2/ NzxNzxob（8.26b）其中 是误差函数，且 是标准正态差的上 分位点。对（8.26）式可)(xerf 2/z/做如下解释：如果用来呈现独立随机观测值并构建相关随机区间 的 Monte xCarlo 程序以较大的 值反复运行，则这些随机区间的 分位点将NNerf2近似等于 。随机区间 称为置信区间， 称为置信度。大多xxxerf数的 Monte Carlo 算法都使
17、用误差 ，这表明 是在实际均值 的标Nx/x准差范围内的。由（8.26）式可得，样本均值 位于区间 内的概率Nx/是 0.6826 或 68.3。若要求更高的置信度，可用两个或三个标准差的值。如,97.054,682.)()(Pr NxMxNxob 321M（8.27）其中 是标准差的个数。M在（8.26）和（8.27）式中均假设总体标准差 为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差 来估计 的值，从而使学生S氏 分布取代正态分布。我们知道当 很大，比如 时， 分布近似趋t N30t于正态分布。 （8.26）式等价于 1Pr ;2/1;2/ txStxobN（8.2
18、8）其中 为自由度是 的学生氏 分布的上 分位点；其值在任何标1;2/Nt t/准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出上限 NStx1;2/（8.29）下限 tx1;2/（8.30）Monte Carlo 算法中关于误差估计的进一步讨论，可参见18,19。例 8.2 利用中心极限定理产生一高斯（正态）分布的随机变量 。根据中X心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数 ，均值NiY,.21,为 ，方差为 ，Y)(YVarYVarNZi1（8.31）渐渐与 合并为零均值、单位标准差的正态分布。若
19、 是均匀分布的随机变量Ni（即 ） ，则 ，故而iiUY12/)(,2/1YVar12/NUZii（8.32）且变量X（8.33）近似等于均值为 、方差为 的正态变量。 小于 3 时近似为我们所熟知的钟2N形高斯分布。为简化计算，通常实际中设 ，因为这样可消除（8.32）式12中的平方根项。然而 的这种取值截掉了 边界处的分布，且无法产生超过N6的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯3分布（参见20 22 ） 。这样，要产生一个均值为 、标准差为 的高斯随机变量 ，就要遵循以X下步骤：（1） 生成 12 个均匀分布的随机数 121,.U（2） 求得 612iiUZ（
20、3） 令 X8.4 数值积分对于一维积分，现已有一些求积公式，如 3.10 节中所述。而对多重积分的公式则相对较少。Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常复杂，而多重积分的 MCM 几乎保持不变。Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法23。这里将讨论两类 Monte Carlo 积分方法，即简易MCM 和具有对立变量的 MCM。其它类型，如成功失败法和控制变量法，可参见24 26。本节还将简要介绍 MCM
21、 在不规范积分中的应用。8.4.1 简易 Monte Carlo 积分设要计算积分RfI（8.34）其中 是 维空间。令 是在 内均匀分布的随机变量。则RnnX,.21也是随机变量，其均值由下式给出27,28)(XfRIff1)(8.35)方差为221RRffXfVar（8.36）其中 RdX（8.37）如果取 的 个独立样本，即 ，它们与 具有相同的分布，且构XNNX,.21成平均项iiffff 121.（8.38）我们便会期望该平均项接近于 的均值。这样，由（8.35）和（8.38）式，)(Xf即可得NiifRI1（8.39）Monte Carlo 公式可用于有限区域 上的任何积分。为了举
22、例说明，现将R（8.39）式应用于一维和二维积分中。对于一维积分，设badxfI)(（8.40）由（8.39）式可得NiiXfbI1（8.41）其中 是区间 内随机变量，即iXba,10,UabXi（8.42）对于二维积分，有badcdXfI2121,（8.43）则相应的 Monte Carlo 公式为NiifdbI12,)(（8.44）其中10,)(2211 UcdXabi（8.45）（8.39）式中无偏估计值 收敛的很慢，这是由于估计值方差的级数是I。一种改良的方法，即控制变量法，可通过减小估计值方差来提高其准确N/1性和收敛性。8.4.2 具有对立变量的 Monte Carlo 积分方法
23、术语 对立变量 29,30用来描述任一系列估计值，它们能互相抵消掉彼此的方差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值10)(dUgI（8.46）我们期望表达式 与 相比具有更小的方差。如果 太)1()2Ug)( )(Ug小，那反过来 就很可能太大。因此，我们定义估计值(Ni iiUgI112（8.47）其中 是 0 和 1 之间的随机数。此估计值方差的级数是 ，这是在（8.39）iU 4N式基础上的一个极大的进步。对于两维积分10212,dUgI（8.48）则相应的估计值为 Ni iiiiii ggUgI1 21212121 ,4（8.49）根据相似性，可将该
24、方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如（8.41）式到（8.45）式的转换。例如，10dUgabdxfbaNi ii112（8.50）其中 ，且 。由（8.47）和（8.49）式可得，当积分)(XfUgUab维数增加时，每一维用来在简易 Monte Carlo 方法上提高效率的对立变量的最小值也会增加。这样使得简易 Monte Carlo 方法在多维运算中更具优越性。8.4.3 不规范积分积分式0)(dxgI（8.51）可用 Monte Carlo 仿真法进行估计31。对于概率密度为 的随机变量 ，)(xfX其中 在区间(0, )上的积分等于 1，则有)(xf00)(dx
25、gxf（8.52）因此，为计算（8.51）式中的 值，首先要得到 个服从同一在区间(0, )上的IN积分等于 1 的概率密度函数 的独立随机变量。其样本均值)(xfiixfg1（8.53）便是 的估计值。I例 8.3 用 Monte Carlo 方法计算积分102cosdeIj解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得 Monte Carlo 解的精确性。其闭合解为12JI其中 是一阶 Bessel 函数。1J图 8.3 给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其
26、中以及 。该程序在 Vax 11/780 中调用子程序 RANDU 来产0,cba2d生随机数 和 。对于不同的 值，简易和对立变量 Monte Carlo 方法都可1U2N用于计算辐射积分。表 8.1 中对 时的结果进行了精确的比较。在应用5（8.49）式时，用到了下面的对应式：11121, UabXUX222cd表 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.46
27、82-j0.00824000 -0.3999-j0.0345 -0.4216-j0.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+j08.5 位势问题的解位势理论与布朗运动（或随机游动）的关系是由 Kakutani 在 1944 年首次提出的32。自此，所谓的概率位势理论便在诸如热传导3338、静电学3946以及电力工程等许多学科领域得到广泛应用。不同方程式的概率解或 MC解的一个基本概念就
28、是随机游动。不同类型的随机游动应用不同的 Monte Carlo方法。最常见的类型是 固定随机游动 和 自动定位随机游动 。其它不常见类型有迁离法 、 缩减边界法 、 内切形法 以及 表面密度法 。0001 C*0002 C INTEGRATION USING CRUDE MONTE CARLO 0003 C AND ANTITHETIC METHODS0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM0005 C FOR ANY MULTIDIMENSIONAL INTEGRATION0006 C*00070008 DATA I
29、S1,IS2,IS3,IS4/1234,5678,9012,3456/0009 DATA A,B,C/0.0,1.0,0.0/00100011 C SPECIFY THE INTEGRAND0012 F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI)00130014 J=(0.0,1.0)0015 ALPHA=5.00016 PIE=3.14159270017 D=2.0*PIE0018 DO 30 NRUN=500,10000,500 0019 SUM1=(0.0,0.0)0020 SUM2=(0.0,0.0)0021 DO 10 I=1,NRUN0022 CA
30、LL RANDU(IS1,IS2,U1)0023 CALL RANDU(IS3,IS4,U2)0024 X1=A+(B-A)*U10025 X2=C+(D-C)*U20026 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 CONTINUE0031 AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN)0032 AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN)0
31、033 PRINT *,NRUN,AREA1,AREA20034 WRITE(6,*) NRUN,AREA1,AREA20035 WRITE(6,20) NRUN,AREA1,AREA20036 20 FORMAT(2X,NRUN=,I5,3X,AREA1=,F12.6,3X,F12.6,AREA2=,0037 1 F12.6,3X,F12.6,/)0038 30 CONTINUE0039 STOP0040 END图 8.3 例 8.3 中用 Monte Carlo 方法计算二维积分的程序8.5.1 固定随机游动为具体起见，假设用固定随机游动的 MCM 解拉普拉斯方程（在区域 内） 02VR（
32、8.54a）满足 Dirichlet 边界条件（在边界 上） P B（8.54b）首先将 分成网状结构，并用其有限差当量替代 。 （8.54a）在二维 内的有R2R限差表达式可由（3.26）式给出，即 yxVpyxVpyxVpyxVpyx ,（8.55a）其中41yx（8.55b）（8.55）式中，假设网络的一个方格尺寸是 ，如图 8.4 所示。下面给出该方程的概率解释。若随机游动粒子在某一瞬时处于 位置处，它从此点向yx,及 移动的概率分别是 和 。),(),(yx)(yx, yxp,决定粒子移动方式的方法是产生一个随机数 ，并按下面的方式指导10,U粒子的随机游动：yx, 175.0.2U
33、（8.56）如果不用方格而用矩形格，则有 且 ，但 。在用xpyyxp立方体晶格表示的三维问题中， 。两种情61zzyx况中都依据概率对区间 进行了分组。10U图 8.4 固定随机游动示意图为了计算 处的位势，让一随机游动粒子在该点出发开始游动。粒子便yx,开始从一点到另一点在网格中蜿蜒前行直到到达边界。此时，粒子终止游动，记下边界点的指定位势 。令第一个粒子游动结束时 的值记为 ，如图 8.4PVPV1P所示。然后将第二个粒子从 释放，使其游动直到到达边界点，游动终止且yx,该点 的值相应的记为 。依次第三个，第四个， ，第 个粒子从PV2P N释放重复此过程，并记下相应的指定位势 ， 。根
34、据yx, 4,3P PKakutani 32， ， 的期望值是 Dirichlet 问题在 的解，,1PV NP yx,即iPVyx1,（8.57）其中 较大，为游动总次数。其收敛速度随 而改变，所以需要保证许多随NN机游动都有精确的结果。若要解 Poisson 方程在 内 yxgV,2R（8.58a）且 满足条件V在 上 PVB（8.58b）则式（8.55）中的有限差分表示为yxpyxpyxV, 4, 2gVyy （8.59）其中概率仍为式（8.55b）所示。 （8.59）式的概率解释也近似于（8.55）式。但随机游动的每一步都要记录下（8.59）式中的 项。如果第 次随机游动从/2gi出发
35、至到达边界需要 步，则记录下yx, imimjjPyxgiV12,4（8.60）由此可得， 的 Monte Carlo 解为yxV, NimjjNiPiyxgV121 ,4,（8.61）对刚才所描述 MCM 的一个有趣的分析就是游走醉鬼问题15,35。我们将随机游动粒子当作一个“醉鬼” ，将网格方块当作“城市街区” ，结点当作“十字路口” ，边界 当作“城市边界” ，而且把边界 上的结点当作“警察” 。尽管BB醉鬼想步行走回家，但他却醉的一塌糊涂以至于在整个城市中随机游走。警察的工作是在城市边界上一看见醉鬼便将其抓获，并勒令醉鬼交付罚金 。那么PV醉鬼要支付的预期罚金将会是多少呢？其答案便是（8.57）式。在电介质边界上，指定边界条件为 。nD21yx,xxxXXXX1.R.Hersch and R.J.Griego,“Brownian motion and potential theory,”Sci.Amer.,Mar.1969,pp.67-74.
