蒙特卡罗计算多重积分

1、5第一作者姓名：文章标题关于计算多重积分的拟蒙特卡罗方法张 涛，周仲礼，安素珍，张 军(成都理工大学数学地质四川省重点实验室，四川成都 610059)摘 要：介绍了拟蒙特卡罗方法计算多重积分的基本原理，对halton序列和rand函数产生的序列的均匀性进行了比较，并给出了拟蒙特卡罗方法计算多重积分的步骤、实例、和matlab语言编写的计算程序实例充分体现了采用拟蒙特卡罗方法计算多重积分的有效性、精确性和优越性关键词：拟蒙特卡罗方法；halton序列；多重积分中图分类号：o242 文献标志码：a 文章编号：1674-3563(2012)01-0001-06doi：10.3875/j.issn.1

2、674-3563.2012.01.001 本文的pdf文件可以从获得在数学和工程科学计算中，求解多重积分的近似值是一个至关重要的环节通常人们所采用的方法有3种，即蒙特卡罗方法、拟蒙特卡罗方法和数论网格方法然而，在实际应用中，数论网格法很难解决高维问题，如对20维这样的问题，就要有一百万多个点1，计算量大体上随维数的幂次增加，几乎是不可能计算的蒙特卡罗方法是采用单和来得到多重积分的近似值的2，计算跟问题的维数无关，克服了“维数灾难”，但其误差却是概率性的，得到解的精度不高本文研究了拟蒙特卡罗方法在计算多重积分中的应用，并给出matlab程序的实现实例表明，采用拟蒙特卡罗方法计算多重积分，既能克服

3、维数的制约，又能得到确定性的误差估计和更高精度的解1 拟蒙特卡罗方法基本原理及其优越性拟蒙特卡罗方法也称低差异序列法，是在蒙特卡罗方法基础上发展起来的一种模拟方法拟蒙特卡罗方法与蒙特卡罗方法相似，但所用的理论基础却不同拟蒙特卡罗方法是通过构造所谓的低差异序列，即用的是确定性的点，然后用koksma-hlawka不等式来确定误差阶的，而不是根据大数定律3对于单位超立方体上的积分，拟蒙特卡罗估计的形式和蒙特卡罗类似，即：， （1）只是此处的点列为确定性的点因此，拟蒙特卡罗方法的基本思想是用精选的确定性点来代替蒙特卡罗方法中的随机点拟蒙特卡罗方法积分的误差估计是根据koksma-hlawka不等式来

4、给出，即令f为区域上hardy krause意义变分有界的实函数，对任意n个点，拟蒙特卡罗方法模拟解积分问题的误差界可以表示为3-4：， （2）含n个点的拟随机序列的偏差为，所以拟蒙特卡罗积分的误差阶5是，而蒙特卡罗积分的误差的阶是对于拟蒙特卡罗方法积分，我们有确定性的误差估计在计算积分中还可以明智地选择这些确定的点来减小式（1）中的误差，得到高精度的解因此，就确定性和高精度性这两点，拟蒙特卡罗积分要优于蒙特卡罗积分对于更高维的积分问题，拟蒙特卡罗法优于数论网格法在过去的几十年里，拟蒙特卡罗法得到了快速发展，人们已经构造出大量优质的低差异序列，如halton序列、sobol'序列、ni

5、derreiter的（t,m,s）网格和（t,s）序列等6本文主要研究halton序列在计算多重积分中的应用2 halton序列halton序列是van der corput序列的推广halton序列是通过将一系列整数表示成某个基的数位的形式，然后将这些数位按反序排列再在前面加小数点而得到的值本文将s维halton序列表示成，其中每一个随机数是一个维向量，即生成halton序列的简单易行的步骤如下：首先选择个基比如选择前个素数，然后对某个整数，将表示成以为基的数位，将这些数位按反序排列再在前面加小数点得到新的值，便是序列中某个随机数向量的第个元素即对某个整数有以下步骤：1）选择适当大的，对每个

6、基将表示成6：，；2）将数位按反序排列再在前面加小数点得到新的值，；3）置，然后重复1）和2）两步举个例子，比如3，即维数是3，设15，选择2,3和5作为基，有1511112，151203，15305，于是便有第一个随机变量(0.11112, 0.0213, 0.035)，即(0.937500, 0.259259, 0.120000)将加1，如此循环，便得到整个序列下面将给出产生halton序列的matlab程序7，并对halton序列与rand函数产生的随机数序列的均匀性进行比较，见图1function seq=gethalton(n,base)seq = zeros(n,1);numbit

