关于计算多重积分的拟蒙特卡罗方法
第7讲多重积分(2007)(蒙特卡洛)
毫无困难地推广到多重积分，这正是Monte
Carlo方法的优点。
计算积分
f (x1,...,xs )dx1 dxs
1.产生给定分布随机变量的方法 如果f(x1,…,xs)是s维随机变量(X1,…,Xs)的概
率密度，则上述积分等价于
f
*(x1i ,...,
xsi )
f (x1i ,..., g(x1i ,...,
xsi ) xsi )
体积
f
(x1i ,...,
xsi), i
1,...,
n
积分的估计值为
1 n
n
i1
f
*(x1i ,...,
xsi )
体积 n
n
i1
f
(x1i ,...,
xsi )
值得注意的是，用随机点法和平均值法计算多 重积分时，误差分析与单重积分类似，误差的 阶依然为O(n -1/2)，它与积分重数无关。
xs )dx1
dxs
Ef
*(X1,...,Xs )
于是我们可以根据s维随机变量(X1,…,Xs)的 概率密度g(x1,…,xs) ，抽取n个随机样本
(x11,..., xs1),..., (x1n,..., xsn)
计算 f *(x11,..., xs1),..., f *(x1n,..., xsn), 然后用它们的算术
弹着点的纵向误差. (X, Y)服从二维正态分
布, 联合概率密度为
f (x, y)
1
1 ( x2 y2 )
e2
2
求任意发射一枚弹，弹着点 (X, Y)落入单位园内
的概率。
2.简单随机数方法
关于计算多重积分的拟蒙特卡罗方法
Ｍｅｔｈｏｄｓ［ＤＩ．Ｍａｓｓａｃｈｕｓｅｔｔｓ：Ｗｏｒｃｅｓｔｅｒ
Ｐｏｌｙｔｅｃｈｕｉｃ Ｉｎｓｔｉｔｕｔｅ，２００３：７－３０．
［７］Ｂｒａｎｄｉｍａｒｔｅ Ｐａｏｌｏ．Ｎｕｍｅｒｉｃａｌ Ｍｅｔｈｏｄ ｉｎ Ｆｉｎａｎｃｅ：Ａ ＭＡＴＬＡＢ—ｂａｓｅｄ Ｉｎｔｒｏｄｕｃｔｉｏｎ【Ｍ］．Ｎｅｗ Ｙｏｒｋ：Ｊｏｈｎ Ｗｉｌｅｙ ａｎｄ
的估计值．
表１
２１５ ６６４ ９
取Ｎ＝２ ０００，在计算机上分别用Ｈａｌｔｏｎ序列和Ｒａｎｄ函数计算的结果见．表１．
分别用Ｈａｌｔｏｎ序列和ｒａｎｄ函数计算的积分结果
本例中仅给出Ｓ＝３时采用Ｈａｌｔｏｎ序列的Ｍａｔｌａｂ程序：
Ｓ＝３；Ｎ＝２０００；
ａ（：，１）２ＧｅｔＨａｌｔｏｎ（Ｎ，２）；ａ（：，２）＝ＧｅｔＨａｌｔｏｎ（Ｎ，３）；ａ（：，３）＝ＧｅｔＨａｌｔｏｎ（Ｎ，５）；
万方数据
３６
＋
（ｆｏ≠）曩，…而如，删禁ｆ、
温州大学学报・自然科学版（２０１２）第３３卷第５期 ”ｐ”
ｘｌ，Ｘ２，…，≮）＝｛ｇ（ｘ１，工２，…，墨）’
ｌ
０，
”。∥。，
ｇ（ｘ１，Ｘ２，…，ｘｓ）＝０
则
，＝』．．．ｐ（誓，Ｘ２，…，ｘ）ｄｘ，…咄＝『Ｉ．．矽（五，Ｘ２，…，≮）ｇ（五，ｘ２，…，ｔ）咄…咄：
数越大，相关性也越大，均匀性也会越差，这将会影响结果的精度，因此，如何选择更优质的低 差异序列有待进一步学习和研究．
参考文献 ［１］雷桂圆．关于蒙特卡罗及拟蒙特卡罗方法的若干研究［Ｄ】．杭州：浙江人学理学院，２００３：３－５． ［２］杜绍洪．高维积分的新型求积公式［Ｄ】．成都：四川大学数学学院，２００４：１１—１５． ［３］ＭｏｒｏｋｏｆｆＷ
计算多重广义积分的一种有效方法和程序
( - LOG( Z( K&) ) )
© 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved.
PRINT ERR ,ERL END IF
DEFDBL A - Z CONST Pai = 3. 141592653589793 # Integrand = EXP ( - ( X (1) ∧2 + X(2) ∧2) Π 2) Π 2
宁德师专学报 ( 自然科学版) 1998 年 12 月 ・270 ・
1 n ζ ζ ζ 定理 1 (柯尔莫哥洛夫) 设 ζ 是独立同分布随机变量序列 ,lim ∑ 1 , 2 , …, n , … i = a
[3] + ∞ - ∞
…∫
+ ∞ - ∞
h ( x1 , x2 , …, x n ) g ( x1 ) g ( x2 ) …g ( x n ) d x1 … d x n . 设 X1 , X2 , …,
+ ∞
Xn 是独立同分布于 I (0 ,1) 的随机变量 ,记 Y = h ( X1 , X2 , …, Xn ) ,则 E ( Y) = ∫-
∞
…∫-
+ ∞
∞
h ( x1 ,
…, x n ) g ( x1 ) …g ( x n ) d x1 … d x n = I . 根据以下两个定理 ,很容易得到 I 的近似值 .
收稿日期 :1997 - 11 - 30 李笃群 ,男 ,1974 年 9 月出生 ,研究生
Ξ
© 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved.
Monte Carlo积分
Sample Mean Method
Start Set N : large integer
s1 = 0, s2 = 0
Loop N times
xn = (b-a) un + a yn = (xn)
s1 = s1 + yn , s2 = s2 + yn2 Estimate mean m’=s1/N Estimate variance V’ = s2/N – m’2 End
Monte Carlo method WINS, when d >> 3
第九章 Monte Carlo积分
Hit-or-Miss Method Sample Mean Method Variance Reduction Technique Variance Reduction using Rejection Technique
E ( M ) Np 2 ( M ) Np(1 P)
M ˆ ˆ E ( I ) E ph(b a) E h(b a) N h(b a) E ( M ) h(b a) p I N
ˆ ˆ 2 ( I ) 2 ph(b a) 2
Hit-or-Miss Method
I p : h(b - a)
So, if we can estimate p, we can estimate I:
ˆ ˆ I ph(b - a)
ˆ where p is our estimate of p
Hit-or-Miss Method
We can easily estimate p:
write this as
1 I 3 e dx 3E[e X ] 3 0
数值分析21求积分的蒙特卡罗方法
) I 1(
h 2
) 2(
2
h
) k(
4
h
)
2k
2k 2
所以
[ 4T (
h 2
) T ( h )] / 3 I O ( h )
4
记
4T ( ) T ( h ) 2 T1 ( h ) 41
4 6 2k
h
T1 ( h ) I 2 h 3 h k h
m
0
1
2
3
4
5
h 3
[ f ( x 2 k 2 ) 4 f ( x 2 k 1 ) f ( x 2 k )]
k 1
复合Simpson公式
Sm
S1
h 3
h1 3
[ f (a ) f (b ) 2 f ( x 2 k ) 4 f ( x 2 k 1 )]
3/16
R[ f ]
f
(4)
( )
4!
x2 x0
( x x 0 )( x x 1 ) ( x x 2 ) dx
2
令 h =(b – a)/2, x = x0+ t h ,则
x2 x0
( x x 0 )( x x 1 ) ( x x 2 ) dx h
2
5
2 0
《数值分析》 21
Simpson公式的误差
格林公式中曲线积分处理 右矩形公式应用 求积分的蒙特卡罗方法
龙贝格外推计算公式
I[ f ]
S[ f ]
b
R[ f ]=I[ f ] – S[ f ]
f ( x ) dx
MonteCarlo方法及其简单应用（图文）
MonteCarlo方法及其简单应用(图文)论文导读：本文介绍了Monte Carlo方法的思想，主要从在定积分计算方面介绍了随机投点法和平均值法，并将其推广到二重积分、三重积分和多重积分情形，最后以棋手分奖金问题介绍了Monte Carlo方法在古典概率问题中的应用.分析了误差，介绍了减少误差的方法. 给出这些方法的实例及其Mathematica实现程序.关键词：MonteCarlo方法，积分计算，古典概率，模拟1 引言Monte Carlo方法，源于第二次世界大战美国关于研制原子弹的“曼哈顿计划”.该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城——摩纳哥的Monte Carlo——来命名这种方法,为它蒙上了一层神秘色彩.Monte Carlo方法的基本思想很早以前就被人们所发现和利用.19世纪人们用投针试验的方法来确定圆周率.20世纪40年代电子计算机的出现,特别是近年来高速电子计算机的出现,使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能.Monte Carlo方法研究的问题大致可分为两种类型:一种是问题本身就是随机的,另一种本身属于确定性问题,但可以建立它的解与特定随机变量或随机过程的数字特征或分布函数之间的联系,因而也可用随机模拟方法解决.文[1]-[7] 介绍了Monte Carlo方法的思想，但没有给出具体的实例及实现过程。
发表论文。
本文介绍了MonteCarlo方法的思想，从计算定积分和古典概率两方面的应用进行研究，给出了实例及其Mathematica实现程序.2 Monte Carlo方法2.1 Monte Carlo方法思想概述Monte Carlo方法，有时也称随机模拟(RandomSimulation)方法或统计试验(Statistical Testing)方法.它的基本思想是:首先建立一个概率模型或随机过程,使它的参数等于问题的解；然后通过对模型或过程的观察、抽样来计算所求参数的统计特征；最后给出所求解的近似值,而解的精度可用估计值的标准误差来表示.假设所求的量是随机变量的数学期望,那么近似确定的方法是对进行重复抽样,产生相互独立的值的序列并计算其算术平均值:根据大数定理，当充分大时,以概率1成立，即可用作为的估计值.Monte Carlo方法以概率统计理论为基础,以随机抽样(随机变量的抽样)为手段，在很多方面有重要的应用.它的优点表现在三个方面：方法和程序的结构简单，易分析、易理解；收敛的概率性和收敛速度与问题的维数无关，很好的避免了维数问题；受问题条件限制的影响较小，很好的提高可行性.使用Monte Carlo方法的步骤如下：(l)构造或描述概率过程(2)实现从已知概率分布中抽样(3)建立各种估计量2.2 Monte Carlo方法的可行性从Monte Carlo方法的基本思想可以得到它通常的做法，利用数学或物理方法产生[0，1]中均匀分布的随机数，在变换得到任意分布的随机数.随机数个数很大时，可以由大数定理，求出事件的概率值.这种做法的可行性主要依据下面的事实：（1）如果随机变量的分布函数是，由于非降.对于任意的,(),可以定义：作为的反函数.我们考虑随机变量的分布，这里假定是连续函数，则对于有：（1）即服从上的均匀分布.（2）反之，如果服从上的均匀分布，则对于任意的分布函数，令，则：（2）因此是服从分布函数的随机变量.所以我们只要能够产生中均匀分布的随机变量的子样，那么通过（2）式我们就可以得到任意分布函数的随机变量的子样.再结合大数定理、就可以运用Monte Carlo方法进行随机模拟，解决一些实际的问题.3 Monte Carlo方法在定积分中的应用3.1随机投点法对于定积分.为使计算机模拟简单起见，设，有限，，令，并设是在上均匀分布的二维随机变量，其联合密度函数为.则是中曲线下方的面积（如图2）.图2假设我们向中进行随机投点.若点落在下方（即）称为中的，否则称不中.则点中的概率为，若我们进行次投点，其中次中的.则可以得到的一个估计（3）该方法的具体计算步骤为：①独立地产生2个随机数,,i=1,…,n；②计算，，和；③统计的个数；④用（3）估计.例1 1777年，法国学者Buffon提出用试验方法求圆周率的值.原理如下：假设平面上有无数条距离为1的等距平行线，现向该平面随机地投掷一根长度为的针.则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解：针中心与最近的平行线间的距离x均匀地分布在区间上，针与平行线间的夹角（不管相交与否）均匀地分布在区间上（如图1）.于是，针与线相交的充要条件是，从而针线相交概率为：图1而由大数定律可以估计出针线相交的概率，其中为掷针次数，为针线相交次数，从而圆周率.其mathematica实现语句见附录1.3．2 样本平均值法对积分，设是上的一个密度函数，改写（4）由矩法，若有个来自的观测值，则可给出的一个矩估计，这便是样本平均值法的基本原理.若，有限，可取.设是来自的随机数，则的一个估计为（5）该方法的具体计算步骤为：①独立地产生个随机数；②计算和，；③用（5）估计.后面将给出一个例子说明此方法的应用.4 Monte Carlo方法在计算多重积分中的应用方法一：（重积分）（7）其中为S维单位立方体，，在上有：.很明显.此时积分（5）可以看作为求维空间长方体V：的体积.即：（8）对于这种较为一般形式的多重积分计算问题，采用的还是随机投点.具体步骤如下：首先产生个随机数（i=1,2,…，）及，构造维随机向量，然后检验是否落后在V中，同理可以推论.检验是否成立，如果在构成的个随机向量中，有个随机向量落于V中，那么取作为积分的近似值，即，如果积分区域及被积函数不满足上述条件，那么可以通过变换便可达到所希望的条件.方法二：其中积分区域包含在维多面体中，此多面体决定于个不等式.设函数在内连续且满足条件：，是在维多面体中均匀分布的随机质点的个数，是在个随机点之中落入以维区域V为底以为顶之曲顶柱体内的随机点的个数.这里表示由不等式和决定的维多面体.则重积分的Monte Carlo近似计算公式为：=（9）例 2 在三维空间中，由三个圆柱面：，，围成一个立体，利用Monte Carlo方法求它的体积.分析：据题意，所求体积,其中{，，且，，}.记，，},考虑在空间内随机的产生个点，落在空间内有个，则.在Mathematica中模拟程序见附录2.5 在古典概率问题中的应用下面的例子说明了Monte Carlo方法在古典概率中的应用.例3 甲乙两位棋手棋艺相当，现他们在一项奖金为1000元的比赛中相遇，比赛为五局三胜制，已经进行了三局的比赛，结果为甲三胜一负，现因故要停止比赛，问应该如何分配这1000元比赛奖金才算公平？分析：平均分对甲欠公平，全归甲则对乙欠公平.合理的分法是按一定的比例分配.现在我们用计算机模拟两位棋手后面的比赛，是否就可以知道奖金分配方案.由于两位棋手的棋艺相当，可以假定他们在以下每局的比赛胜负的机会各半.Mathematica中函数产生随机数0或1，0与1出现的机会各占一半，可以用随机数1表示甲棋手胜，而随机数0表示乙胜.（也可以用中的随机实数来模拟两人的胜负，随机数大于0.5表示甲胜，否则乙胜）连续模拟1000次（或更多次数）每次模拟到甲乙两方乙有一方胜了三局为止.按所说方案分配奖金，1000次模拟结束后，计算两棋手每次的平均奖金，就是该棋手应得的奖金.模拟结果：甲：750，乙：250（程序见附录1）最终以甲分到；乙分到.即甲750元，乙250元.实际上，因为比赛只需进行两局.则可分出胜负.结果无非是以下四种情况之一：甲甲、甲乙、乙甲、乙乙.上面四种情况可看出，甲获胜的概率为，乙获胜的概率为.在Mathematica 中模拟程序见附录3.6 误差分析6.1 收敛性蒙特卡罗方法是由随机变量的简单子样的算术平均值：作为所求解的近似值.由大数定律可知，如独立同分布，且具有有限期望值（<∞），则.即随机变量的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值.6.2 误差蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案.该定理指出，如果随机变量序列，，…，独立同分布，且具有有限非零的方差，即是的分布密度函数.则当N充分大时，有如下的近似式其中称为置信度，1－称为置信水平.这表明，不等式近似地以概率1－成立，且误差收敛速度的阶为.通常，Monte Carlo方法的误差ε定义为上式中与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出.关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的.第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出.例4 求用平均值法估计圆周率，并考虑置信度为5%，精度要求为0.01的情况下所需的试验次数.解：易知，故考虑令～，令，其期望值为，因此=，其中是[0，1]区间上均匀分布的随机数.此时，，，，所以（次）.6.3 减小方差的各种技巧显然，当给定置信度α后，误差ε由σ和N决定.要减小ε，或者是增大N，或者是减小方差.在固定的情况下，要把精度提高一个数量级，试验次数N 需增加两个数量级.因此，单纯增大N不是一个有效的办法.另一方面，如能减小估计的均方差σ，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果.因此降低方差的各种技巧，引起了人们的普遍注意.一般来说，降低方差的技巧，往往会使观察一个子样的时间增加.在固定时间内，使观察的样本数减少.所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量.这就是蒙特卡罗方法中效率的概念.它定义为，其中c 是观察一个子样的平均费用.显然越小，方法越有效.总的来说，增大样本的值对计算机要求较高；减小方差的技巧都只具有指导思想上的意义.对于实际的计算问题，往往要求对涉及的随机变量有先验的了解，或者对发生的物理过程的性态有一定的认识.通过利用这些预知的信息采取相应的手段减小误差，提高精度.附录1.（1）n=1000;p={}Do[m=0;Do[x=Random[];y=Random[];If[x+y<=1,m++],{k,1,n}];AppendTo[p,N[4m/n]],{t,1,10}];Print[p];Sum[p[[t]],{t,1,10}]/10（2）n=10000;p={}Do[m=0;Do[x=Random[];y=Random[];If[x+y<=1,m++],{k,1,n}];AppendTo[p,N[4m/n]],{t,1,10}];Print[p];Sum[p[[t]],{t,1,10}]/10（3）n=100000;p={}Do[m=0;Do[x=Random[];y=Random[];If[x+y<=1,m++],{k,1,n}];AppendTo[p,N[4m/n]],{t,1,10}];Print[p];Sum[p[[t]],{t,1,10}]/102. n=1000;p={}Do[m=0;Do[x=Random[];y=Random[];z=Random[];If[x+y<=1&&x+z<=1&&y+z<=1,m++],{k,1,n}]; AppendTo[p,N[8m/n]],{t,1,10}];Print[p];Sum[p[[t]],{t,1,10}]/103. n=1000;p={}Do[m=0;Do[x=Random[Integer]+2;y=Random[Integer]+1;If[x>y,m++],{k,1,n}];AppendTo[p,N[m]],{t,1,20}]Print[m];{Sum[p[[t]],{t,1,20}]/20,1000-Sum[p[[t]],{t1,20}]/20}参考文献[1] 徐钟济．蒙特卡罗方法[M]．上海：上海科学技术出版社，1985：171-188.[2] 茆诗松，王静龙，濮晓龙．高等数理统计[M].北京：高等教育出版社，2006：415－454．[3] 周铁，徐树方，张平文等．计算方法[M]．吉林：清华大学出版社，2006：299－353．[4] 李尚志，陈发来，张韵华等．数学实验[M]．北京：高等教育出版社，2004：23－30．[5] 王岩．Monte Carlo方法应用研究[J].云南大学学报(自然科学版)，2006，28(S1): 23－26.[6] 薛毅，陈立萍.统计建模与R软件[M].北京：清华大学出版社，2008：476-485.[7] 杨自强.你也需要蒙特卡罗方法———提高应用水平的若干技巧[J]. 数理统计与管理, 2007,27(2):355-376.。
计算多重积分的均匀随机数蒙特卡罗法的实现
时，当有 N 充分大时，有
(7)
将 (7)式 推 广 到 多 重 积 分 ，有
其中 s 为积分变量的个数，N 的取值要充分大。 采用(8)式 进 行 积 分 计 算 时 ，首 先 需 要 获 取 随 机 数 列 ，（t1,t2, …,tk）代 入 被 积 函 数，然后乘上对应积分变量上下限之差的累积，重复 N 次，并累加求和，当 N 充分 大时，所得结果即为积分的近似结果。 使用(8)式编写程序计算积分时，其算法流 程如图 1 所示。 图 1 所示流程可分解为以下 8 个步骤： ① 设 置 变 量 Sum，并 令 其 初 值 为 0，设 置 变 量 N 并 对 其 赋 初 值 ，N 的 取 值 要 求充分大； ② 取 0 到 1 之间均匀分布的随机数列，（t1,t2, …,ts），其中 s 为积分变量的个数； ③ 计算积分变量 xj 的取值 uj，uj 取 值 介 于 aj 和 bj 之 间 ， 由 uj＝aj+(bj-aj)*tj 计 算可实现，j=0,1,2, …,k，aj,bj 为积分变量 xj 的 积 分 上 下 限 ，重 复 s 次 ，求 得 随 机 数 列 （u1,u2,… ,us）； ④ 将 第③步 求 得 随 机 数 列 依 次 代 入 计 算 被 积 函 数 f(x0,x1,…,xs)，求 得 被 积 函 数的值 V；
计算多重积分的均匀随机数蒙特卡罗法的实现
刘辉玲 1，叶锋 2
（1.武汉船舶职业技术学院 电子系，湖北 武汉 430051；2.江汉大学 数学与计算机学院，湖北 武汉 430056）
摘要：蒙特卡罗法是用一系列随机数来近似解决问题的一种方法。 采用均匀随机数蒙特卡罗法计算多重积分是一种简单而有效的 方法，其程序结构简单，易于编制和调试。 本文给出了采用均匀随机数蒙特卡罗法计算多重积分的步骤、算法流程，并给出了实例的 具体实现过程。
计算多重积分的均匀随机数蒙特卡罗法的实现
计算多重积分的均匀随机数蒙特卡罗法的实现作者：刘辉玲叶锋来源：《电脑知识与技术》2008年第35期摘要：蒙特卡罗法是用一系列随机数来近似解决问题的一种方法。
采用均匀随机数蒙特卡罗法计算多重积分是一种简单而有效的方法，其程序结构简单，易于编制和调试。
本文给出了采用均匀随机数蒙特卡罗法计算多重积分的步骤、算法流程，并给出了实例的具体实现过程。
关键词：蒙特卡罗法；多重积分；均匀随机数中图分类号：TP393文献标识码：A文章编号：1009-3044(2008)35-2289-03Realization of Monte Carlo Method by Using Uniform Random Sequence for Calculating Multiple IntegralsLIU Hui-ling1,YE Feng2(1.Department of Electronic and Electrical Engineering, Wuhan Institute of Shipbuilding Technology, Wuhan 430050,China; 2. School of Mathematics and Computer Science, Jianghan University, Wuhan 430056, China)Abstract: Monte Carlo method is a approximate way to solve problems by using a series of random sequence. Monte Carlo method by adopting uniform random number is a simple and effective way to calculate multiple integrals, its structure is simple and easy to program and debug. This paper describes the realization steps and algorithm using Monte Carlo method of by using uniform random sequence to calculate multiple integrals, and gives the realization of concrete examples.Key words: monte carlo method; multiple integrals; uniform random sequence蒙特卡罗法[1-2]，也称为统计试验方法，是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。
7.蒙特卡罗方法在积分计算中的应用——【数学建模 蒙特卡罗算法】
3. 俄国轮盘赌和分裂
1) 分裂
•
设整数 n≥1，令
gi (P) g(P) n
•
则
i Vs gi (P) f (P)dP
n
•
gi(P)于为是原计来算θ的θ的估问计题Vgs，(gP)(可P的化) 1f为/(Pn计)，d算P这n就个i是1θ分i 的i 裂和技来巧得。到，而每个
2) 俄国轮盘赌
N
g(xi , yi )
i 1
•
其方差为
2 gˆ N
1 N
2 x
f1
(
x)dx
•
与通常蒙特卡罗方法相比，方差减少了约
1
N
( x )2 f1(x)dx
6. 分层抽样
•
考虑积分
1
•
g(x) f (x)dx
•
特别地，当 g(P)≥0 时，有
g(P) f (P) g(P) f (P)
f1(P)
g(P) f (P)dP
Vs
•
这时
2 g1
0
•
即 g1的方差为零。实际上，这时有
g1(P) Vs g1(P) f1(P)dP
•
不管那种情况，我们称从最优分布 称函数 | g(P) | 为重要函数。
fl(P)的抽样为重要抽样，
•
我们知道，由f (x,y)抽样 (x,y)的步骤是：
•
从 fl(x) 中抽取 xi，
xi
f1 ( x)dx
1i
•
再由 f2(y|xi) 中抽样确定 yi，
yi
f2(y
xi )dy 2i
•
现在改变 xi 的抽样方法如下：
xi
f1(x)dx i
蒙特卡罗方法与拟蒙特卡罗方法的积分估计
通信系统仿真原理与无线应用课题设计题目蒙特卡洛方法与拟蒙特卡洛方法的几分估计学院电子信息学院专业电路与系统学生姓名贾伟学号 S101832教师徐家品成绩二Ο年月日蒙特卡罗方法与拟蒙特卡罗方法的积分估计贾伟（电路与系统，学号s101832）摘要：在本文中，利用计算机分别产生了伪随机数序列和低差异数序列在此基础上，研究了蒙特卡罗积分与拟蒙特卡罗积分区别关键词：蒙特卡罗积分；拟蒙特卡罗积分；积分估计1 蒙特卡罗方法与拟蒙特卡罗方法简介蒙特卡罗（Monte Carlo）法亦称为随机仿真（random simulation）方法，有时也称作随机抽样（random sampling）技术或统计试验（statistical testing）方法。
蒙特卡罗方法是一种与一般数值计算方法有本质区别的计算方法，属于试验数学的一个分支，起源于早期的用几率近似概率的数学思想，它利用随机数学进行统计试验，以求得的统计特征值（如均值、概率等）作为待解问题的数值解。
这一方法源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯•诺依曼把他和乌拉姆所从事的与研制原子弹有关的秘密工作——对裂变物质的中子随机扩散进行直接模拟，并以摩纳哥国的世界闻名赌城蒙特卡罗作为秘密代号来称呼。
用赌城名比喻随机仿真，风趣又贴切，很快得到广泛接受，此后，人们便把这种计算机随机仿真方法称为蒙特卡罗方法，该方法的基本思想很早以前就被人们所发现和利用。
早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”，而在19世纪人们用投针试验的方法来确定圆周率π。
随着现代计算机技术的飞速发展，用计算机仿真随机过程，实现多次仿真试验并统计计算结果，进而可获得所求问题的近似结果。
蒙特卡罗方法已经在原子弹工程的科学研究中发挥了极其重要的作用，并正在日益广泛地应用于物理、工程、经济、金融的各个方面。
它的基本思想是：为了求解数学、物理、工程技术以及生产管理等方面的问题，首先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求随机参数的统计特征，最后给出所求解的近似值，解的精度可用估计值的标准误差来表示。
蒙特卡罗计算三重积分
蒙特卡罗方法计算三重积分平均值摘要：本文是我读过文献[1]后学习到的一些知识和技巧。
通过阅读文献我对蒙特卡罗方法有了一些新的认识，并对利用蒙特卡罗方法求三重积分有了较好的掌握。
1、蒙特卡罗方法人类通过算法使计算机完成复杂的工作。
人们通过算法设计和程序编制是计算机具有了“智能”。
一般的算法都是确定性算法，即每一工作步骤都是确定的。
还有一些算法叫概率算法，比如说舍伍德算法、拉斯维加斯算法，以及本文将要介绍的蒙特卡洛算法。
这种算法允许在执行过程中随机选择下一个计算步骤。
在某些问题中，随机性算法比确定性算法更有效率。
这种算法的一个特征时，对于一个例子，用该算法计算几次得出的结果可能是不想同的，但是如果我们多计算几次，反复求解，就回降低误差，得到满意的结果。
本文介绍蒙特卡罗方法（随机抽样法，统计实验法）。
它使用概率模型进行近似计算。
蒙特卡罗方法最早起源于18世纪法国学者蒲丰提出的头针问题：…….。
后来，从事原子弹研制工作的冯·诺依曼和乌拉姆将其以摩纳哥的城市名称“蒙特卡罗”作为秘密代号称呼。
2、蒙特卡罗方法的应用蒙特卡罗方法在数学、机械、金融、医学等领域的有着广泛的应用，成功解决了许多经典的数学和物理问题（前文提到的诺依曼从事的就是对裂变物质的中子随机扩散进行直接模拟的工作）。
蒙特卡罗方法还可以用来求解多重积分，解线性方程组和非线性方程组[2]。
下面将主要阐述利用蒙特卡罗方法解三重积分，并给出几个例子。
利用蒙特卡罗方法计算实际问题的基本过程[2]:1）构造并描述随机概率过程。
定积分计算问题本身并不是随机性质的问题，在用蒙特卡罗方法时要首先构造认为的概率过程，将非随机问题转化为随机问题。
2）建立估量，根据随机过程建立某些估计值，即为所求问题的目标解。
3）获取结果，对构造的概率模型进行程序设计，获得预期结果。
第一个过程是很关键的，这个转化过程3、用蒙特卡罗法解三重积分[1]均匀随机数计算三重积分可分为平均值法和掷点法两种，在这里我们只考虑三重积分的平均值法。
多重积分的蒙特卡罗算法编程
s2 =
0.411 839 625 630 18
Elapsed time is 0.047 000 seconds. 在运用了向量化编程后运行时间从 40.094 s 缩 短到 0.047 s, 是原来的 1/853, 效果好; 误差为 0.000 089 920 546 120 003 66, 也比前者好. 若将上面程序稍作如下修改, 即输入参数减少 2 个, 那么当 n 为 1 百万时, 运行时间大约增加 0.45 s, 精度不变.
function s=mtcxg(f, fai1, fai2, a, b, n)
if nargin<6 n=10 000; end
(下转第 6 页)
6
湖 南 文 理 学 院 学 报（自 然 科 学 版）
2011 年
b . G7 × Sn 包 含 一 个 子 图 同 胚 于 Hn , 所 以
cr(G7
quad(fun, a, b) , 这里 x, y 分别是自变量的列向量
和相应被积函数列向量, fun 是被积函数名, 或用
inline()函数直接定义. 虽然上述近似计算能得到比
较满意的结果, 但随着积分重数的增加, 该算法的
计算量显著增加, 以至于计算机也很难完成, 另外
算法的复杂程度要比蒙特卡罗更困难.
x=unifrnd(a, b, 1, n); y=unifrnd(c, d, 1, n);
s=0;
for k=1: n
if and(y(k)>=feval(fai1, x(k)), y(k)<=feval
(fai2, x(k)))
s=s+feval(f, x(k), y(k));
end
end
蒙特卡罗方法计算三重积分
蒙特卡罗方法计算三重积分作者：王洪涛李满枝沈有建来源：《科技视界》 2014年第15期王洪涛李满枝沈有建（海南师范大学数学与统计学院，海南海口 571158）【摘要】应用蒙特卡罗方法的随机投点法原理求解三重积分时，传统的方法是构造四维长方体，然后在该长方体内投点求解，本文创新提出根据被积函数的要求构造出包含积分区域的四维圆柱体并计算。
数值结果表明：四维圆柱区域更贴近于原积分区域，此方法计算效果更好，在三重积分的计算上是个行之有效的方法。
【关键词】蒙特卡罗法；三重积分；随机投点法；四维圆柱体方法Calculation of Triple Integral with Monte-Carlo MethodWANG Hong-tao LI Man-zhi SHEN You-jian(School of Mathematics and Statistics,Hainan Normal University, Haikou Hainan 571158, China)【Abstract】This paper briefly introduces the application of Monte Carlo(MC) method on computing triple integral.it gives one new method by structuring 4-dimension cylinder containing the integral domain. Results of the numerical calculations show that the method is faster and effective than the classic one.【Key words】Monte Carlo method;Triple integral;Randomly cast-point method;4-dimention cylinder me1 蒙特卡罗法的基本原理蒙特卡罗方法以随机模拟和统计试验为手段，从随机变量的概率分布中，通过选择随机数的方法产生一种符合该随机变量概率分布特性的随机数值序列，作为输入变量序列进行特定的模拟试验、求解的方法[1]。
蒙特卡罗方法计算三重积分
Science &Technology Vision 科技视界1蒙特卡罗法的基本原理蒙特卡罗方法以随机模拟和统计试验为手段,从随机变量的概率分布中,通过选择随机数的方法产生一种符合该随机变量概率分布特性的随机数值序列,作为输入变量序列进行特定的模拟试验、求解的方法[1]。
在应用蒙特卡罗方法时,必须要产生非均匀的随机数,而产生特定的、非均匀的随机数序列,可行的方法是先产生一种均匀分布的随机数序列,然后设法转换成我们需要的随机数序列并以此作为数字模拟试验的输入变量序列进行模拟求解[2-6]。
根据三重积分的数学定义,三重积分相当于密度函数为f (x ,y ,z )的空间物体Ω的质量,可将质量假设为四维空间区域的体积。
因为随机投点法需要构造一个能包含积分区域的四维空间区域,我们需要根据该空间区域的体积来计算,因此该方法又叫取体积法。
2蒙特卡罗法计算三重积分的原理及数值算法2.1原理定理[2,6-13]设f (x ,y ,z )为区域Ω上的有界函数,对于三重积分I=f (x ,y ,z )dv ,其中Ω∶a 1≤x ≤b 1,a 2(x )≤y ≤b 2(x ),a 3(x ,y )≤z ≤b 3(x ,y ),满足1)a 2(x ),b 2(x )在区间[a 1,b 1]上连续,a 3(x ,y ),b 3(x ,y )在区域D ∶a 1≤x ≤b 1,a 2(x )≤y ≤b 2(x )上连续,且满足不等式a 2(x )≥c ,b 2(x )≤d,a 3(x ,y )≥e,b 3(x ,y )≤h 2)设M ≥max f (x ,y ,z );3)设四维超长方体Ψ由a 1≤x ≤b 1,c ≤y ≤d,p ≤z ≤q,0≤g ≤M 所围成,体积为V Ψ;4)区域Γ由a 1≤x ≤b 1,a 2(x )≤y ≤b 2(x ),a 3(x ,y )≤z ≤b 3(x ,y ),0≤g≤f (x,y,z )所围;5)在四维超长方体Ψ内产生N 个均匀随机点(x i ,y i ,z i ,g i ),i =1,2,…,N ,设(x i ,y i ,z i ,g i ),i =1,2,…,n 为落在Γ中的n 个随机数,则当N 充分大时,有2.2数值算法1)赋初值:令落在区域Γ中的点的个数n =0;规定投点试验的总次数N ;2)产生四组相互独立的均匀随机数列ξ,η,ω,ψ~U (0,1)由x =a 1+(b 1-a 1)ξ,y =c +(d-c )η,z =p +(q-p )ω,g =Mψ求x i ,y i ,z i ,g i (i =1,2,…,N );3)选出满足a 2(x i )≤y i ≤b 2(x i )的y i (i =1,2,…,N ),设它们是y ki (i =1,2,…,m );4)选出满足a 3(x i ,y i )≤z ki ≤b 3(x i ,y i )且与y ki (i =1,2,…,m )相对应的z ki (i =1,2,…,m ),设它们是z ki (i =1,2,…,l );5)选出满足g ki ≤f (x i ,y i ,z i )且与z ki (i =1,2,…,l )相对应的g ki (i =1,2,…,l ),设它们是g ki (i =1,2,…,n ),则选出的值的个数n 即为所求的落在区域Γ内点的个数。
蒙特卡罗方法
蒙特卡罗方法蒙特卡罗（Monte-Carlo ,简写为M-C ）方法属于计算数学的一个分支， 它是在二十世纪四十年代中期 为了适应当时原子能事业的发展而发展起来的， 但它与一般计算方法有很大区别， 一般计算方法对于解决多维或因素复杂的问题非常困难， 而蒙特卡罗方法对于解决这方面的问题却比较简单。
因而蒙特卡罗方法在近十年来发展很快，特别是随着快速电子计算机的发展，蒙特卡罗方法得到了迅速发展与广泛应用。
蒙特卡罗方法也称随机抽样技术（Random Sampling Technique ）或统计试验方法（Method ofStatistical Test ）。
蒙特卡罗是欧洲摩纳哥国的一个重要城市， 以赌博著称。
蒙特卡罗方法是以概率论与数理统计学为基础的，是通过统计试验达到计算某个量的目的。
而赌博时，概率论是一种有力的手段。
所以，以蒙特 卡罗作为方法的名字，原因大概于此。
由于蒙特卡罗方法是利用一连串的随机数来求解问题的，因此求解随机过程，放射性衰变和布朗运动等问题，它是很有效的。
它除了在原子能工业广泛应用外，在物理、化学、地质、石油、线性规划、 计算机研制、计算机模拟试验、解决多体问题等领域中都有不同程度上的应用。
第一节. 蒙持卡罗方法的基本思想、特点及其局限性一、 蒙特卡罗方法的基本思想用下述三个例子，说明蒙特卡罗方法的基本思想。
例1产品合格率的计算 某工厂生产一批产品，其合格率表示是：为了确定合格率，应该检查这批产品的全部，确定其中合格的数目。
但是，由于产品数量多，检查全部 产品花费的代价大。
因此，通常采取抽取部分产品，在这部分产品中确定其合格的数目。
然后用这部分 产品的合格率F （部分产品合格率） 1 - ■ ™N （部分产品的总数）来代替所要计算的合格率 P 。
例如，检查某批产品，当被检查的产品长度介于 13. 60cm —13. 90cm 内时，则认为是合格的，否则是次品。
分别抽取5件，10件，60件，150件，600件，900件，1200件,1800件来检查，其情况如下表和图 20所示。
第八章 Monte Carlo积分
1. Hit-or-Miss Method
推广到多重积分：
1 2 d I g ( x )dx g ( x1 , x2 ,, xd )dx1dx2 dxd a a1 , a2 ,, ad b b1 , b2 ,, bd g ( x ) 0, h g Method
I Vd E[ g ( X )]
X 在积分域Vd上均匀分布
产生容量为n的X的随机样本Xi,并计算g(Xi),则根据大数定 理，当n足够大时
1 E[ g ( X )] n Vd I In n
g( X )
i 1 i
n
g( X )
2
2. Sample Mean Method 一维积分的情况：
1 I g ( x)dx (b a) g ( x) dx a a ba (b a) E[ g ( X )]
b b
(b a) ˆ I I n (b a) E[ g ( X )] n
2 In
n
where X1, X2, …, Xn are n independent unif(0,3)’s.
2. Sample Mean Method
Simulation Results:
Simulation
1
true = 19.08554, n=100,000
ˆ I
19.10724
2
3 4 5
ˆp ˆ h(b - a) I
ˆ is our estimate of p where p
1. Hit-or-Miss Method
We can easily estimate p:
throw N “uniform darts” at the rectangle
关于计算多重积分的拟蒙特卡罗方法
关于计算多重积分的拟蒙特卡罗方法
张涛;周仲礼;安素珍;张军
【期刊名称】《温州大学学报（自然科学版）》
【年(卷),期】2012(033)005
【摘要】介绍了拟蒙特卡罗方法计算多重积分的基本原理,对Halton序列和Rand 函数产生的序列的均匀性进行了比较,并给出了拟蒙特卡罗方法计算多重积分的步骤、实例、和Matlab语言编写的计算程序.实例充分体现了采用拟蒙特卡罗方法计算多重积分的有效性、精确性和优越性.
【总页数】6页(P33-38)
【作者】张涛;周仲礼;安素珍;张军
【作者单位】成都理工大学数学地质四川省重点实验室,四川成都 610059;成都理工大学数学地质四川省重点实验室,四川成都 610059;成都理工大学数学地质四川省重点实验室,四川成都 610059;成都理工大学数学地质四川省重点实验室,四川成都 610059
【正文语种】中文
【中图分类】O242
【相关文献】
1.使用拟蒙特卡罗方法计算点模型的体积 [J], 刘玉身;雍俊海;张慧;杜明翠;孙家广
2.计算多重积分的蒙特卡罗方法与数论网格法 [J], 宫野
3.利用蒙特卡罗方法与拟蒙特卡罗方法计算定积分 [J], 韩俊林;任薇
4.基于拟蒙特卡罗方法的Copula-GARCH类模型在外汇风险计算中的应用 [J], 王梦媛;邹玉梅
5.蒙特卡罗方法的原理及在数值计算方面的应用\r——以定积分为例 [J], 李姣娜因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
【８】宫野．计算多重积分的蒙特卡罗方法与数论网格法【Ｊ】．大连理工大学学报，２００１，４ｌ（１）：２１—２２， ［９】Ｗａｎｇ
Ｘ
ｅ，Ｆａｎｇ
Ｋ Ｔ Ｔｈｅ Ｅｆｆｅｃｔｉｖｅ Ｄｉｍｅｎｓｉｏｎ ａｎｄ Ｑｕａｓｉ—Ｍｏｎｔｅ Ｃａｒｌｏ
［２０１２－０１—０８１．ｈｔｔｐ：／／ｃｉｔｅｓｅｅｒｘ．ｉｓｔ．ｐｓｕ．ｅｄｕ／ｖｉｅｗｄｏｃ／ｄｏｗｎｌｏａｄ？ｄｏｉ＝１０．１．１．３７．９３９６＆ｒｅｐ＝ｒｅｐｌ＆ｔｙｐｅ２ ［６］Ｋｒｙｋｏｖａ
Ｉ．Ｅｖａｌｕｔｉｎｇ ｏｆ Ｐａｔｈ・ｄｅｐｅｎｄｅｎｔ Ｓｅｃｕｒｉｔｉｅｓ ｗｉｔｈ Ｌｏｗ Ｄｉｓｃｒｅｐａｎｃｙ
万方数据
张涛等：关于计算多重积分的拟蒙特卡罗方法
ｔ’ｏｒｉ＿１：ｎ
３５
ｊ＝ｌ；
ｏｋ＝０： ｗｈｉｌｅ ｏｋ＝＝０
ＷｏｒｋＶｅｔ（ｊ）＝ＷｏｒｋＶｅｔ（ｊ）＋ｌ； ｉｆＷｏｒｋＶｅｔ（ｊ）＜ｂａｓｅ
ｏｋ＝ｌ： ｅｌｓｅ
ＷｏｒｋＶｅｔ（ｊ）＝０； Ｊ＝ｊ＋ｌ；
ｅｎｄ ｅｎｄ
Ｓｅｑ（ｉ）＝ｄｏｔ（ＷｏｒｋＶｅｔ，ＶｅｔＢａｓｅ）；
如Ｈａｌｔｏｎ序列、Ｓｏｂｏｌ’序列、Ｎｉｄｅｒｒｅｉｔｅｒ的（ｆ，ｍ，ｓ）网格和（ｆ，ｓ）序列等［６１．本文主要研究Ｈａｌｔｏｎ 序列在计算多重积分中的应用． ２
Ｈａｌｔｏｎ序列
Ｈａｌｔｏｎ序列是Ｖａｎ
ｄｅｒ
Ｃｏｒｐｕｔ序列的推广．Ｈａｌｔｏｎ序列是通过将一系列整数表示成某个基的
数位的形式，然后将这些数位按反序排列再在前面加小数点而得到的值．本文将Ｓ维Ｈａｌｔｏｎ序列 表示成ｘ１，ｘ２，…，其中每一个随机数是一个Ｓ维向量，即ｔ＝（Ｘｉｌ，Ｘｉ２，…，Ｘｉｓ）．生成Ｈａｌｔｏｎ序列
对于单位超立方体１５＝［０，１］５上的积分，拟蒙特卡罗估计的形式和蒙特卡罗类似，即：
弘胁专善他，，
只是此处的点列葺，…，ｘＮ∈１５为确定性的点． 因此，拟蒙特卡罗方法的基本思想是用精选的确定 性点来代替蒙特卡罗方法中的随机点． 拟蒙特卡罗方法积分的误差估计是根据ＫｏＫｓｍａ．Ｈｌａｗｋａ不等式来给出，即令厂为区域
Ｍｅｔｈｏｄｓ［ＤＩ．Ｍａｓｓａｃｈｕｓｅｔｔｓ：Ｗｏｒｃｅｓｔｅｒ
Ｐｏｌｙｔｅｃｈｕｉｃ Ｉｎｓｔｉｔｕｔｅ，２００３：７－３０．
［７］Ｂｒａｎｄｉｍａｒｔｅ Ｐａｏｌｏ．Ｎｕｍｅｒｉｃａｌ Ｍｅｔｈｏｄ ｉｎ Ｆｉｎａｎｃｅ：Ａ ＭＡＴＬＡＢ—ｂａｓｅｄ Ｉｎｔｒｏｄｕｃｔｉｏｎ【Ｍ］．Ｎｅｗ Ｙｏｒｋ：Ｊｏｈｎ Ｗｉｌｅｙ ａｎｄ
的估计值．
表１
２１５ ６６４ ９
取Ｎ＝２ ０００，在计算机上分别用Ｈａｌｔｏｎ序列和Ｒａｎｄ函数计算的结果见．表１．
分别用Ｈａｌｔｏｎ序列和ｒａｎｄ函数计算的积分结果
本例中仅给出Ｓ＝３时采用Ｈａｌｔｏｎ序列的Ｍａｔｌａｂ程序：
Ｓ＝３；Ｎ＝２０００；
ａ（：，１）２ＧｅｔＨａｌｔｏｎ（Ｎ，２）；ａ（：，２）＝ＧｅｔＨａｌｔｏｎ（Ｎ，３）；ａ（：，３）＝ＧｅｔＨａｌｔｏｎ（Ｎ，５）；
ｒ
１
ｇ（＿，屯，…，ｔ）：Ｊ寺，ｇ（五，Ｘ２，…，Ｘｓ）∈Ｄ，
１ ４算例及其实现
下面给出两个特殊数值算例验证拟蒙特卡罗方法在计算多重积分中的有效性和优越性． 例１用Ｈａｌｔｏｎ序列计算Ｓ重积分
０，ｇ（ｘｌ，Ｘ２，…，ｔ）萑Ｄ
这里Ｄ也表示积分区域的体积．
Ｇ＝ｆｆ．一ｆｉｎ Ｓｌｎｍａｘ（ｘ。州Ｘ２一，ｔ）ｂ…咄＝一０．５７７
Ｉ。＝［Ｏ，Ｉ】。上Ｈａｒｄｙ Ｋｒａｕｓｅ意义变分Ｖ（ｆ）有界的实函数，对任意Ⅳ个点Ｘ张涛（１９８６），男，河南信阳人，硕士研究生，研究方向：数学计算方法
万方数据
３４
温州大学学报・自然科学版（２０１２）第３３卷第５期
卡罗方法模拟解积分问题的误差界可以表示为‘３４】：
ｅｎｄ
瞳机函数ｒａｎｄ产生的＋雉随机数点 Ｈａｌｔｏｒｌ序列产生的＋雉睫机披点
图ｌ
Ｈａｌｔｏｎ序歹Ｕ与Ｒａｎｄ函数ｊ“生的十维随机数点的均匀性比较
由图ｌ可明显看出，Ｒａｎｄ函数产牛的伪随机数序列有抱团现象，分布不均匀，而Ｈａｌｔｏｎ序 列则要均匀得多．有了这样优质的均匀随机序列，用拟蒙特卡罗方法求解得到的将是确定性的误
ＤＯＩ：１０．３８７５牙．ｉｓｓｎ．１６７４—３５６３．２０１２．０５．００６
在数学和工程科学计算中，求解多重积分的近似值是一个至关重要的环节．通常人们所采用 的方法有３种，即蒙特卡罗方法、拟蒙特卡罗方法和数论网格方法．然而，在实际应用中，数论 网格法很难解决高维问题，如对２０维这样的问题，就要有一百万多个点［１］，计算量大体上随维数 的幂次增加，几乎是不可能计算的．蒙特卡罗方法是采用单和来得到多重积分的近似值的［２１，计 算跟问题的维数无关，克服了“维数灾难”，但其误差却是概率性的，得到解的精度不高．本文 研究了拟蒙特卡罗方法在计算多重积分中的应用，并给出Ｍａｔｌａｂ程序的实现．实例表明，采用拟 蒙特卡罗方法计算多重积分，既能克服维数的制约，又能得到确定性的误差估计和更高精度的解．
ｆｕｎｃｔｉｏｎ ５００，０．２５９ ２５９，Ｏ．１２０
ｏｏｏ）．将ｍ加
Ｓｅｑ＝ＧｅｔＨａｌｔｏｎ（ｎ，ｂａｓｅ）
Ｓｅｑ＝ｚｅｒｏｓ（ｎ，１）；
ＮｕｍＢｉｔｓ
２
１＋ｃｅｉｌ（１０９（ｎ）／１０９（ｂａｓｅ））；
ＶｅｔＢａｓｅ＝ｂａｓｅ．＾（一（１：ＮｕｍＢｉｔｓ））； ＷｏｒｋＶｅｔ＝ｚｅｒｏｓ（１．ＮｕｍＢｉｔｓ）；
Ｊ，Ｃａｆｌｉｓｃｈ Ｒ
Ｅ．Ｑｕａｓｉ—ｒａｎｄｏｍ
Ｓｅｑｕｅｎｃｅｓ ａｎｄ Ｔｈｅｉｒ
Ｄｉｓｃｒｅｐａｎｃｉｅｓ［Ｊ】ＳＩＡＭ Ｊ
ＳＣＩ Ｃｏｍｐｕｔ，１ ９９４，１
５（６）
万方数据
３８
１２５ｌ一１２７９．
温州大学学报－自然科学版（２０１２）第３３卷第５期
［４］Ｎｉｅｄｅｒｒｅｉｔｅｒ ｆ５］Ｍｏｒｏｂｏｓｉ
ｆｏｒ ｉ＝ｌ：２０００
ａ＿ｍａｘ（ｉ）＝ｍａｘ（ａ（ｉ，：）；
万方数据
张涛等：＆５－计算多重积分的拟蒙特卡罗方法
３７
ｂ＝Ｓ丰ｌｏｇ（ａ＿ｍａｘ）；Ｃ Ｎ＝ａｂｓ（ｂ）；Ｄ Ｎ＝ｌｏｇ（Ｃ－Ｎ）； Ｇ＿Ｎ＝ｍｅａｎ（Ｄ＿Ｎ）．
靴月Ｈａｌｔｏｎ蒯计夥重积分Ｇ牛聪等警”咄＝１的估襻． …’ｏ
０ ０
取Ｎ＝１ ０００，ａ１＝ａ２＝…＝臼，＝０，在计算机上分别用Ｈａｌｔｏｎ序列和Ｒａｎｄ函数计算的结果 她表２．
ｌ！ｆ（ｘ）ｄｘ－专善讹）ｈ门成（矿＾），
阶‘５１是Ｏ（Ｎ。１），而蒙特卡罗积分的误差的阶是Ｏ（Ｎ。１佗）． 特卡罗积分要优于蒙特卡罗积分．对于更高维的积分问题，拟蒙特卡罗法优于数论网格法．
㈣
含Ⅳ个点的拟随机序列的偏差球（Ｘｌ，…，ＸＮ）为Ｏ（Ｎ一（１０９Ｎ）”１），所以拟蒙特卡罗积分的误差
对于拟蒙特卡罗方法积分，我们有确定性的误差估计．在计算积分中还可以明智地选择这些 确定的点来减小式（１）中的误差，得到高精度的解．因此，就确定性和高精度性这两点，拟蒙 在过去的几十年里，拟蒙特卡罗法得到了快速发展，人们已经构造出大量优质的低差异序列，
２）将数位按反序排列再在前面加小数点得到新的值
吻＝Ｚａｍ。ｂ；叫＂一，／＝１，…，Ｓ；
３）置ｍ＝ｍ＋１，然后重复１）和２）两步．
举例，如Ｓ＝３，即维数是３，设ｍ＝１５，选择２，３和５作为基，有１５＝１１１１２，１５＝１２０３，１５＝３０５，
于是便有第一个随机变量ｘ＝（０．１１ １１２，０．０２１ ３，Ｏ．０３５），即（Ｏ．９３７ １，如此循环，便得到整个序列． 下面将给出产生Ｈａｌｔｏｎ序列的Ｍａｔｌａｂ程序［７］，并对Ｈａｌｔｏｎ序列与Ｒａｎｄ函数产生的随机数序 列的均匀性进行比较，见图１．
第３３卷第５期
Ｖｂｌ ３３．ＮＯ ５
温州大学学报・自然科学版
Ｊｏｕｒｎａｌ ｏｆＷｅｎｚｈｏｕ Ｕｎｉｖｅｒｓｉｔｙ’Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅｓ
２０１２年１０月
０ｃｔ．２０１２
关于计算多重积分的拟蒙特卡罗方法
张涛，周仲礼，安素珍，张军
（成都理工大学数学地质四川省重点实验室，四川成都
６１００５９）
的简单易行的步骤如下：首先选择Ｓ个基岛，ｂ２，…，瓦，比如选择前ｓ个素数，然后对某个整数ｍ，
将ｍ表示成以ｂｉ为基的数位，将这些数位按反序排列再在前面３ｎ４，数点得到新的值，便是序列
中某个随机数向量的第，个元素．即对某个整数ｍ有以下步骤［６１：
１）选择适当大的ｔ。ｉ，对每个基将ｍ表示成：
撤＝∑ａ。々ｂ，ｋ，ｊ＝ｌ，－－－，ｓ；
ＪＦ ｓｕｍ＝０： ｆｏｒ ｉ＝ｌ：Ｎ ＬＣ＝Ｉ；
ｆｏｒｊ＿１：Ｓ ＬＣ＝ＬＣ水ａｂｓ（４术ａ（ｉ，ｊ）一２）；
ｅｎｄ
ＪＦ＿ｓｕｍ＝ＪＦ＿ｓｕｍ＋ＬＣ；
ｅｎｄ ＪＦ＝ＪＦ—ｓｕｍ／Ｎ．
５结
语
通过计算可以看出，拟蒙特卡罗方法是有效的，积分重数对积分结果的误差无明显影响，精
度较高，计算机程序实现也是很简单的．本文只是用了原始的Ｈａｌｔｏｎ序列，由于Ｈａｌｔｏｎ序列基
１拟蒙特卡罗方法基本原理及其优越性
拟蒙特卡罗方法也称低差异序列法，是在蒙特卡罗方法基础上发展起来的一种模拟方法．拟 蒙特卡罗方法与蒙特卡罗方法相似，但所用的理论基础却不同．拟蒙特卡罗方法是通过构造所谓 的低差异序列，即用的是确定性的点，然后用Ｋｏｋｓｍａ．Ｈｌａｗｋａ不等式来确定误差阶的，而不是根 据大数定律ｐ】．
Ｈ．ＲａｎｄｏｍＮｕｍｂｅｒ Ｇｅｎｅｒａｔｉｏｎ ａｎｄ
Ｑｕａｓｉ Ｍｏｎｔｅ Ｃａｒｌｏ Ｍｅｔｈｏｄｓ［Ｍ］．Ｐｈｉｌａｄｅｌｐｈｉａ：ＳＩＡＭ，１９９２：１０—２１． Ｉｎｔｅｇｒａｔｉｏｎｓ【ＥＢ／ＯＬ］．
ｐｄｆ．
Ｈ，Ｆｕｓｂｉｍｉ Ｍ，Ａ Ｐｒａｃｔｉｃａｌ Ａｐｐｒｏａｃｈ ｔｏ ｔｈｅ Ｅｒｒｏｒ Ｅｓｔｉｍａｔｉｏｎ ｏｆ Ｑｕａｓｉ—Ｍｏｎｔｅ Ｃａｒｌｏ