蒙特卡罗模拟

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

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

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

4、论述了产生随机数的各类方法。 其中较为普遍应用的产生随机数的方法 是选取一个函数g (x)，使其将整数变换为随机数。以某种方法选取X。,并按照xkA = g (xk) 产生下一个随机数。最一般的方程 g (x)具有如下形式：g (x) = (ax c) mod m( 8.1)其中x。二初始值或种子。0 )a 乘法器9 - 0) c =t曾值(C - 0 ) m二模数对于t数位的二进制整数，其模数通常为2t。例如，对于31位的计算机m即可取231,o这里x0，a和 c都是整数，且具有相同的取值范围m a,m c, m x。所需的随机数序：x便可由下式得该序列称为 线性同余序列。例如，若X。=a

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

6、其随机性的方法是看其均值是否为0.5。其它检测方法可参见3,6。)产生区间(a,b)内均匀分布的随机数X，可用下式X 二 a (b -a)U( 8.5)用计算机编码产生的随机数(利用式(8.2)和(8.4)并不是完全随机的；事实上，给定序列种 子，序列的所有数字U都是完全可预测的。一些作者为强调这一点，将这种计算机产生的序列称为伪随机数。但如果适当选取a,c和m，序列的随机性便足以通过一系列的统计检测。它们相对于真随机数具有可快速产生、需要时可再生的优点，尤其对于程序调试。Monte Carlo程序中通常需要产生服从给定概率分布F (x)的随机变量X。该步可用)，13 - 15中的几种方法加以

7、实现，其中包括直接法和舍去法。直接法(也称反演法或变换法)，需要转换与随机变量X相关的累积概率函数F(x)= prob(X - x)(即：F (x)为X x的概率)。0 F (x) 1显然表明，通过产生(0,1)内均匀分布随机数U，经转换我们可得服从F(x)分布的随机样本x%y得到这样的具有概率分布F(x)的随机数X，不妨设U = F(x)，即可得X 二 F U)( 8.6)其中X具有分布函数F(x)。例如，若X是均值为呈指数分布的随机变量，且.c 项 L1(8.7)F(x) =1 e ,0 : (8.7)在U二F(x)中解出X可得(8.8)由于(1 -U )本身就是区间(0,1)内的随机数，

8、故可简写为(8.9)(8.10)(8.9)(8.10)有时(8.6)式所需的反函数F J(x)不存在或很难获得。这种情况可用舍去法来处理。令 f(x) =dF凶为随机变量X的概率密度函数。令a x b时的f (x) = 0,且f (x)上界 dx 为M (即:f(x)兰M)，如图8.1所示。我们产生区间(0,1)内的两个随机数(Ui，U2)，则X.= a (b-a) f1 =U2Mfdf(X1)分别为在(a,b )和fdf(X1)(8.11)则Xi为X的可选值，否则被舍去，然后再试新的一组(山2 2)。如此运用舍去法，所有位于f (x)以上的点都被舍去，而位于f (x)上或以下的点都由XAa

9、- (b-a)U来产生X。图8.1舍去法产生概率密度函数为f(x)的随机变量例8.1设计一子程序使之产生0,1之间呈均匀分布的随机数 U。用该程序产生随机变0，其概率分布由下式给定)2 (1 COST ),0 : V :解：生成U的子程序如图8.2所示。该子程序中，m =221-1二2147483647, 0且a =7=16807。 应用种子数(如1234)，主程序中每调用一次子程序，就会生成一个随机 数。种子数可取1到m间 的任一整数。00010002*CPROGRAM FOR GENERA TING RANDOM V ARIABLES WITH0003CA GIVEN PROBABILIT

10、Y DISTRIBUTION0004AC*00050006000700080009001000110012001310DOUBLE PRECISION ISEEDISEED=1234.D0DO 10I=1,100CALL RANSOM(ISEED,R)THETA=ACOSD(1.0-2.0*R)WRITE(6,*)I,THETACONTINUE001400150001ACSTOP END*0002CSUBROUTINE FOR GENERATING RANDOM NUMBERS IN0003CTHE INTERV AL (0,1)0004AC*00050006SUBROUTINE RANDOM

11、 (ISEED,R)0007DOUBLE PRECISION ISDDE,DEL,A0008DATA DEL,A/2147483647.D0,16807.D0/00090010ISDDE=DMOD(A*ISDDE,DEL)0011R=ISDDE/DEL0012RETURN0013END图8.2例8.1的随机数生成器图8.2的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有生成随机数的子程序。为了生成随机变量0，令U =T(0)(1-cos0)2则有-T -(u) =cos(1 -2u )据此，一系列具有给定分布的随机变量0便可由图8.2所示主程序中生成。8.3误差计算Monte Car

