蒙特卡罗随机模拟投点法

蒙特卡罗随机模拟投点法在数字积分中的

应用

数学与应用数学0901班：张瑞宸

指导老师：任明慧

摘要：本文首先介绍了蒙特卡罗方法的产生和发展，然后分析了蒙特卡罗方法计算数值积分的理论原理，最后给出了蒙特卡罗方法计算数值积分的MATLAB编程实现，全文主要是讨论了蒙特卡罗方法在定积分计算的应用。而蒙特卡罗的优点：可以计算被积函数非常复杂的定积分、重积分，并且维数没有限制，这是别的数值积分方法还未达到的。蒙特卡罗的缺点：收敛速度慢，误差一般较大，且是概率的误差，不是真正的误差。

关键词：蒙特卡罗方法，均值估计法，数值积分，Matlab编程

Abstract：This paper first introduces the emergence and development of the Monte Carlo method, and then analyze the theoretical principles of Monte Carlo numerical integration method, Full-text mainly discussed the application of the Monte Carlo method in the definite integral. The advantages of Monte Carlo: can be calculated the integrable functions very complex definite integral, Multiple integrals, and dimension no limit, other numerical integration methods have not yet reached. Monte Carlo Disadvantages: slow convergence speed, error generally higher, and the probability of error, not a real error.

Keywords: Monte Carlo method，Mean estimation method，numerical integral，Matlab programming

0 引言

历史上有记载的蒙特卡罗试验始于十八世纪末期（约1777年），当时布丰（Buffon）为了计算圆周率，设计了一个“投针试验”，后文会给出。虽然方法已经存在了200多年，此方法命名为蒙特卡罗则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用蒙特卡罗方法模拟中子对某种特殊材料的穿透作用。出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一von Neumann灵犀一点选择摩洛哥著名赌城蒙特卡罗（Monte Carlo）作为该项目名称，自此这种方法也就被命名为Monte Carlo方法广为流传。

蒙特卡罗方法，又名随机模拟法或统计实验法它是以概率统计理论为基础，依据大数定律（样本均值替代总体均值）利用电子计算机数字模拟技术，解

决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的求出数值解的方法。一般蒙特卡罗法在数学中最常见的应用就是蒙特卡罗随机投点和蒙特卡罗数值积分。

1 蒙特卡罗方法的产生与发展

蒙特卡罗方法是在二战期间产生和发展起来的他的奠基者是美籍匈牙利人数学家冯诺伊曼（J.Von Neumann 1903-1957）由于通常计算量相当大而电子计算机在当时还没有出现，所有运算只能用手工进行，故而相当长的时间里蒙特卡罗方法难以推广。

1.1 蒙特卡罗方法的产生

作为蒙特卡罗方法的最初应用，是解决蒲丰氏问题1777 年，法国数学家Buffon 提出利用投针实验求解的问题：

设平面上有无数多条距离为 1 的等距平行线，现向该平面随机的投掷一根长度为()1l l ≤的针。随机投针是指针的中心点于最近的平行线间的距离x 均匀分布在[0,1/2]上，针与平行线的夹角?（不管相交与否）均匀分布在[]0,π上，如图一所示，故而，故其1~(0,)2

x U ，~(0,)U ?π概率密度函数分别为11(),()2p x p ?π

==。 故我们得到针与线相交的充要条件是2

sin l x ≤?

图1.1 针与线相交的几种情况

则针与线相交的概率是(sin )2

l p x ?≤ =sin sin 2200002()()l l p x p dxd π

?π???π=???? =2

22.sin 2l l d π??ππ=? 所以得到圆周率222(sin )l l l p x p

π?==≤。 假如我们能做大量的投针实验并记录下针与线的相交次数，则可以根据大数定律估计出针线相交的概率P 。投针实验N 次可能有n 次使针与任意平行线相交, 那么n p N

≈

，显然，试验次数N 越多，P 的近似程度越好。历史上曾有几位学者做过这样的投针试验，并用手工计算出π值，结果参见表1。 表 1

应该指出，上述试验的精度一般不会很高。譬如，假设1,p l π==

则。则由中心极限定理，如果试验次数为N ，则p 的估计值?p 渐进服从(1)(,)p p N p N

-，近似为0.2313(0.6366,)N N

