关于Monte－Carlo方法应用的三点注记

2013-07-25吕书龙梁飞豹刘文丽

福州大学学报（自然科学版） 2013年3期

吕书龙，梁飞豹，刘文丽

(福州大学数学与计算机科学学院，福建福州 350116)

Monte－Carlo方法是利用随机数通过计算机进行随机模拟试验的一种数值方法，该方法能够很好地应付高维数据的“维数祸根”(course dimensionality)问题.但应用该方法需要注意3个问题:①需要大量统计性质优良的均匀随机数;②计算精度只是概率精度;③收敛速度慢［1］.本文探讨这3个基本问题，尝试改进其不足，并以定积分计算为例加以说明.

1 Monte－Carlo方法计算积分

对于低维的积分计算问题，使用传统的数值积分方法基本能得到满意的结果;对于高维的多重积分问题，传统的数值积分方法存在较大局限性.Monte－Carlo方法是将积分问题转化成某个事件出现的概率或者某个随机变量或其函数的数学期望问题，然后通过随机数及计算机模拟，得到事件的频率或随机变量取值的算术平均值来近似积分结果［2］.

则积分变为:

即把求解积分I的问题转化成求解随机变量X的函数G(X)的数学期望问题.

取随机变量X的容量为n的简单随机样本Xi，i=1，2，…，n，由独立同分布大数定律有:

其方差估计式为:

从上述Monte－Carlo方法计算积分的基本原理和过程可以看出，该方法与积分问题的维数无关，其计算依赖于大量的m维随机数，且只有概率上的收敛速度O(n－1/2)和估计精度.显然要提高估计精度，就必须减少方差σ2或增大抽样容量n.这既是Monte－Carlo方法的优点也是该方法不足的根源所在.

通常取m维随机变量X服从区域D上的均匀分布，设区域D的几何测度为M(D)，则式(2)变为:

2 随机数问题

由式(2)知随机数是Monte－Carlo方法的基础.而近年来发展迅速的拟蒙特卡罗方法(Quasi Monte－Carlo)提出用“确定性的超均匀分布序列”代替一般的随机数序列进行计算，此法有更高的计算精度和更快的速度［3－5］.随机数优良的统计性质决定了Monte－Carlo方法的性能和精度.

目前大量使用的随机数都是通过递推算法产生的“伪随机数”，其统计性质主要集中在:周期性、均匀性和独立性这3个方面.因此构造超长周期、超均匀的独立随机数序列是Monte－Carlo方法的关键.

2.1 混沌映射法

An研究了混沌映射(chaotic mapping)产生随机数的方法［2］，其公式如下:

其经验分布的极限分布为:

由产生随机数的定理可知，若随机变量Y～F(y)，则x=F(Y)～U［0，1］，则通过

产生的序列{xi}可看作服从U［0，1］的随机数序列.因为序列{yi}无限不循环，所以产生的随机序列{xi}的周期是无限的.(注:文献［2］中给出的公式xi=ln(yi+1/2)/ln 3应是个笔误)

以下给出在R统计软件中实现递推算法的程序，并对混沌随机数的均匀性和独立性作直观检验.

执行:r1=rmap(500);r2=rmap(500);plot(1:500，r1);hist(r1);plot(r1，r2).得到如图1所示结果:

图1 混沌映射随机数Fig.1 Chaotic mapping random number

从图1(a)和1(b)可知，混沌随机数的均匀性相当好;由图1(c)可知，任意产生两组混沌随机数并画成点图，形成两条分段直线，可见不独立;从图1(a)和1(c)可知，该法产生的随机数的随机性和独立性相当差.

2.2 低差异序列

低差异序列(low－discrepancy sequences)即确定性的超均匀分布序列，它是拟蒙特卡罗方法的基础.如Halton序列［3－5］，其构造方法是:令b是一个素数，则任一个正整数k都可写成以b为基数的形式，

其中:di∈{0，1，2，…，b－1}，且bj－1≤k＜bj.定义基于b的根式逆函数φb(k):

