蒙特卡洛算法详讲

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 g 具有如下形式：

m c ax x g mod )()(+=

（8.1）

其中

=0x 初始值或种子（00>x ）

=a 乘法器（0≥a ）

=c 增值（0≥c ）

=m 模数

对于t 数位的二进制整数，其模数通常为t 2。例如，对于31位的计算机m 即可取1312-。这里a x ,0和c 都是整数，且具有相同的取值范围0,,x m c m a m >>>。所需的随机数序{}n x 便可由下式得

m c ax x n n m od )(1+=+ （8.2）

该序列称为线性同余序列。例如，若70===c a x 且10=m ，则该序列为

7，6，9，0，7，6，9，0…… （8.3）

可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。（8.3）式中序列周期长度为4。当然，一个有用的序列必是具有相对较长周期的序列。许多作者都用术语乘同余法和混合同余法分别指代

0=c 和0≠c 时的线性同余法。选取c a x ,,0和m 的法则可参见[6,10]。

这里我们只关心在区间)1,0(内服从均匀分布的随机数的产生。用字符U 来表示这些数字，则由式（8.2）可得

m

x U n 1

-= （8.4）

这样U 仅在数组{}m m m m /)1(,......,/2,/1,0-中取值。（对于区间（0,1）内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。其它检测方法可参见[3,6]。）产生区间),(b a 内均匀分布的随机数X ，可用下式

U 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 可得

)1l n (U X --=μ （8.8）

由于)1(U -本身就是区间(0,1)内的随机数，故可简写为

U X ln μ-= （8.9）

有时(8.6)式所需的反函数)(1x F -不存在或很难获得。这种情况可用舍去法来

处理。令dx

x dF x f )

()(=

为随机变量X 的概率密度函数。令b x a >>时的0)(=x f ，且)(x f 上界为M （即：M x f ≤)(），如图8.1所示。我们产生区间(0,1)内的两个随机数),(21U U ，则

1

1)(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 。