。因此，若要以95%的概率保证p 的精确到三位有效数字，即?p 与p 的差距小于0.001，则N 必须满足N ≥2251.960.2313/0.0018.8910?≈?。 重复进行上千次的投针实验和手工计算，要消耗大量的人力、财力，蒙特卡罗方法虽然能解决此类问题，但得不到推广应用。

1.2 蒙特卡罗方法的发展

20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算

处理，于是Monte Carlo 方法重新被提起，引起世人重视，应用日渐广泛。实际上，采用Monte Carlo 方法在计算机上建立模型来解决Buffon 问题是非常简单的。

我们在计算机上进行模拟试验，给定l ，我们可以在计算机上随机产生x 和φ，然后判断2

sin l x ≤?是否成立。若成立，则针线相交，否则不交。假如我们在计算机上独立的产生N 对这样的x 和?，并记录下

2sin l x ≤?成立的次数,记为n 。则π估计值可取为2

2l N ??，这就是随机模拟计算的结果。 借助Matlab 在计算机上进行投针实验，（相应程序见附录1）

所取得的计算结果参见表2。

表 2

不难看出，利用计算机运算不仅结果精确，而且迅速。

随着电子计算机的普及，蒙特卡罗方法作为一种独到的方法得到开发，并首先应用到核武器的试验与研制中。尤其是各种可是编程方法的不断涌现，更显示出蒙特卡罗方法的最独到的优点，即形象直观地用数学方法在电子计算机上实现数字模拟实验。

2 蒙特卡罗方法计算数值积分的理论原理

随机模拟是一种随机试验的方法，也称为蒙特卡罗方法（Monte Carlo ）]1[。这种方法利用随机试验，根据频率与概率、平均值与期望值等之间的关系，推断出预期的结果。

2.1 蒙特卡罗算法的理论原理

随机模拟法的基本原理非常简单，下面给以直观的说明。 首先构造一个定积分π=-?10214dx x ，其中积分?-=1

021dx x S （即1/4单位圆的面积）用随机模拟方法求得。

图2.1中粗线是1/4单位圆，如果向图2.1

中边长为1的正方形里随机投n 块小石头，当n

很大时小石头会大致均匀地分布在正方形中，

数一下落在1/4单位圆内的小石头，假定有k

个，那么k /n 就可以看作1/4单位圆面积π/4

的近似值。 小石头的位置坐标可以用产生均匀分布随机数的程序得到，记为),...,2,1(,n i y x i i =是[0,1]区间均匀分布随机数（i i y x ,相互独立），记录满足),...,2,1(122n i y x i i =≤+的数量k ，即得π=4k /n 。

我们用概率论中的大数定律来说明这个直观认识的原理。

大数定律（伯努利（Bernoull ）定理）]2[ 设k 是n 次独立重复试验中事件A 发生的次数，p 是事件A 每次试验中发生的概率，则对任意的正数ε，有

1lim =?

?????<-∞→εp n k P n （1） 若规定“向图中正方形随机投一块小石头落在四分之一单位圆里”为事件A 发生，则A 发生的概率p 应该等于四分之一单位圆面积，随机投n 块石头就是独立重复做n 次试验，事件A 发生k 次，由（1）式，n 无限变大时k /n 与p 之差小于任意一个数ε的概率趋于1。

用计算机算一下就会发现，即使n 很大结果也不好，并且很不稳定，远不如常规数值积分的几种方法。实际上，随机模拟法很少用来做定积分，它的特点是能够方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。

2.2 蒙特卡罗方法计算定积分的理论原理

通常有两类办法计算定积分。

（1）随机投点法

在“投石算面积”的例子中，事件A 在每次试验中发生的概率p 是四分之一单位圆面积，即

dx x dydx p x ???-==-1021

01012 （2）

n 次试验由计算机完成，采用[0,1]区间上的均匀分布产生相互独立的随机数。记这样产生的n 个点的坐标为),...,2,1),(n i y x i i =。事件A 发生等价于),(i i y x 满足21i i x y -≤，A 发生的次数是满足21i i x y -≤),...,2,1(n i =点的个数k 。由伯努利定理，p 可以用k /n 近似替代。

这种方法可以推广如下。对任意区间],[b a 内的连续函数)(x f ，满足d x f c ≤≤)(,(c )0≥d ，为计算定积分?b

a dx x f )(，先由计算机产生n 个点的坐标),(i i y x ),...,2,1(n i =，

其中i i y x ,分别为],[b a 和],[d c 区间上的均匀分布随机数，然后记录n 个点中满足)(i i x f y ≤的数目k ，则

c a b n

k c d a b dx x f b a )())(()(-+--≈?

（3）

（2）均值估计法 这种方法依据概率论的一下两个定理：

大数定理（辛钦定理）]2[ 若随机变量n Y Y Y ,...,,21相互独立，服从同一个分布，且具有数学期望μ=i EY ),...,2,1(n i =,则对任意的正数ε，有

11lim 1=?

?????<-∑=∞→εμn i i n Y n P （4） 随机变量函数的期望]2[ 若随机变量X 的概率分布密度是))((b x a x p ≤≤，则随机变量的函数)(X f Y =的期望为

dx x p x f X f E b

a ?=)()())(( （5）

于是，当X 为],[b a 区间均匀分布的随机变量时，)/(1)(a b x p -=)(b x a ≤≤，（5）式给出

))(()()(x f E a b dx x f b

a -=? （6）

只要产生],[b a 区间相互独立、均匀分布的随机数i x ),...,2,1(n i =，)(i i x f y =就相互独立。根据辛钦定理，当n 很大时期望))((X f E EY =就可以用)(i x f 的平均值近似，所以由（6）式得到

∑?=-≈n i i b

a x f n a

b dx x f 1

)()()( （7） 与随机投点法相比，均值估计法没有d x f c ≤≤)(,(c )0≥d 的限制，只需计算)(i x f ，不需要产生随机数i y ，也不需要作i y ≤)(i x f 的比较，显然大为方便。

2.3 蒙特卡罗方法计算重积分的理论原理

用随机模拟作数值积分的优点不仅在于计算简单，尤其是它可以方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。而通过上文对蒙特卡罗计算定积分的两种方法的比较，我们可知均值估计法比投点法更好操作，经过计算还可知均值估计法精确度（若选择的密度函数适合）会更高。重积分本文讨论均值估计法。

考虑k 维积分：

k k dx dx dx x x x f I ...),...,,(...2121??Ω

=

其中Ω为k 维积分区域，V ?Ω(i i i b x a V ≤≤:，k i ,...,2,1=)。

选取Ω上的一个概率密度函数),...,,(21k x x x p ，V ?Ω(i i i b x a V ≤≤:，

k i ,...,2,1=)，且),...,,(21k x x x p 满足1...

),...,,(...2121=??Ω

k k dx dx dx x x x p ，Ω∈),...,,(21k x x x ，故：

k k k k dx dx dx x x x p x x x p x x x f I ...),...,,()

,...,,(),...,,(...21212121??Ω= ))