显然，对于任意正整数k，φb(k)∈［0，1］，因此它可作为Halton序列的第k个元素，k=1，2，…［3－4］.当然也可有更一般的取法，即Halton序列的第k个元素可取为φb(k0+k)，k=0，1，2，…;k0∈Z+.

执行命令:plot(1:500，r_halton(500，1，2));plot(1:200，runif(500，0，1))得到图2和图3.

表1 创建Halton序列的R函数Tab.1 R －functions of halton sequence

图2 Halton序列(500个点)Fig.2 Halton sequence(500 points)

图3 均匀随机数序列(500个点)Fig.3 Uniform random sequence(500 points)

从图2可知，Halton序列的均匀性高于均匀随机数序列.与图1(a)相比，图2中的Halton序列是随机的，且不同的基数b产生的序列之间具有较强的独立性，比较结果如图4.

虽然Halton序列产生的随机数均匀性较好，独立性比混沌随机数好，但Halton序列自身存在一个相邻符号问题.记符号函数:

则Halton序列产生的相邻符号序列定义如下:

图4 两个基数不同的Halton序列点图Fig.4 Two plots of Halton sequences with different base number

这说明Halton序列本身在数值大小关系的分布上具有非参数型自相关性，对序列自身的独立性产生不利的影响，虽然两组基数不同的Halton序列独立性较好，关于这点在应用中要注意.

易知其精确值为［Φ(1)－Φ(0)］3≈0.039 772 181 953…，取3组独立的均匀、混沌、Halton随机向量序列(三维)，各序列长度均为500，随机截取5次模拟的计算结果见表2.

表2 模拟结果Tab.2 Simulation results

由表2可知，超均匀随机数序列Halton的表现最稳定，且每次模拟结果和真实值都较其它两种方法好，5次平均也最接近真实值，这正是超均匀特性所带来的好处.

3 精度和收敛速度问题

由(3)式知，要提高MC方法的精度，就要减少抽样方差或增加随机数容量.上述采用Halton序列就是减少抽样方差的一种方法;另外，增加随机数容量意味着增大计算量，但无法同速提高计算精度，因为其收敛速度阶是O(n－1/2).在提高MC方法的精度上，文献［3］提出拟蒙特卡罗法，文献［6］提出改进平均值法，特别是低维空间上，文献［1，8］提出数值法与MC法的结合，文献［7］引入数论网格法.这些都是有益的尝试和改进.文献［8］中提到的“类矩形”和“类梯形”法，将误差阶从O(n－1/2)提升到O(n－1)和O(n－2)，大大提升了收敛速度和计算精度.文献［8］采用的是均匀随机，如果采用的是低差异的随机数，结果可能会更好.基于这层考虑，下面给出拟蒙特卡罗思想与文献［8］中的“类矩形”和“类梯形”法相结合的方法.

类矩形法［8］:

类梯形法(取2n个随机数)［8］:

式(2)和式(9)的计算次数一样，但精度差别较大;而式(10)使用的随机数和计算量均是式(2)和式(9)的两倍，可是精度却未得到如文献［8］所示的跳跃式的提升，这点将在例2的计算中加以说明.

表3 均匀随机数与Halton随机数计算结果误差比较Tab.3 Error comparation of integrals by uniform random numbers and Halton random numbers

采用均匀随机数计算的结果与文献［8］给出的结果有较大的偏差，“类矩形法”并未一致地呈现出10－4级的误差阶，而是更高;同样，“类梯形”法也没有一致地表现出10－8级的误差阶，而是和“类矩形”法相差不大，可能是因为均匀随机数产生的机制不同.“类梯形”法计算虽然采用了2n个随机数，但其计算的误差阶相对于“类矩形”法并未得到较大的提升.而采用Halton随机序列的计算结果却一致地表现出比均匀随机数更高的误差阶.

另外，本文提出的采用Halton随机序列代替均匀随机数进行Monte－Carlo模拟，其计算量并没有显著增加，可计算精度却得到较大的提升，不失为一种提升Monte－Carlo方法精度和收敛速度的有效方法.

为说明模拟结果的真实性，将实现上述计算过程的R程序提供如下，以便验证.