Monte Carlo 法§8.1 概述Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。

Monte Carlo 方法（MCM ），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态[1]。

它是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。

运用该近似方法所获得的问题的解in spirit 更接近于物理实验结果，而不是经典数值计算结果。

普遍认为我们当前所应用的MC 技术，其发展约可追溯至1944年，尽管在早些时候仍有许多未解决的实例。

MCM 的发展归功于核武器早期工作期间Los Alamos （美国国家实验室中子散射研究中心）的一批科学家。

Los Alamos 小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MCM 在各种问题中的应用[2]-[4]。

“Monte Carlo ”的名称取自于Monaco （摩纳哥）内以赌博娱乐而闻名的一座城市。

Monte Carlo 方法的应用有两种途径：仿真和取样。

仿真是指提供实际随机现象的数学上的模仿的方法。

一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动来模仿中子的锯齿形路径。

取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。

例如，在上的平均值可以)(x f b x a <<通过间歇性随机选取的有限个数的点的平均值来进行估计。

这就是数值积分的Monte Carlo 方法。

MCM 已被成功地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。

任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。

这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以及其它竞赛活动中都会出现。

Monte Carlo 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。

§8.2 随机数和随机变量的产生[5]－[10]全面的论述了产生随机数的各类方法。

其中较为普遍应用的产生随机数的方法是选取一个函数，使其将整数变换为随机数。

以某种方法选取)(x g 0x ，并按照产生下一个随机数。

最一般的方程具有如下形式：)(1k k x g x =+)(x gm c ax x g mod )()(+=（8.1）其中 初始值或种子（）=0x 00>x 乘法器（）=a 0≥a 增值（）=c 0≥c模数=m对于数位的二进制整数，其模数通常为。

例如，对于31位的计算机即可取t t 2m 。

这里和都是整数，且具有相同的取值范围。

所1312-a x ,0c 0,,x m c m a m >>>需的随机数序便可由下式得 {}n xm c ax x n n mod )(1+=+（8.2）该序列称为线性同余序列。

例如，若且，则该序列为70===c a x 10=m 7，6，9，0，7，6，9，0…… （8.3）可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。

（8.3）式中序列周期长度为4。

当然，一个有用的序列必是具有相对较长周期的序列。

许多作者都用术语乘同余法和混合同余法分别指代和时的线性同余法。

选取和的法则可参见[6,10]。

0=c 0≠c c a x ,,0m 这里我们只关心在区间内服从均匀分布的随机数的产生。

用字符来)1,0(U 表示这些数字，则由式（8.2）可得mx U n 1-=（8.4）这样仅在数组中取值。

（对于区间（0,1）内的随机数，U {}m m m m /)1(,......,/2,/1,0-一种快速检测其随机性的方法是看其均值是否为0.5。

其它检测方法可参见[3,6]。

）产生区间内均匀分布的随机数，可用下式),(b a XU a b a X )(-+=（8.5）用计算机编码产生的随机数（利用式（8.2）和（8.4））并不是完全随机的；事实上，给定序列种子，序列的所有数字都是完全可预测的。

一些作者为强调U 这一点，将这种计算机产生的序列称为伪随机数。

但如果适当选取和，序c a ,m 列的随机性便足以通过一系列的统计检测。

它们相对于真随机数具有可快速产U 生、需要时可再生的优点，尤其对于程序调试。

Monte Carlo 程序中通常需要产生服从给定概率分布的随机变量。

)(x F X 该步可用[6]，[13]－[15]中的几种方法加以实现，其中包括直接法和舍去法。

直接法（也称反演法或变换法），需要转换与随机变量相关的累积概率函X 数（即：为的概率）。

显然表明，通)()(x X prob x F ≤=)(x F x X ≤1)(0≤≤x F过产生(0,1)内均匀分布随机数，经转换我们可得服从分布的随机样本。

U )(x F X 为了得到这样的具有概率分布的随机数，不妨设，即可得)(x F X )(x F U =)(1U F X -=（8.6）其中具有分布函数。

例如，若是均值为呈指数分布的随机变量，且X )(x F X μ∞<<-=-x e x F x 0,1)(/μ（8.7）在中解出可得)(x F U =X)1ln(U X --=μ（8.8）由于本身就是区间(0,1)内的随机数，故可简写为)1(U -U X ln μ-=（8.9）有时(8.6)式所需的反函数不存在或很难获得。

这种情况可用舍去法来)(1x F -处理。

令为随机变量的概率密度函数。

令时的，dxx dF x f )()(=X b x a >>0)(=x f 且上界为（即：），如图8.1所示。

我们产生区间(0,1)内的两个)(x f M M x f ≤)(随机数，则),(21U U11)(U a b a X -+=M U f 21=（8.10）分别为在(a,b)和(0,M)内均匀分布的随机数。

若)(11X f f ≤（8.11）则为的可选值，否则被舍去，然后再试新的一组。

如此运用舍去法，1X X ),(21U U 所有位于以上的点都被舍去，而位于上或以下的点都由)(x f )(x f 来产生。

11)(U a b a X -+=1X图8.1 舍去法产生概率密度函数为的随机变量)(x f 例8.1 设计一子程序使之产生0,1之间呈均匀分布的随机数。

用该程序产生U 随机变，其概率分布由下式给定Θπθθθ<<-=0),cos 1(21)(T 解：生成的子程序如图8.2所示。

该子程序中，U ,0,21474836471221==-=c m 且。

应用种子数（如1234），主程序中每调用一次子程序，就会生1680775==a 成一个随机数。

种子数可取1到间的任一整数。

U m0001C**********************************************************0002 C PROGRAM FOR GENERATING RANDOM VARIABLESWITH0003 C A GIVEN PROBABILITY DISTRIBUTION 0004C********************************************************** 0005 0006 DOUBLE PRECISION ISEED 0007 0008 ISEED=1234.D0 0009 DO 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 END0001C********************************************************** 0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN0003 C THE INTERVAL (0,1) 0004C********************************************************** 0005 0006 SUBROUTINE RANDOM (ISEED,R) 0007 DOUBLE PRECISION ISDDE,DEL,A 0008 DATA 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的子程序只是为了说明本章所介绍的一些概念。

大多数计算机都有生成随机数的子程序。

为了生成随机变量，令Θ)cos 1(21)(Θ-=Θ=T U 则有)21(cos )(11U U T -==Θ--据此，一系列具有给定分布的随机变量便可由图8.2所示主程序中生成。

Θ§8.3 误差计算Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。

因此，该解中包含了平均值附近的浮动量，而且不可能达到100％的置信度。

要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。

我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计[13,16]。

设是随机变量。

则的期望值或均值定义为X X x⎰∞∞-=dx x xf x )(（8.12）这里是的概率密度分布函数。

如果从中取些独立的随机样本)(x f X )(x f ，那么的估计值就表现为个样本值的均值。

N x x x ,...,,21x N∑==Nn nxNx11ˆ（8.13）是的真正的平均值，而只是的有着准确期望值的无偏估计。

虽然x X xˆx x ˆ的期望值等于，但。

因此，我们还需要的值在附近的分布测度。

x x x≠ˆx ˆx 为了估计以及在附近的的值的分布，我们需要引入的方差，其定义X xˆx X 为与差的平方的期望值，即X x⎰∞∞--=-==dx x f x x x x x Var )(()()(222σ（8.14）由，故有2222)(x x x x x x +-=-⎰⎰⎰∞∞-∞∞-∞∞-+-=dx x f x dx x xf x dx x f x x )()(2)()(222σ（8.15）或者 222)(x x x -=σ（8.16）方差的平方根称为标准差，即 2/122)()(x x x -=σ（8.17）标准差给出了在均值附近的分布测度，并由此给出了误差幅度的阶数。