12、lo程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平均值附近的浮动量，而且不可能达到100%的置信度。要计算Monte Carlo算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计13,16。设X是随机变量。则X的期望值或均值x定义为_ O0(8.12)X = xf (x) dx j-=0(8.12)这里f (X)是X的概率密度分布函数。如果从f(x)中取些独立的随机样本XN，那么的x估计值就表现为N个样本值的均值。(8.13)x是X的真正的平均值，而5?只是x的有着准确期望值的无偏估计。虽然于x，但？ -

13、x。因此，我们还需要？的值在x附近的分布测度。为了估计X以及刃在x(8.13)x是X的真正的平均值，而5?只是x的有着准确期望值的无偏估计。虽然于x，但？ -x。因此，我们还需要？的值在x附近的分布测度。为了估计X以及刃在x附近的的值的分布，我们需要引入x的方差，其定义为X与x差的平 方的期望值，即5?的期望值2Z2、一 2Var(x) (x - x) (x_x) f (x)dx(8.14)由板 _x)2 = X2-2xx - 乂2，故有2(x) X2 f (x)dx2x xf(x)dx X2 f (x)dx(8.15)或者二 2(X)=X2X2(8.16)方差的平方根称为标准差，即匚(x)

14、=(x2 - X2)1/2(8.17)标准差给出了 x在均值x附近的分布测度，并由此给出了误差幅度的阶数。 的标准差的关系表示为?的标准差与X-(x)5 N(8.18)该式表明，如果用根据(8.13)式由N个xn值构成的？来估计乂，那么结果中？在x附近的扩散范围便与二(X)成比例，且随着样本数 N的增加而降低。为了估计？的扩散范围，我们定义 样本方差为S1 优-幻2N - 1 n d(8.19)由此式还可看出S2的期望值等于2(x)。因此样本方差是二2(x)的无偏估计。将(8.19)式中平方 项乘出来，便可得样本标准差为/ N、1/ 2 1/2IN 1 J(8.20)当N较大时，系数N/(N-

15、1) 可设为1。(8.13)作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数B(M )二N!M!(N -M)!(8.21B(M )二N!M!(N -M)!(8.21)该式表明N次独立随机试验中有M次成功的概率。（8.21）中，P是一次试验中成功的概率,且q=1-p。当M和N-M都较大时，便可用Stirli ng公式n! nne.2 二 n（8.22）因此（8.21）式可近似为正态分布171 厂 exp （?-刃2B（M ）: f （X =（?） .一 2 二 tUU（8.23）CT（幻其中X二N p且二（X） Npq。因此当N：时，中心极限定理表明，描述由N点MonteCar

16、lo算法获得的刃的分布的概率密度函数是（Carlo算法获得的刃的分布的概率密度函数是（8.23）式所示的正态分布函数f（X）。也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得(8.24)N（） ? X）(8.24)2心）正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。由于样本数N是有限的，所以Monte Carlo算法不可能达到绝对的确定性。故而在x附近估计出某 一范围或

17、区间以确保我们估计的？落入该范围内。假设要得到？位于X和X :之间的概率。由定义ProbX -： ：）? : X f （）c） d?（ 8.25）*Q令一 lx可得.2/N 二 xPr obx 一；(8.26a)j d(8.26a)N/2x )J_ CTProbx _ CTProbx -乙仞？兰？兰 x + ZG/2kVN(8.26b)其中erf （x）是误差函数，且z./2是标准正态差的上:/2分位点。对（8.26）式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间x _ :.的Monte Carlo程序以较大的N值反复运行，则这些随机区间的e rfN？。随机区间x_、：称 2 b

18、（x）J为置信区间，曰f称为置信度。大多数的释：如果用来呈现独立随机观测值并构建相关随机区间x _ :.的Monte Carlo程序以较大的N值反复运行，则这些随机区间的e rfN？。随机区间x_、：称 2 b（x）J为置信区间，曰f称为置信度。大多数的Monte Carlo算法都使用误差侦 2 o（x ）J、二匚X/N，这表明？是在实际均值x的标准差范围内的。由（8.26）式可得，样本均I值外位于区间、二X / , N内的概率是0.6826或68.3%。若要求更高的置信度，可用两个或三个标准差的值。如ProbX M -（x） .N其中M是标准差的个数。0.6826,M=1:、坯M= 2095

19、4,M=30.997,(8.27)在（8.26）和（8.27）式中均假设总体标准差二为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差S来估计匚的值，从而使学生氏t 一分布取代正态分布。我们知道当N很大，比如N -30时，t 分布近似趋于正态分布。（8.26）式等价于Probj x _ SL,2；n j(8.28)其中t-./2； N为自由度是N -1的学生氏t 一分布的上:/2分位点；其值在任何标准统计学课本中都 有表列可查。这样置信区间的上、下限便可由下式给出S2；(8.29)IN 4、NSt：./2;N下限二x1V NMonte Carlo算法中关于误差估计的进一步讨

