蒙特卡罗方法.doc

蒙特卡罗方法
蒙特卡罗(Monte-Carlo，简写为M-C)方法属于计算数学的一个分支，它是在二十世纪四十年代中期为了适应当时原子能事业的发展而发展起来的，但它与一般计算方法有很大区别，一般计算方法对于解决多维或因素复杂的问题非常困难，而蒙特卡罗方法对于解决这方面的问题却比较简单。因而蒙特卡罗方法在近十年来发展很快，特别是随着快速电子计算机的发展，蒙特卡罗方法得到了迅速发展与广泛应用。蒙特卡罗方法也称随机抽样技术(Random Sampling Technique)或统计试验方法(Method of Statistical Test)。
蒙特卡罗是欧洲摩纳哥国的一个重要城市，以赌博著称。蒙特卡罗方法是以概率论与数理统计学为基础的，是通过统计试验达到计算某个量的目的。而赌博时，概率论是一种有力的手段。所以，以蒙特卡罗作为方法的名字，原因大概于此。
由于蒙特卡罗方法是利用一连串的随机数来求解问题的，因此求解随机过程，放射性衰变和布朗运动等问题，它是很有效的。它除了在原子能工业广泛应用外，在物理、化学、地质、石油、线性规划、计算机研制、计算机模拟试验、解决多体问题等领域中都有不同程度上的应用。
第一节． 蒙持卡罗方法的基本思想、特点及其局限性
一、 蒙特卡罗方法的基本思想
用下述三个例子，说明蒙特卡罗方法的基本思想。
例1．产品合格率的计算
某工厂生产一批产品，其合格率表示是：
（1.1）
为了确定合格率，应该检查这批产品的全部，确定其中合格的数目。但是，由于产品数量多，检查全部产品花费的代价大。因此，通常采取抽取部分产品，在这部分产品中确定其合格的数目。然后用这部分产品的合格率
(1.2)
来代替所要计算的合格率P。例如，检查某批产品，当被检查的产品长度介于13．60cm—13．90cm内时，则认为是合格的，否则是次品。分别抽取5件，10件，60件，150件，600件，900件，1200件，1800件来检查，其情况如下表和图20所示。

图1.1 产品抽样统计图

上表可看作八次试验，从结果中看出，随着抽取件数的增多，合格率愈来愈趋于一个稳定值o.9。如果定义随机变量w1

则做了N次试验后，正品个数共为

这样， (1.2)式可进一步写为

人们从经验中还知道，当N数目越大，r/N作为正品率的估计值就越准确的可能性也越大。类似这种把观察正品出现的频率作为近似概率的例子在生产中是很常见的。
例2．射击问题(打靶游戏)
设r表示射击运动员的弹着点到靶心的距离，g(r)表示击中r处相应的得分数(环数)，分布密度函数f(r)表示该运动员的弹着点分布，它反映运动员射击水平。积分
（1.3）
表示这个运动员的射击成绩。用概率语言说，＜g＞就是随机变量g(r)的数学期望值，即＜g＞＝Eg(r)。 现在，假设这个射击运动员射击N次，弹着点依次是r1….rz，，则自然地认为N次射击得分的平均值
（1.4）
相当好地代表了这个射击运动员的成绩。换句话说，gN是积分(1．3)式的一个估计值(或近似值)。这个例子通常称为打靶游戏，它直观地说明了蒙特卡罗方法计算定积分的基本思想。为进一步阐明这个思想，我们再举个例子：计算积分

直观上，就是在边长为1的正方形里随机投点，当点落在y＝f(x)曲线下面(见图21)，对积分值有“贡献”，否则对积分值无“贡献”。为此，设x是[0，1]上均匀分布的随机变数，定义

（1.5）
就是积分I的一个估计值。

例3．求解三维椭圆型偏微分方程的边值问题考虑三维空间中的一个体积V上的Laplace方程：
（1.6）
其中S为边界，要求f在V中的数值。
用蒙特卡罗方法求解是：将V用等距dx＝Jy＝J2的网格剖分，在扩中任一内点Q，要求f(Q)之值。用一个骰子投掷，其点数l，2，3，4，5及6设分别表示向X轴正向，负向，y轴正向，负向，2轴正向，负向移动一步。由Q点出发，按投骰子的办法，在严内进行移动，当此点到达边界时，记下边值之数值；又回头由Q点按投段子之点数移动，如此继续下去，……，设得到一系列的边值：