7、s = 1+ceil(log(n)/log(base);vetbase = base.(-(1:numbits);workvet = zeros(1,numbits);for i=1:nj=1;ok = 0;while ok = 0workvet(j) = workvet(j)+1;if workvet(j)<baseok = 1;elseworkvet(j) = 0;j = j+1;endendseq(i) = dot(workvet,vetbase);图1 halton序列与rand函数产生的十维随机数点的均匀性比较end由图1可明显看出，rand函数产生的伪随机数序列有抱团现象，分

8、布不均匀，而halton序列则要均匀得多有了这样优质的均匀随机序列，用拟蒙特卡罗方法求解得到的将是确定性的误差，从而避免了蒙特卡罗方法得到概率误差的缺陷3拟蒙特卡罗方法计算多重积分的均匀随机数平均值法考虑重积分：，其中为维积分域拟蒙特卡罗方法求解步骤如下8：1）在所求积分区域上构造一个分布密度函数取上任一概率密度函数使它满足当时，令，则=，即是随机变量函数的数学期望2）用算术平均值来近似的数学期望抽选服从的个样本点，则算术平均值就是积分值的一个近似估计值选取最简单的方法是取上的均匀分布，即，这里也表示积分区域的体积4 算例及其实现下面给出两个特殊数值算例验证拟蒙特卡罗方法在计算多重积分中的有效

9、性和优越性例1 用halton序列计算重积分的估计值取2000，在计算机上分别用halton序列和rand函数计算的结果见表1表1 分别用halton序列和rand函数计算的积分结果sg(h)g(r)3-0.566 6-0.567 86-0.549 6-0.673 68-0.566 2-0.585 310-0.571 3-0.590 9本例中仅给出3时采用halton序列的matlab程序：3;n=2000;a(:,1)=gethalton(n,2);a(:,2)=gethalton(n,3);a(:,3)=gethalton(n,5);for i=1:2000a_max(i)=max(a(i

10、,:);endb=s*log(a_max); c_n=abs(b);d_n=log(c_n);g_n=mean(d_n).例2 用halton序列计算多重积分的估计值9取1000，0，在计算机上分别用halton序列和rand函数计算的结果见表2表2 分别用halton序列和rand函数计算的结果sg(h)g(r)40.989 50.911 090.893 01.053 3131.116 11.409 9181.286 90.720 3本例中仅给出4采用halton序列的matlab程序：4;n=1000;a(:,1)=gethalton(n,2);a(:,2)=gethalton(n,3);

11、a(:,3)=gethalton(n,5);a(:,4)=gethalton(n,7);jf_sum=0;for i=1:nlc=1;for j=1:slc=lc*abs(4*a(i,j)-2);endjf_sum=jf_sum+lc;endjf=jf_sum/n.5 结 语通过计算，可以看出拟蒙特卡罗方法是有效的，积分重数对积分结果的误差无明显影响，精度较高，计算机程序实现也是很简单的本文只是用了原始的halton序列，由于halton序列基数越大，相关性也变大，均匀性也会越差，这将会影响结果的精度，因此，如何选择更优质的低差异序列有待进一步学习和研究参考文献1 雷桂圆. 关于蒙特卡罗及拟蒙

12、特卡罗方法的若干研究d. 浙江: 浙江大学理学院, 2003: 3-5.2 杜绍洪. 高维积分的新型求积公式d.四川：四川大学数学学院，2004：11-15. 3 morokoff w j, caflisch r e. quasi-random sequences and their discrepancies j. siam j sci comput, 1994, 15(6): 1251-1279.4 niederreiter h. random number generation and quasi monte carlo methods m. philadelphia: siam, 19

13、92:,pages 10-21.5 hozumi morohosi,masanori fushimi. a practical approach to the error estimation of quasi-monte carlo integrationsm. tokyo:mathematical engineering ,1998:1-10.6 inna krykova. evaluting of path-dependent securities with low discrepancy methodsd. massachusetts:worcester polytechnic ins

14、titute, december 2003, pages 7-30.7 paolo brandimarte . numerical method in finance:a matlab-based introductionm. 2002 by john wiley &sons.inc.new york. pages 265-280.8 宫野. 计算多重积分的蒙特卡罗方法与数论网格法. 大连理工大学学报j. 2001.1. 第41卷第1期，21-22.9 xiaoqun wang,kai-tai fang. the effective dimension and quasi-monte

15、carlo integrationj. journal of complexity volume 19 ,issue 2,april 2003,pages 101-124.quasi monte carlo method for calculating multiple integralszhang tao, zhou zhongli, an suzhen, zhang jun(key laboratory of mathematical geology of sichuan province, chengdu university of technology, chengdu, china

16、610059)abstract: basic principles of quasi monte-carlo method for calculating multiple integrals were introduced. the uniformities of halton sequences and of sequences produced by the rand function were compared then. later, steps, examples, and computer program written in matlab language of quasi monte- carlo method