,...,,(),...,,((2121k k x x x p x x x f E = 其中,Ω∈),...,,(21k x x x ∑=≈m i k k x x x p x x x f m 12121)

,...,,(),...,,(1 其中,Ω∈),...,,(21k x x x （8） 这里m 是落在Ω区域中的随机投点数。

对于二重积分：

??Ω

=dxdy y x f I ),(

其中b x a ≤≤Ω:，d x g y x g c ≤≤≤≤)()(21。记],[],[d c b a V ?=，设ΩS ，V S 分

别是Ω，V 的测度，若选取Ω是),(y x g 上的均匀分布，则Ω∈=Ω

y)(x, , 1),(S y x g 。由几何概率得n

m S S V ≈Ω，即n S m S V ≈Ω，于是： ??Ω=dxdy y x g y x g y x f I ),()

,(),( ∑==m i y x g y x f m 1)

,(),(1 ∑=Ω

=m i y x f m S 1

),(∑==m i V y x f n S 1),( 综上可得， ∑==m

i V y x f n

S I 1),( （9）

3蒙特卡罗方法MATLAB编程实例

用随机模拟作数值积分主要利用均匀分布产生的随机数及相应的判断、计算和简单运算。MATLAB提供的命令

unifrnd(a,b,m,n)

产生m行n列]

a区间上均匀分布随机数。当0

[b

,

=

b时，可用rand(m,n)。

=

a，1

例1用蒙特卡罗方法计算π。

解(1) 随机投点法。MATLAB程序如下：

n=10000;

x=rand(2,n);

k=0;

for i=1:n

if x(1,i)^2+x(2,i)^2<=1

k=k+1;

end

end

p=4*k/n

重复计算4次，计算结果分别为：

p=3.1244 P=3.1420 P=3.1780 P= 3.1592

当n提高到50000时，重复计算4次，计算结果为：

p=3.1372 P=3.1456 P=3.1442 P=3.1426

解(2) 均值估计法，此时2

x

f-

=，MA TLAB程序如下：

(x

)

1

n=50000;

x=rand(1,n);

y=0;

for i=1:n

y=y+sqrt(1-x(i)^2);

end

p=4*y/n

重复计算4次，计算结果分别为：

p= 3.1439 P= 3.1510 P= 3.1474 P= 3.1449

当n提高到50000时，重复计算4次，计算结果为：

p= 3.1471 P=3.1419 P=3.1423 P= 3.1430

可以看出它的缺点是计算量大，结果具有随机性，精度较低。一般说来精度为2/1-

n阶，n增加时精度提高较慢。但是用随机模拟方法可以计算被积函数非常复杂的定积分、重积分，并且维数没有限制。如下介绍多重积分的MALTAB实现。

3.1 多重积分的蒙特卡罗方法MATLAB 编程实现

以二重积分为例，运用算法（8）编程，一个是通常的for 循环编程，文件名为mtc.m ，另一个是向量化编程，文件名为mtcx.m ，后者运行速度要快很多。

function s=mtc(f,fai1,fai2,a,b,c,d,n)

if nargin<8 n=10000;end

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

s=s/n*(b-a)*(d-c);

return

function s=mtcx(f,fai1,fai2,a,b,c,d,n)

if nargin<8 n=10000;end

x=unifrnd(a,b,1,n);y=unifrnd(c,d,1,n);

s=sum(feval(f,x,y)).*and(y>=feval(fai1,x),y<=feval(fai2,x)));

s=s/n*(b-a)*(d-c);

return

三重或三重以上积分只须增加若干个输入参数，即新增积分变量的上、下限，方法不变，对程序稍作修改即可。

例 2 计算??-----+1212121222

22)sin(dydx y x e x x x

积分的精确值为0.41192954617630，在命令窗口调用函数M 文件mtc.m 和mtcx.m 情况如下：

>>tic;s1=mtc(inline(‘exp(-x.^2/2).*sin(x.^2+y)’),

inline(‘-sqrt(1-x.^2/2)’),inline(‘sqrt(1-x.^2/2)’),-0.5,1,-1,1,100000),toc

s1=

0.409 010 343 471 83

Elapsed time is 40.094 000 seconds.

>>tic;s2=mtcx(inline(‘exp(-x.^2/2).*sin(x.^2+y)’),

inline(‘-sqrt(1-x.^2/2)’),inline(‘sqrt(1-x.^2/2)’),-0.5,1,-1,1,100000),toc s2=

0.411 839 625 630 18

Elapsed time is 0.047 000 seconds.

在运用了向量化编程后运行时间从40.094s缩短到0.047s，是原来的1/853，效果好；误差为0.00 008 992 054 612 000 366，也比前者好。

若将上面程序稍作如下修改，即输入参数减少2个，那么当n=1 000 000时，运行时间大约增加0.45s，精度不变。

function s=mtcxg(f,fai1,fai2,a,b,c,d,n)

if nargin<6 n=10000;end

x=unifrnd(a,b,1,n);

c=min(feval(fai1,x));d=max(feval(fai2,x));

y=unifrnd(c,d,1,n);

s=sum(feval(f,x,y)).*and(y>=feval(fai1,x),y<=feval(fai2,x)));