则f(Q)的近似值可取

或

注意，这种作法的来源是依据两个基本假设，其一为方程(1．6)可以用差分格式

来近似，其二是假设骰子为均匀的，即出现6个数机会均等。
从上述三个例子可以看到，当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。
由上面例子还可看出，用蒙特卡罗方法求解问题时，首先要建立一个随机模型，然后要制造一系列的随机数用以模拟这个过程，最后要作统计性的处理。关于建立随机模型，因问题而异。下面来研究随机数的产生方法。
二、随机数和伪随机数的产生
1．单位矩形分布
最简单、最重要、最基本的一个概率分布是单位区间(0，1)上均匀分布(或称单位矩形分布)，其分布密度函数是：

2．随机数
由具有单位矩形分布的总体内抽取的简单子样：X1，X2，…，XN，简称为随机数序列。其中的每一个体称为随机数。由于它在蒙特卡罗方法中占有极其重要的地位，我们将用专门的符号x1,x2,，…，xN来表示。 连续产生随机数的例子如投骰子，有奖储蓄开奖时的摇号码机等。
3．伪随机数及其产生的方法
计算机不会掷骰子，它是利用数论的方法来产生随机数的。由于这种办法属于半经验性质，因此只能近似地具备随机数的性质，所以称为伪随机数。最初冯·诺伊曼(Von Neumann)建议的“平方取中法”如下；以十进制数为例，平方取中法是把一个28位的十进制数自乘后，去头裁尾保留出间23个数字，然后用102S来除。例如

相应自6伪随机数序列便是0.6406，0.0368，0.1354，o.8333，0.4388，．．．。
平方取中法的一般形式是

其中“mod”表示取模运算，Xl为任意给定的初始值，由2s位十进制数组成。在电子计算机上采用二进制的数，此时平方取中法如下

X1为任意给定的2s位二进制数。 同平方取中法相类似，乘积取中法的一般形式是

X1,X2为任意两个初始值，由2s位二进制数组成。 乘同余方法的一般形式是

X1为任意初值。 乘加同余法的一般形式是

其中

关于产生伪随机数的方法，性能以及检查已有很多工作。其实，在大多数计算中，产生均匀分布的伪随机数已作为内部函数或库文件存在计算机内，随时可以调用。
三、 蒙特卡罗方法的特点及其局限性
对于多种介质，几何上不规则，维数高等类型的问题，
蒙特卡罗方法照样可用，不须修改。因此，在一定意义上说，它是一种很普遍适用的方法。而它的误差是由概率论中的大数定律

所控制的，并由中心极限定理

这表明，不等式
（1.7）
以概率1-a成立。通常a称为置信度，1-a称为置信水平。 如果s¹0，则由(1.7)式知，蒙特卡罗方法的误差Î为
（1.8）
a和Xa的关系可从正态分布N(0，1)积分表中查到。常用的几组a和Xa如下

特别称a＝0.5时的误差0.6745a／为概然误差。再如，取置信水平为95％，则Xa＝1.96，此时表明误差不等式：以95％的可能性具有精确度为E＝1.960a／。所以，M—C方法对于误差的估计具有概率性质。即对于这个方法不能断言误差不超过某值，而只能指出误差以某种(如接近1)的概率不超过某值。还可看出，当给定置信度a后，误差E由s和决定。要减小E，或者是增大N，或者是减小方差s2。在固定s下，要提高精度一位数字，就要增加100倍工作量，因此，单纯增大N，不是一个有效的办法。改进的方法之一是减少方差s。这里有许多技巧，如“重要抽样”、“系统抽样”、 “相关、“对偶变数“、“半解析方法”、“统计估计”、 “分裂与俄国轮盘赌”等等。
蒙特卡罗方法主要有如下三个特点：
1．收敛速度与问题维数无关
从误差表达式(1．8)中，明显看出，在一定的置信水平下，M—C方法的误差，除了与方差s2有关外，只取决于于样容量N，而与于样中的元素所在的集合空间**的组成无关。问题的维数变化，除了引起抽样时间及计算估计量的时间有变动外，不影响问题的误差。换言之，要达到同一精度，用M—C方法选取的点数与维数无关；计算时间仅与维数成比例。但一般数值方法，比如在计算多重积分时，达到同样的误差，点数与维数的幂次成比例，即计算量要随维数的幂次方增加。这一特性，决定了M—C方法适宜于多维问题。例如，计算积分：

