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 随机数和随机变量的产生5
4、10全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数 ，使其将整数变换为随机数。以某种方法选取 ，并按照)(xg 0x产生下一个随机数。最一般的方程 具有如下形式：)(1kkx )(xgmcaxod)()（8.1）其中 初始值或种子（ ）0x0x乘法器（ ）a增值（ ）c模数m对于 数位的二进制整数，其模数通常为 。例如，对于 31 位的计算机 即可取 。t t2m132这里 和 都是整数，且具有相同的取值范围 。所需的随机数序ax,0c 0,xcam便可由下式得 nxmcaxnnod)(1（8.2）该序列称为 线性同余序列 。例如，若 且 ，则该序列为701
5、07，6，9，0，7，6，9，0 （8.3）可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。 （8.3）式中序列周期长度为 4。当然，一个有用的序列必是具有相对较长周期的序列。许多作者都用术语 乘同余法 和 混合同余法 分别指代 和 时的线性同余法。0c选取 和 的法则可参见6,10。cax,0m这里我们只关心在区间 内服从均匀分布的随机数的产生。用字符 来表示这些)1,0( U数字，则由式（8.2）可得mxUn1（8.4）这样 仅在数组 中取值。 （对于区间（0,1）内的随机数，一U/)(,./21,0种快速检测其随机性的方法是看其均值是否为 0.5。
6、其它检测方法可参见3,6。 ）产生区间内均匀分布的随机数 ，可用下式),(baXUab)(（8.5）用计算机编码产生的随机数（利用式（8.2）和（8.4） ）并不是完全随机的；事实上，给定序列种子，序列的所有数字 都是完全可预测的。一些作者为强调这一点，将这种计算机产生的序列称为 伪随机数 。但如果适当选取 和 ，序列 的随机性便足以通过一ca,m系列的统计检测。它们相对于真随机数具有可快速产生、需要时可再生的优点，尤其对于程序调试。Monte Carlo 程序中通常需要产生服从给定概率分布 的随机变量 。该步可用)(xFX6，1315 中的几种方法加以实现，其中包括 直接法 和 舍去法 。直
7、接法（也称反演法或变换法） ，需要转换与随机变量 相关的累积概率函数（即： 为 的概率） 。 显然表明，通过产生)()(xXprobxF)(xFX1)(0x(0,1)内均匀分布随机数 ，经转换我们可得服从 分布的随机样本 。为了得到这样U)(FX的具有概率分布 的随机数 ，不妨设 ，即可得)(xxU)(1UFX（8.6）其中 具有分布函数 。例如，若 是均值为 呈指数分布的随机变量，且X)(xFxex0,1/（8.7）在 中解出 可得)(xU)1ln(UX（8.8）由于 本身就是区间(0,1)内的随机数，故可简写为)1(ln（8.9）有时(8.6)式所需的反函数 不存在或很难获得。这种情况可用
8、舍去法来处理。令)(1xF为随机变量 的概率密度函数。令 时的 ，且 上界dxFf)()(Xbxa0)(xf)(xf为 （即： ） ，如图 8.1 所示。我们产生区间(0,1)内的两个随机数 ，Mf ,21U则11)(Uab Mf21（8.10）分别为在(a,b)和(0,M)内均匀分布的随机数。若)(11Xf（8.11）则 为 的可选值，否则被舍去，然后再试新的一组 。如此运用舍去法，所有1X ),(21U位于 以上的点都被舍去，而位于 上或以下的点都由 来产生)(xf )(xf 1)(abX。 1图 8.1 舍去法产生概率密度函数为 的随机变量 )(xf例 8.1 设计一子程序使之产生 0,
9、1 之间呈均匀分布的随机数 。用该程序产生随机变 ，U其概率分布由下式给定0),cos1(2)T解：生成 的子程序如图 8.2 所示。该子程序中， 且U ,021478361cm。应用种子数（如 1234） ，主程序中每调用一次子程序，就会生成一个随168075a机数 。种子数可取 1 到 间的任一整数。m0001 C*0002 C PROGRAM FOR GENERATING RANDOM VARIABLES WITH0003 C A GIVEN PROBABILITY DISTRIBUTION0004 C*0005 0006 DOUBLE PRECISION ISEED0007 0008
10、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,A0008
11、 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 程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平均值附近的浮动量，而且不可能达到 1
12、00的置信度。要计算 Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计13,16。设 是随机变量。则 的 期望值 或 均值 定义为XXxdfx)(（8.12）这里 是 的概率密度分布函数。如果从 中取些独立的随机样本 ，)(xf )(xf Nx,.21那么的 估计值就表现为 个样本值的均值。NNnx1（8.13）是 的真正的平均值，而 只是 的有着准确期望值的无偏估计。虽然 的期望值xX x等于 ，但 。因此，我们还需要 的值在 附近的分布测度。x为了估计 以及 在 附近的的值的分布，我们需要引入
13、 的方差，其定义为 与x XX差的平方的期望值，即dxfxVar )()()( 222（8.14）由 ，故有22)(xx dxfdxfdf )()()( 2（8.15）或者 22)(x（8.16）方差的平方根称为 标准差 ，即 2/12)()x（8.17）标准差给出了 在均值 附近的分布测度，并由此给出了误差幅度的阶数。 的标准差与x 的标准差的关系表示为xNx)(（8.18）该式表明，如果用根据(8.13)式由 个 值构成的 来估计 ，那么结果中 在 附近的nxxx扩散范围便与 成比例，且随着样本数 的增加而降低。)(x为了估计 的扩散范围，我们定义 样本方差 为NnxS122)(（8.19
14、）由此式还可看出 的期望值等于 。因此样本方差是 的无偏估计。将2)(2)(2x（8.19）式中平方项乘出来，便可得样本标准差为2/1122/1NnxS（8.20）当 较大时，系数 可设为 1。N)/(N作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数MNqpMB)!(!)(（8.21）该式表明 次独立随机试验中有 次成功的概率。(8.21)中， 是一次试验中成功的概率，Np且 。当 和 都较大时，便可用 Stirling 公式pq1nen2!（8.22）因此(8.21)式可近似为正态分布17(8.23)(2exp)(1)(xfMB其中 且 。因此当 时，中心极限定理表明，描
15、述由 点pNxpqxNNMonte Carlo 算法获得的 的分布的概率密度函数是（8.23）式所示的正态分布函数 。)(xf也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得)(2exp)(12)( xxf （8.24）正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例 8.2 中我们给出了根据中心极限定理产生高斯随机变量的算法。由于样本数 是有限的，所以 Monte Carlo 算法不可能达到绝对的确定性。故而在N附近估计出某一
16、范围或区间以确保我们估计的 落入该范围内。假设要得到 位于x x x和 之间的概率。由定义xdfxob)(Pr（8.25）令 可得xN/2/202Pr Ndexobxerf/（8.26a）或 1Pr 2/2/ NzxNzxob（8.26b）其中 是误差函数，且 是标准正态差的上 分位点。对（8.26）式可做如下解)(xerf 2/z/释：如果用来呈现独立随机观测值并构建相关随机区间 的 Monte Carlo 程序以较大x的 值反复运行，则这些随机区间的 分位点将近似等于 。随机区间NNerfx称为置信区间， 称为置信度。大多数的 Monte Carlo 算法都使用误xxNerf2差 ，这表明
17、 是在实际均值 的标准差范围内的。由（8.26）式可得，样本/均值 位于区间 内的概率是 0.6826 或 68.3。若要求更高的置信度，可用xx/两个或三个标准差的值。如,97.054,682.)()(Pr NxMxNxob 321M（8.27）其中 是标准差的个数。M在（8.26）和（8.27）式中均假设总体标准差 为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差 来估计 的值，从而使学生氏 分布取代正态St分布。我们知道当 很大，比如 时， 分布近似趋于正态分布。 （8.26）式等价N30t于 1Pr ;2/1;2/ NtxtxobN（8.28）其中 为自由度是 的
18、学生氏 分布的上 分位点；其值在任何标准统计学1;2/Nt t/课本中都有表列可查。这样置信区间的上、下限便可由下式给出上限 NStx1;2/（8.29）下限 tx1;2/（8.30）Monte Carlo 算法中关于误差估计的进一步讨论，可参见18,19 。例 8.2 利用中心极限定理产生一高斯（正态）分布的随机变量 。根据中心极限定理，X大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数 ，均值为 ，方差为 ，NiY,.21,Y)(YVarVarZi1（8.31）渐渐与 合并为零均值、单位标准差的正态分布。若 是均匀分布的随机变量（即
19、NiY） ，则 ，故而iiUY12/)(,2/1YVar12/NUZii（8.32）且变量X（8.33）近似等于均值为 、方差为 的正态变量。 小于 3 时近似为我们所熟知的钟形高斯分2N布。为简化计算，通常实际中设 ，因为这样可消除（8.32）式中的平方根项。然而1的这种取值截掉了 边界处的分布，且无法产生超过 的值。对于分布曲线的两端N6比较重要的仿真，就必须用其它方案来产生高斯分布（参见2022 ） 。这样，要产生一个均值为 、标准差为 的高斯随机变量 ，就要遵循以下步骤：X（1） 生成 12 个均匀分布的随机数 121,.U（2） 求得 612iiUZ（3） 令 X8.4 数值积分对于
20、一维积分，现已有一些求积公式，如 3.10 节中所述。而对多重积分的公式则相对较少。Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常复杂，而多重积分的 MCM 几乎保持不变。Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法23。这里将讨论两类 Monte Carlo 积分方法，即简易 MCM 和具有对立变量的 MCM。其它类型，如成功失败法和控制变量法，可参见2426。本节还将简要介绍 MCM 在不规范积分中的应用。8.4.
21、1 简易 Monte Carlo 积分设要计算积分RfI（8.34）其中 是 维空间。令 是在 内均匀分布的随机变量。则 也RnnX,.21 )(Xf是随机变量，其均值由下式给出27,28(8.35)RIfXf1)(方差为221RRfffVar（8.36）其中 RdX（8.37）如果取 的 个独立样本，即 ，它们与 具有相同的分布，且构成平均XNNX,.21项iiffff 121.（8.38）我们便会期望该平均项接近于 的均值。这样，由（8.35）和（8.38）式，即可得)(XfNiifRI1（8.39）Monte Carlo 公式可用于有限区域 上的任何积分。为了举例说明，现将（8.39）式应用于一维和二维积分中。对于一维积分，设badxfI)(（8.40）由（8.39）式可得NiiXfbI1（8.41）其中 是区间 内随机变量，即iXba,10,UabXi（8.42）对于二维积分，有badcdXfI2121,（8.43）