s=s/n*(b-a)*(d-c);

return

>>tic;s3=mtcxg(inline(‘exp(-x.^2/2).*sin(x.^2+y)’),

inline(‘-sqrt(1-x.^2/2)’),inline(‘sqrt(1-x.^2/2)’),-0.5,1,-1,1,100000),toc s3=

0.412 005 502 840 08

Elapsed time is 0.922 000 seconds.

对比常规数值积分和蒙特卡罗方法数值积分，同样数量的n值——也就意味这几乎相同的计算量——常规数值积分结果的精确度要高于蒙特卡罗数值积分的结果。那么，我们为何还需要用蒙特卡罗来算数值积分呢？答案的关键在于，常规数值积分的精度直接取决于每个维度上取点数量，维度增加了，但是每个维

度上要取的点却不能减少。在多重积分中，随着被积函数维度增加，需要计算的函数值数量以指数速度递增。

例如在一重积分?a

b dx x x f )()(ψ中，只要沿着x 轴取n 个点；要达到相同大小的精确度，在k 重积分k n k dx dx dx x x x x x x f ...)...()...(...212121???ψ中，仍然需要在每个维度上取n 个点，k 个纬度的坐标相组合，共需要计算n ×k 个坐标对应的)(x f 函数值。取点越多，会占用计算机大量内存，也需要更长运算时间，最终导致这种计算方法不可行！

蒙特卡罗方法却不同，不管是积分有多少重，取n 个点计算的结果精确度都差不多。因此，即使在一重积分的情形下，蒙特卡罗方法的效率比不过常规数值积分，但随着积分维度增加，常规数值积分的速度呈指数下降，蒙特卡罗方法的效率却基本不变。经验表明，当积分重数达到4重积分甚至更高时，蒙特卡罗方法将远远优于常规数值积分方法。

4 总结

归纳起来本论文主要论述了：蒙特卡罗算法在数值积分中的运用。本篇论文给出了蒙特卡罗方法应用在数值积分中的原理，讨论蒙特卡罗算法运用在数值积分中的两种计算方法，随机投点法和均值估计法，具体蒙特卡罗算法的实现给出了MATLAB 编程的实例。最后对蒙特卡罗算法的优缺点进行了讨论。

蒙特卡罗理论的应用非常广泛,两种蒙特卡罗数值计算方法的应用也很全面,由于作者的知识水平,研究能力有限,蒙特卡罗理论在各行各业中都有广泛使用,思想、理论等方面仍存在欠缺之处,但这些都是以后的努力研究的方向,有侍继续发掘。

参考文献：

[1]曲双石，王会娟.MonteCarlo方法及其应用[J].统计教育，2009.1.

[2]茆诗松，程依明，濮晓龙.概率论与数理统计教程[M].高等教育出版社，2004年7月.

[3]高雷阜.Monte-Carlo理论与优化方法的研究[J].辽宁工程技术大学学报(自然科学版).2002年03期.

[4]李世刚,刘辉,陈标华.素数的一个特殊性质及其用于伪随机数生成的方法[J].北京化工大学学报(自然科学版).2003年03期.

[5]张韵华，奚梅成，陈效群.数值计算方法与计算[J]---2版.--北京科学出版社，2006.

[6]孙维君，秦华.Monte Carlo方法在三重积分中的应用[J].山东理工大学学报（自然科学报），2008第1期.

[7]柴中林，银俊成.蒙特卡罗方法计算定积分的进一步讨论[J].应用数学与计算数学学报.2008.第1期.

[8]张韵华，奚梅成，陈效群.数值计算方法与算法[M].科学出版社2006.第2版.

[9] JONATHANJ.FORSTER ,JOHNW.M CDONALD and PETERW.F.SMITH.Monte Carlo methods to calculate[J].2005.42,37-46.spring-verlag.

致谢:首先要感谢以梅院长等领导对我们的精心统筹安排,科学培养,使我们能够完成所学的知识积累,完成论文。我更要感谢在论文完成过程中一直给与悉心指导的任老师，任老师在指导论文的过程中,耐心的讲解和教诲，使我能掌握写作本文的基础知识并理解掌握参考文献的写作的框架思路,尤其感谢姚芬芳等多位同学，是你们的帮助让我顺利完成论文的完善和校对。

再一次对给予我帮助的老师、同学们表示最诚挚的感谢！

附录

附录1：

a=1; % 设置两条平行线之间的距离

l=0.6; % 投针的长度

counter=0; % 针与平行线相交的次数

n=10000000; % 投掷次数

x=unifrnd(0,a/2,1,n); %产生n个(0,a/2)之间均匀分布的随机数，这里

a/2是投针的中点到最近的平行线的距离phi=unifrnd(0,pi,1,n); % 产生n个(0,pi)之间均匀分布的随机数，这里

pi是投针到最近的平行线的角度for i=1:n

if x(i)

end

end

frequency=counter/n; % 计算相交的频率，即相交次数比总次数

Pi=2*l/(a*frequency) % 从相交的频率总求的pi

%运行结果