用M—C方法是抽取随机数x1,0和x1,1，令

于是
（1.9）
就是积分P1的一个估计值。由于均方差，，根据误差不等式，知

于是相对误差
（1.10）
将上述积分推广到S维的情形，取

积分变成
（1.11）
此时，M—C方法的估计同样是
（1.12）
其中
（1.13）
由此可见，用M—C方法计算单重积分P1和S重积分Ps，除了确定w1值时与问题的维数有关外(前者是x1,0£x1,1，后者，仅多产生S—1个随机数) ，其它与维数无关，尤其重要的是误差与维数无关。
2．受问题的条件限制的影响小
比如计算S维中任意一个区域Ds上积分
（1.14）
无论Ds的形状如何特殊，只要能给出描述Ds特性的几何条
件，那么总可以给出类似于(1．4)式的估计如下
（1.15）
其由å是对所有满足()ÎD s的点取的。而其它数值方法受问题的条件限制影响比较大。可见， M—C方法是解决复杂的几何模型的一种有效的方法。
3．程序结构简单，占用内存单元较少
最后，在电子计算机上实现M—C计算时，程序结构清晰简单。例如，用式(1．15)计算积分(1．14)的M—C方法程序结构如下面框图所示。其中假设D s居于S维单位超立方体，DN表示每组抽样的数目，为随机数，Xa为常数，根据置信水平要求确定，E表示对计算相对误差的要求。
4．蒙特卡罗方法的局限性
前面已经指出，蒙特卡罗方法的弱点是收敛速度慢，误差的概率性质，误差与点数(抽样数N)的平方根成反比，而其它数值方法的误差是确切的，通常情况下，误差与点数成反比。因此，对于维数不高的问题(一维、二维)，数值方法可以结出误差很小的结果，而且计算量小。也就是说，M—C方法计算量相对说来比较大，又给出的是概率误差，

通常给出二至三位有效数字。除此之外，经证明，只有当系统的大小与粒子的平均自由程可以相比较时，一般在10个平均自由程左右，M—C方法算出的结果较为满意。而对于大系统、深穿透问题(一般指系统大小为16—20个平均自由程以上)，算出的结果往往偏低。而对于大系统，数值方法往往很适应，能算出较好的结果。因此，现在已有人将数值方法与蒙特卡罗方法联合起来使用，相互取长补短，克服这种局限性，取得了一定的效果。
第二节 蒙持卡罗方法计算积分
计算多重积分是蒙特卡罗方法的重要应用领域之一。为任一积分，都可看作某个随机变量的数学期望。因此，以用这个随机变量的算术平均值来近似。
一、蒙特卡罗方法求积分
例如求积分
（2.1）
选取分布密度函数

则式(2．1)可写成
（2.2）
即q是随机变量f(x)的数学期望。现从p(x)抽取随机变量X的N个样本：Xi，i＝1，2，…，N，则算术平均值
（2.3）
就是q的一个近似估计。对于一般的S重积分
（2.4）
其中，P＝P(X1，X 2，…，Xs)表示S维空间的点，VS表示积分区域。
用M—C方法求解的步骤是：
1．在所求解约区域上构造一个分布密度函数。取Vs上任一概率密度函数f(P)，它满足如下条件：

令
（2.5）
则式(2.4)可改写为
（2.6）
即q是随机变量g(P)的数学期望。
2．将g(P)的数学期望值用其算术平均值来近似。现从f(P)抽取随机向量P的N个样本：Pi，i＝1，2，…，N，则