20、论，可参见18,19。(8.30)X。根据中心极限定理,例8.2利用中心极限定理产生一高斯（正态）分布的随机变量大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数Y=1,2,.N，均值为Y，方差为Var （Y）,NY -nY_目(8.31)NVar Y渐渐与N合并为零均值、单位标准差的正态分布。若Y.是均匀分布的随机变量（即Y=Uj）,则 Y-1/2,Var（YA1/ .12，故而Nv UN/2Z = i (8.32)“ N /12且变量X - ；Z(8.33)近似等于均值为、方差为二2的正态变量。N小于3时近似为我们所熟知的钟形高斯分布

21、。为简化计算，通常实际中设N=12,因为这样可消除(8.32)式中的平方根项。然而N的这种取 值截掉了 _6匚边界处的分布，且无法产生超过3二的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布(参见20- 22)。这样，要产生一个均值为J、标准差为二的高斯随机变量X，就要遵循以下步骤：生成12个均匀分布的随机数UR?,山212求得 Z 二 U6令X二；之-j 8.4数值积分对于一维积分，现已有一些求积公式，如3.10节中所述。而对多重积分的公式则相对较少。Monte Carlo技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常 复杂，而多重积分的MCM几

22、乎保持不变。Monte Carlo积分的收敛性与维数无关，但对求积公式并 非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维 或三维积分的很有效的方法23。这里将讨论两类Monte Carlo积分方法，即简易MCM和具有对立变量的MCM。其它类型，如成功一失败法和控制变量法，可参见 24 - 26。本节还将简要介绍MCM在不规范积分中的应用。8.4.1 简易 Monte Carlo 积分设要计算积分(8.34)f (X)也是其中R是n维空间。令X二X1,X2,.Xn是在Rf (X)也是1 f(XLR f 只方差为其中R = 0X(8.37)如果取X的N个独

23、立样本，即X1,X2,.,XN，它们与X具有相同的分布，且构成平均项f Xif X2 . f X其中R = 0X(8.37)如果取X的N个独立样本，即X1,X2,.,XN，它们与X具有相同的分布，且构成平均项f Xif X2 . f XNN(8.38)我们便会期望该平均项接近于f（X）的均值。这样由（8.35）和（8.38 ）式，即可得Rn1 N.* 1i4(8.39)Monte Carlo公式可用于有限区域R上的任何积分。为了举例说明，现将（ 一维和二维积分中。对于一维积分，设8.39 ）式应用于b、I=af (X)dX（8.40）由（8.39）式可得b Ynf Xi(8.41)其中*是区间

24、a,b内随机变量Xj = a b-a U ,0 : U : 1对于二维积分，有(8.42(8.36)Var(f(X)气卜肯 f(8.36)(8.43)I = fff(X1,X2dX1dX2(8.43)a c则相应的Monte Carlo公式为(8.44)(b - a) d f 秘(8.44)其中X；二 X；二 a (b ajUi,。：U i : 12 2 2X.二 c (d c)U ,0 : U : 1（8.39）式中无偏估计值1收敛的很慢，这是由于估计值方差的级数是良的1/ No 一种改 方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛性。8.4.2具有对立变量的Monte Ca

25、rlo积分方法术语对立变量29,30用来描述任一系列估计值，它们能互相抵消掉彼此的方差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值I=jg (U) dU(8.46)我们期望表达式xg(U厂g(1 -U)1与g(U)相比具有更小的方差。如果g (U)太小，那 反 过来g(1 u )就很可能太大。因此，我们定义估计值I 1 N 2 g U.g 1-U. |(8.47)1其中U.是 0和1之间的随机数。此估计值方差的级数是幽 这是在(8.39)式基础上的11 11 “ 1 2 1 2 匕.gu，U dU )(8.48dU(8.48则相应的估计值为gU：，1-U.2g1-

26、U：，U.gU：，1-U.2g1-U：，U.2 g 1-U ； ,1 - U :1财.根据相似性，可将该方法延拓至更高阶的积分。对于不是(0，1)的区间，可应用诸如(8.41)式到(8.45)式的转换。例如，；f xdxh b -a 0g U dUb a n 1讥沙 Ui gi(8.50)其中g (U)二f (X)，且X二a - b -a。由(8.47)和(8.49)式可得，当积分维数增加时，每一维用来在简易Monte Carlo方法上提高效率的对立变量的最小值也会增加。这样使得简易Monte Carlo方法在多维运算中更具优越性。8.4.3 不规范积分积分式0I 二 0 g(x)dx(8.5

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

28、中JN 是一阶Bessel函数。图8.3给出由（8.44 ）和（8.45）式计算该积分的一个简单的程序，其中a=0,b=1,C = 0以及d=2二。该程序在Vax 11/780中调用子程序RANDU来产生随机数U 1和U 2。对于不 同的N值，简易和对立变量Monte Carlo方法都可用于计算辐射积分。表8.1中对=5时的结果进行了精确的比较。在应用（8.49）式时，用到了下面的对应式：UAX15UAX2J-UAAXA b-a 1-U11-U2 三 d -X2 二 d -c 1-U2表8.1例8.3辐射积分的 MC方法积分结果5001000-0.2892-j0.0742-0.5737+j0.

29、080820004000-0.4922-j0.0040-0.3999-j0.0345-0.2887-j0.0585-0.4982-j0.0080-0.4682-j0.0082-0.4216-j0.03236000800010,0005001000-0.2892-j0.0742-0.5737+j0.080820004000-0.4922-j0.0040-0.3999-j0.0345-0.2887-j0.0585-0.4982-j0.0080-0.4682-j0.0082-0.4216-j0.03236000800010,000-0.3608-j0.0270-0.4327-j0.0378-0.42

30、29-j0.0237-0.3787-j0.0440-0.4139-j0.0241-0.4121-j0.0240精确解：-0.4116+j08.5位势问题的解位势理论与布朗运动（或随机游动）的关系是由 Kakutani在1944年首次提出的32。自此，所谓的概率位势理论便在诸如热传导33 - 38、静电学39 - 46以及电力工程等许多学科领域得到广泛应用。不同方程式的概率解或MC解的一个基本概念就是随机游动。不同类型的随机游动应用不同的Monte Carlo方法。最常见的类型是固定随机游动和自动定位随机游动其它不常见类型有迁离法、缩减边界法、内切形法以及表面密度法。位随机游动C*C INTEG

31、RATION USING CRUDE MONTE CARLOC AND ANTITHETIC METHODSC ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAMC FOR ANY MULTI -DIMENSIONAL INTEGRATIONC*DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/DATA A,B,C/0.0,1.0,0.0/C SPECIFY THE INTEGRAND F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI)J=(0.0,1.0) ALPHA=5.0PIE

32、=3.1415927D=2.0*PIEDO 30 NRUN=500,10000,500SUM1=(0.0,0.0)SUM2=(0.0,0.0)DO 10 I=1,NRUNCALL RANDU(IS1,IS2,U1)CALL RANDU(IS3,IS4,U2) X1=A+(B-A)*U1 X2=C+(D-C)*U2 X3=(B-A)*(1.0- U1) X4=(D-C)*(1.0-U2) SUM1=SUM1+F(X1,X2)SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4)10 CONTINUE AREA1=(B-A)*(D-C)*SUM1/FLOA T(

33、NRUN) AREA2=(B-A)*(D- C)*SUM2/(4.0*FLOAT(NRUN) PRINT *,NRUN,AREA1,AREA2 WRITE(6,*) NRUN,AREA1,AREA2 WRITE(6,20) NRUN,AREA1,AREA220 FORMAT(2X, NRUN= ,I5,3X, AREA1= ,F12.6,3X,F12.6, AREA2=1 F12.6,3X,F12.6,/)30 CONTINUESTOPEND00010002000300040005000600070008000900100011001200130014001500160017001800190

34、02000210022002300240025002600270028002900300031003200330034003500360037003800390040图8.3例8.3中用Monte Carlo方法计算二维积分的程序8.5.1固定随机游动为具体起见，假设用固定随机游动的2为具体起见，假设用固定随机游动的2V =0满足Dirichlet边界条件V=VP首先将R分成网状结构，并用其有限差当量替代可由(3.26)式给出，即Vx,y 二 PxV x 入 yPXV x - :,(8.55a)其中Px = Px_= Py = Py(8.55)式中，假设网络的一个方格尺寸是LMCM解拉普拉斯方

35、程(在区域R内)(8.54a)(在边界B上)(854b)、2(8.54a)在二维R内的有限差表达式y PyV x, ypAV x, y 1(8.55b)4、，如图8.4所示。下面给出该方程的概率解释。及x, y -L移动的概率分别是若随机游动粒子在某一瞬时处于x, y位置处，它从此点向 危*二，y) , (x -L y及x, y -L移动的概率分别是Px., Px_Py .和Py_。决定粒子移动方式的方法是产生一个随机数U 随机数U ,0 : U :1，并按下面的方式指导粒子的随机游动:X, y x . X, y x . : , y 0 : U : 0.25x, y x-，,y0.25 :：0.5: x, y x,y 0.5 : U