完整版蒙特卡洛算法详讲.docx

完整版蒙特卡洛算法详讲.docx

《完整版蒙特卡洛算法详讲.docx》由会员分享，可在线阅读，更多相关《完整版蒙特卡洛算法详讲.docx（18页珍藏版）》请在冰豆网上搜索。

完整版蒙特卡洛算法详讲

MonteCarlo法

§8.1概述

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

MonteCarlo方法（MCM），也称为统计试验方法，是理论物理学两大主要学科的合并：

即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态[1]。

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

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

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

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

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

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

MonteCarlo方法的应用有两种途径：

仿真和取样。

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

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

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

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

这就是数值积分的MonteCarlo方法。

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

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

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

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

§8.2随机数和随机变量的产生

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

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

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

最一般的方程具有如下形式：

（8.1）

其中

初始值或种子（）

乘法器（）

增值（）

模数

对于数位的二进制整数，其模数通常为。

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

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

所需的随机数序便可由下式得

（8.2）

该序列称为线性同余序列。

例如，若且，则该序列为

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

可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。

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

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

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

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

这里我们只关心在区间内服从均匀分布的随机数的产生。

用字符来表示这些数字，则由式（8.2）可得

（8.4）

这样仅在数组中取值。

（对于区间（0,1）内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。

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

）产生区间内均匀分布的随机数，可用下式

（8.5）

用计算机编码产生的随机数（利用式（8.2）和（8.4））并不是完全随机的；事实上，给定序列种子，序列的所有数字都是完全可预测的。

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

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

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

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

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

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

为的概率）。

显然表明，通过产生（0,1）内均匀分布随机数，经转换我们可得服从分布的随机样本。

为了得到这样的具有概率分布的随机数，不妨设，即可得

（8.6）

其中具有分布函数。

例如，若是均值为呈指数分布的随机变量，且

（8.7）

在中解出可得

（8.8）

由于本身就是区间（0,1）内的随机数，故可简写为

（8.9）

有时（8.6）式所需的反函数不存在或很难获得。

这种情况可用舍去法来处理。

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

令时的，且上界为（即：

），如图8.1所示。

我们产生区间（0,1）内的两个随机数，则

（8.10）

分别为在（a,b）和（0,M）内均匀分布的随机数。

若

（8.11）

则为的可选值，否则被舍去，然后再试新的一组。

如此运用舍去法，所有位于以上的点都被舍去，而位于上或以下的点都由来产生。

图8.1舍去法产生概率密度函数为的随机变量

例8.1设计一子程序使之产生0,1之间呈均匀分布的随机数。

用该程序产生随机变，其概率分布由下式给定

解：

生成的子程序如图8.2所示。

该子程序中，且。

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

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

0001C**********************************************************

0002CPROGRAMFORGENERATINGRANDOMVARIABLESWITH

0003CAGIVENPROBABILITYDISTRIBUTION

0004C**********************************************************

0005

0006DOUBLEPRECISIONISEED

0007

0008ISEED=1234.D0

0009DO10I=1,100

0010CALLRANSOM（ISEED,R）

0011THETA=ACOSD（1.0-2.0*R）

0012WRITE（6,*）I,THETA

001310CONTINUE

0014STOP

0015END

0001C**********************************************************

0002CSUBROUTINEFORGENERATINGRANDOMNUMBERSIN

0003CTHEINTERVAL（0,1）

0004C**********************************************************

0005

0006SUBROUTINERANDOM（ISEED,R）

0007DOUBLEPRECISIONISDDE,DEL,A

0008DATADEL,A/2147483647.D0,16807.D0/

0009

0010ISDDE=DMOD（A*ISDDE,DEL）

0011R=ISDDE/DEL

0012RETURN

0013END

图8.2例8.1的随机数生成器

图8.2的子程序只是为了说明本章所介绍的一些概念。

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

为了生成随机变量，令

则有

据此，一系列具有给定分布的随机变量便可由图8.2所示主程序中生成。

§8.3误差计算

MonteCarlo程序给出的解按大量的检测统计都达到了平均值。

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

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

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

设是随机变量。

则的期望值或均值定义为

（8.12）

这里是的概率密度分布函数。

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

（8.13）

是的真正的平均值，而只是的有着准确期望值的无偏估计。

虽然的期望值等于，但。

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

为了估计以及在附近的的值的分布，我们需要引入的方差，其定义为与差的平方的期望值，即

（8.14）

由，故有

（8.15）

或者（8.16）

方差的平方根称为标准差，即（8.17）

标准差给出了在均值附近的分布测度，并由此给出了误差幅度的阶数。

的标准差与的标准差的关系表示为

（8.18）

该式表明，如果用根据（8.13）式由个值构成的来估计，那么结果中在附近的扩散范围便与成比例，且随着样本数的增加而降低。

为了估计的扩散范围，我们定义样本方差为

（8.19）

由此式还可看出的期望值等于。

因此样本方差是的无偏估计。

将（8.19）式中平方项乘出来，便可得样本标准差为

（8.20）

当较大时，系数可设为1。

作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数

（8.21）

该式表明次独立随机试验中有次成功的概率。

（8.21）中，是一次试验中成功的概率，且。

当和都较大时，便可用Stirling公式

（8.22）

因此（8.21）式可近似为正态分布[17]

（8.23）

其中且。

因此当时，中心极限定理表明，描述由点MonteCarlo算法获得的的分布的概率密度函数是（8.23）式所示的正态分布函数。

也就是说，大量随机变量的集合趋于呈正态分布。

将（8.18）式代入（8.23）式可得

（8.24）

正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。

高斯模型的显著特性源于中心极限定理。

因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。

例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。

由于样本数是有限的，所以MonteCarlo算法不可能达到绝对的确定性。

故而在附近估计出某一范围或区间以确保我们估计的落入该范围内。

假设要得到位于和之间的概率。

由定义

（8.25）

令可得

（8.26a）

或

（8.26b）

其中是误差函数，且是标准正态差的上分位点。

对（8.26）式可做如下解释：

如果用来呈现独立随机观测值并构建相关随机区间的MonteCarlo程序以较大的值反复运行，则这些随机区间的分位点将近似等于。

随机区间称为置信区间，称为置信度。

大多数的MonteCarlo算法都使用误差，这表明是在实际均值的标准差范围内的。

由（8.26）式可得，样本均值位于区间内的概率是0.6826或68.3％。

若要求更高的置信度，可用两个或三个标准差的值。

如

（8.27）

其中是标准差的个数。

在（8.26）和（8.27）式中均假设总体标准差为已知。

由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差来估计的值，从而使学生氏－分布取代正态分布。

我们知道当很大，比如时，－分布近似趋于正态分布。

（8.26）式等价于

（8.28）

其中为自由度是的学生氏－分布的上分位点；其值在任何标准统计学课本中都有表列可查。

这样置信区间的上、下限便可由下式给出

上限＝（8.29）

下限＝（8.30）

MonteCarlo算法中关于误差估计的进一步讨论，可参见[18,19]。

例8.2利用中心极限定理产生一高斯（正态）分布的随机变量。

根据中心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总