（2.7）
就是积分q的近似估计。
选取f(P)的最简单方法是取VS上的均匀分布：
（2.8）
这里，Vs也表示积分区域的体积。这时
（2.9）
显然，存在一个如何选取f(P)的问题。选取最优的估计形式，将是M—C方法要讨论的问题之一。
二、无偏估计
如果随机变量Y，其期望值为q*，即
（2.10）
则称Y是q*的一个无偏估计。不难验证，g(P)和都是q的无偏估计。
三、 收敛性及误差估计
由上节可得如下结果：
1．当N®¥时，抽样(算术)平均以概率为1收敛到q即
（2.11）
2．如果有不为零的有限方差存在，即
（2.12）
存在，则在1-a置信水平下，成立不等式：
（2.13）
也就是说，用M—C方法求积分，其误差的阶为o(N-1/2)．
最后，给出用M—C方法产生均匀分布随机数相计算多重积分的程序。
产生均匀分布随机数的程序：
SUBROUTINE RAND(YFL)
IF(YFL．EQ．0.0) YFL＝0.314l 5926
Y＝YFL *YFL
1 Y＝Y*10.0
IF(Y-1.0) l，2，2
2 YFL＝Y-FLOAT(IFIX(Y))
RETURN
END
程序功能
本程序用于产生(0，1)区间上均匀分布的随机数。
使用说明
1．子程序语句
SUBROUTINE RAND(YFL)
2．哑元说明
YFL：实变量，输入、输出参数，要求在第一次调用前，对 YFL赋值， 当YFL＝0.0时，程序自带初值：YFL＝0.31415926；当使用者自己选定初值时，可将选定的初值赋于YFL。以后每调用一次便在YFL中产生一个(o，1)区间上的随机数，同时，YFL又作为产生下一个随机数的初值。
方法简介
本程序采用的计算步骤如下：
1．对于结定的一个P位数的初值r。(例如取ro＝0.3415926，0.27l 828l 8，…)，对ro进行平方，得到一个2P位数，实际上由于计算机字长有限，尾数长度是有限制的，当的位数超出计算机尾数长度的时候，机器自动会舍掉多出的尾数，这样截去尾数避免了任意数的平方的末位数只能出现o，1，4，5，6，9而不会出现2，3，7，8的系统偏倚性。
2．判断是否小于1，如果＜1时，则右移的小数点直至³1时为止，然后，截去整数部分，保留小数部分，把保留的小数部分作为(0，1)区间上的均匀分布的随机数r1，这样截去首位避免了小于1的数的平方有对小数目偏倚的现象。
3．以r1为初值，重复l，2过程产生r2。这一过程一直重复下去，就可产生(o，1)区间上均匀分布的随机数序列r1，r2， ....。
例题
例1．取初值YFL＝0.31415926所产生的2000个(0，1）区间上的均匀随机数进行统计特性检验。
1)参数检验
最大值max(rl，r2，…，r2000)= 0.99907696
最小值min(r1，r2，…，r2000)= .17911196×10-3
一阶矩 ＝0.48966279
二阶矩 ＝0.31479022
二阶中心矩 ＝0.84561238
2)均匀性检验
将(0，1)区间等分成39个小区间，对2000个随机数进行均匀性检验，由计算得：
X2(38)＝45.316(45.316＜53.33354)
在显著水平a＝0.05下，接受均匀性假设。
例2．取YFL＝0.27182818产生2000个(0，1)区间上均勾分布的随机数。
计算结果
最大值 max(rl，r2，…，r2000,)＝0.999828589
最小值 min(rl，r2，…，r2000,)＝0.57747960×10-3，
一阶矩 ＝0.49808762．
附注
1．设rl为(o，1）区间上的均匀分布的随机数，则
R1=a+(b-a)r1，
即为(a，b)区间上均匀分布的随机数。
2．本方法程序简单，只要一个存贮单元就够了。但在使用上要注意初值的选取，初值选得好可以产生周期较长、统计特性较好的随机数。
产生均匀分布随机数的方法很多，这里就不一一列举了。下面给出用M—C方法计算多重积分的程序：
SUBROUTINE RJ (N， N1， YFL， RY)
F＝0.0
DO 5 I＝1， N1
5 F＝F十FX(N)
RY＝F/FLOAT(N1)
RETURN
END ·
SUBROUTINE RAND(YFL)
… 随机数程序段体部分
END
程序功能
本程序采用蒙特卡罗(M—C)方法求多重积分

的近似值。
使用说明
1.子程序语句
SUBROUTINE RJ(N,N1,YFL,RY)
2.哑元说明
N：整变量，输入参数，积分重数。
N1：整变量，输入参数，一个充分大的正整数。
YFL： 同子程序RAND的哑元说明。
RY：实变量，输出参数，存放积分结果。
3．所调用的过程
SUBROUTINE RAND(YFL)， 表示产生(0，1) 区 间上均匀分布的随机数的子程序段。
FUNCTION FX(N)，表示计算被积函数的函数段，由使用者自行编写，其自变量由于程序RAND产生的随机数求得，N为积分重数。
例题
取S＝2，3，…，8计算多重积分

的结果。
精确结果

当产生随机数的初值取YFL＝0.31415926时，计算结果如下表：

可见，计算结果的误差小于10％。
11