蒙特卡罗随机模拟投点法.docx

蒙特卡罗随机模拟投点法.docx

《蒙特卡罗随机模拟投点法.docx》由会员分享，可在线阅读，更多相关《蒙特卡罗随机模拟投点法.docx（18页珍藏版）》请在冰豆网上搜索。

蒙特卡罗随机模拟投点法

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

应用

数学与应用数学0901班：

张瑞宸

指导老师：

***

摘要：

本文首先介绍了蒙特卡罗方法的产生和发展，然后分析了蒙特卡罗方法计算数值积分的理论原理，最后给出了蒙特卡罗方法计算数值积分的MATLAB编程实现，全文主要是讨论了蒙特卡罗方法在定积分计算的应用。

而蒙特卡罗的优点：

可以计算被积函数非常复杂的定积分、重积分，并且维数没有限制，这是别的数值积分方法还未达到的。

蒙特卡罗的缺点：

收敛速度慢，误差一般较大，且是概率的误差，不是真正的误差。

关键词：

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

Abstract：

ThispaperfirstintroducestheemergeneeanddevelopmentoftheMonteCarlomethod,andthenanalyzethetheoreticalprinciplesofMonteCarlonumericalintegrationmethod,Full-textmainlydiscussedtheapplicationoftheMonteCarlomethodinthedefiniteintegral.TheadvantagesofMonteCarlo:

canbecalculatedtheintegrablefunctionsverycomplexdefiniteintegral,Multipleintegrals,anddimensionnolimit,othernumericalintegrationmethodshavenotyetreached.MonteCarloDisadvantages:

slowconvergeneespeed,errorgenerallyhigher,andtheprobabilityoferror,notarealerror.

Keywords:

MonteCarlomethod,Meanestimationmethod,numericalintegral,Matlabprogramming

0引言

历史上有记载的蒙特卡罗试验始于十八世纪末期（约1777年），当时布丰

（Buffon）为了计算圆周率，设计了一个“投针试验”，后文会给出。

虽然方

法已经存在了200多年，此方法命名为蒙特卡罗则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用蒙特卡罗方法模拟中子对某种特殊材料的穿透作用。

出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一vonNeumann灵犀一点选择摩洛哥著名赌城蒙特卡罗（MonteCarlo）作为

该项目名称，自此这种方法也就被命名为MonteCarlo方法广为流传。

蒙特卡罗方法，又名随机模拟法或统计实验法它是以概率统计理论为基

础，依据大数定律（样本均值替代总体均值）利用电子计算机数字模拟技术，解

决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似

计算法。

本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。

通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。

对于那些由于

计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一

种有效的求出数值解的方法。

一般蒙特卡罗法在数学中最常见的应用就是蒙特卡罗随机投点和蒙特卡罗数值积分。

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

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

1.1蒙特卡罗方法的产生

作为蒙特卡罗方法的最初应用，是解决蒲丰氏问题1777年，法国数学家

Bufon提出利用投针实验求解的问题：

设平面上有无数多条距离为1的等距平行线，现向该平面随机的投掷一根长度为ll<1的针。

随机投针是指针的中心点于最近的平行线间的距离x均匀分

布在［0,1/2］上，针与平行线的夹角:

（不管相交与否）均匀分布在〔0,二1上，如

1

图一所示，故而，故其x~U（0,-），「~U（0,二）概率密度函数分别为

2

11

P（x）,P（：

）：

2兀

xl

故我们得到针与线相交的充要条件是——--

sin屮2

则针与线相交的概率是px»「）

2i

丁

假如我们能做大量的投针实验并记录下针与线的相交次数，则可以根据大数

定律估计出针线相交的概率P。

投针实验N次可能有n次使针与任意平行线相交那么p：

、巴，显然，试验次数N越多，P的近似程度越好。

历史上曾有几位学者

N

做过这样的投针试验，并用手工计算出二值，结果参见表1

表1

实验者

时间（年）

针长1

投针次数

相交次数

兀估计值

Wolf

1850

0.80

5000

2532

3.15956

Smith

1855

0.60

3024

1218

3.15665

Fox

1884

0.75

1030

489

3.15951

Lazzarini

1925

0.83

3408

1808

3.14159292

应该指出，上述试验的精度一般不会很高。

譬如，假设丨「，则p=Z。

贝u由

3T

中心极限定理，如果试验次数为N,则p的估计值0渐进服从N（p,p（1—p）），近

N

似为N（0.6366,0.2313）。

因此，若要以95%勺概率保证p的精确到三位有效数字，

N

即p与p的差距小于0.001，则N必须满足N=1.962x0.2313/0.001^8.8^105。

重复进行上千次的投针实验和手工计算，要消耗大量的人力、财力，蒙特卡罗方法虽然能解决此类问题，但得不到推广应用。

1.2蒙特卡罗方法的发展

20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。

计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算处理，于是MonteCarlo方法重新被提起，引起世人重视，应用日渐广泛。

实际上，采用MonteCarlo方法在计算机上建立模型来解决Buffon问题是非常简单的。

我们在计算机上进行模拟试验，给定I，我们可以在计算机上随机产生x和

'，然后判断」-是否成立。

若成立，则针线相交，否则不交。

假如我们在

sin®2

计算机上独立的产生N对这样的x和'，并记录下」L成立的次数,记为n<

sin®2

则二估计值可取为，这就是随机模拟计算的结果。

2

借助Matlab在计算机上进行投针实验，（相应程序见附录1）所取得的计算结果参见表2。

表2

N、、

I=0.4

I=0.6

1=0.8

I=0.9

I=1.0

50000

3.1252442

3.1583934

3.1392246

3.14642309

3.1464351

100000

3.1241457

3.1211978

3.1269543

3.1324504

3.1286664

200000

3.1414435

3.1396353

3.1373472

3.1366262

3.1371811

500000

3.1519392

3.1470787

3.1455443

3.1422826

3.1425492

800000

3.1483358

3.1465609

3.1440825

3.1423004

3.1430294

1000000

3.1435296

3.1459072

3.1412646

3.1419370

3.1423368

2000000

3.1497197

3.1462907

3.1418722

3.1410570

3.1407306

5000000

3.1393330

3.1398735

3.1410784

3.1416836

3.1412106

其中二=3.1416。

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

随着电子计算机的普及，蒙特卡罗方法作为一种独到的方法得到开发，并首先应用到核武器的试验与研制中。

尤其是各种可是编程方法的不断涌现，更显示出蒙特卡罗方法的最独到的优点，即形象直观地用数学方法在电子计算机上实现数字模拟实验。

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

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

随机模拟法的基本原理非常简单，下面给以直观的说明。

11

首先构造一个定积分4[Jl-x2dx=^，其中积分S=加1-x2dx（即1/4单

位圆的面积）用随机模拟方法求得。

图2.1中粗线是1/4单位圆，如果向图2.1中边长为1的正方形里随机投n块小石头，当n很大时小石头会大致均匀地分布在正方形中，数一下落在1/4单位圆内的小石头，假定有k个，那么k/n就可以看作1/4单位圆面积n/4的近似值。

小石头的位置坐标可以用产生均匀分布随机数的程序得到，记为

X,yi（i=1,2,...,n）是[0,1]区间均匀分布随机数（x$相互独立），记录满足

x2yi25i=1,2,...,n）的数量k，即得n=4k/n。

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

[2]

大数定律（伯努利（Bernoull）定理）设k是n次独立重复试验中事

件A发生的次数，p是事件A每次试验中发生的概率，则对任意的正数；，有

=1

（1）

若规定“向图中正方形随机投一块小石头落在四分之一单位圆里”为事件

A发生，则A发生的概率p应该等于四分之一单位圆面积，随机投n块石头就是

独立重复做n次试验，事件A发生k次，由

（1）式，n无限变大时k/n与p之差小于任意一个数；的概率趋于1

用计算机算一下就会发现，即使n很大结果也不好，并且很不稳定，远不如常规数值积分的几种方法。

实际上，随机模拟法很少用来做定积分，它的特点是能够方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。

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

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

（1）随机投点法

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

一单位圆面积，即

1初-X21

pdydx1「x2dx

（2）

n次试验由计算机完成，采用［0,1］区间上的均匀分布产生相互独立的随机数。

记这样产生的n个点的坐标为（Xi,yi）i=1,2，…,n）。

事件a发生等价于以』）满足yi—1f2，A发生的次数是满足yi「1-X：

（i=1,2,…,n）点的个数k。

由

伯努利定理，p可以用k/n近似替代

这种方法可以推广如下。

对任意区间［a,b］内的连续函数f（x），满足

b

C乞f（x）乞d（c,d—0），为计算定积分f（x）dx，先由计算机产生n个点的坐标

La

（Xi,yJ（i=1,2,...,n），其中人$分别为［a,b］和［c,d］区间上的均匀分布随机数，然

后记录n个点中满足yj_f（xj的数目k，则

（3）

bk

f（x）dx：

（b-a）（d-c）（b-a）c

an

（2）均值估计法

这种方法依据概率论的一下两个定理:

［2］

大数定理（辛钦定理）若随机变量¥,丫2,..・,£相互独立，服从同一个

分布，且具有数学期望EYf二（i=1,2,...,n）,则对任意的正数；，有

［2］

随机变量函数的期望若随机变量X的概率分布密度是

p（x）（a岂x空b），则随机变量的函数丫二f（X）的期望为

b

E（f（X））=Jf（x）p（x）dx（5）

a

于是，当X为［a,b］区间均匀分布的随机变量时，

p（x）=1/（b-a）（alx乞b），（5）式给出

b

af（x）dx二（b-a）E（f（x））（6）

只要产生［a,b］区间相互独立、均匀分布的随机数Xi（i=1,2,…,n），f（xj就相互独立。

根据辛钦定理，当n很大时期望EY二E（f（X））就可以用f（xj的平均值近似，所以由（6）式得到

b（b-a）n

ff（x）dx^无f（x）（7）

any

与随机投点法相比，均值估计法没有C乞f（x）岂d（c,d—0）的限制，只需计算f（xj，不需要产生随机数yi，也不需要作yi乞f（xj的比较，显然大为方便。

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

用随机模拟作数值积分的优点不仅在于计算简单，尤其是它可以方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。

而通

过上文对蒙特卡罗计算定积分的两种方法的比较，我们可知均值估计法比投点法

更好操作，经过计算还可知均值估计法精确度（若选择的密度函数适合）会更高。

重积分本文讨论均值估计法。

考虑k维积分：

\=…f（x1,x2,...,xk）dx1dx2...dxk

其中。

为k维积分区域，0匸V（V:

q兰Xj兰b，i=1,2,...,k）

选取门上的一个概率密度函数P（X1,X2，…,Xk），门V（V：

a

i=1,2,...,k），且p（XnX2,...,Xk）满足...p（Xi,X2,...,Xk）dx』X2...dxk=1，Q

（Xi，X2,…，Xk）门，故:

...f（Xi,X2,...,Xk）p（Xi,X2,...,Xk）dX.,dX2...dXk

I

Q'P（Xi，X2，...,Xk）

二E（f（Xi,X2,..X"）p（Xi,X2,..Xk）

其中，（Xi,X2，…，Xk）-''-1

1「f（Xi，X2,...Xk）

P（Xi，x2,...Xk）其中，（Xi，X2，...，Xk）11

（8）

这里m是落在门区域中的随机投点数。

对于二重积分：

\=f（x,y）dXdy

Q

其中「：

a空x^b，cmq（x）乞y乞g2（x）乞d。

记V=[a,b][c,d]，设s，Sv分

1

别是'1,V的测度，若选取g（x,y）是门上的均匀分布，则g（x,y）—,（x,yf'■10

由几何概率得S「m，即鱼，于是:

SVnmn

I=芈単g（x,y）dxdy门g（x,y）

=1、mf（x,y）

myg（x,y）

Smsm

二f（x,y）二鱼、f（x,y）

m^ny

m

综上可得，If（x,y）（9）

ni#

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

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

MATLAB提供的命令

unifrnd（a,b,m,n）

产生m行n列［a,b］区间上均匀分布随机数。

当a=0,b=1时，可用rand（m,n）。

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

解

（1）随机投点法。

MATLAB程序如下：

n=10000;

x=rand（2,n）;

k=0;fori=1:

n

ifx（1,i）A2+x（2,i）A2<=1

k=k+1;

end

end

p=4*k/n

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

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

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

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

解

（2）均值估计法，此时f（X）「1-X2,MATLAB程序如下：

n=50000;

x=rand（1,n）;

y=0;fori=1:

n

y=y+sqrt（1-x（i）A2）;

endp=4*y/n

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

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

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

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

可以看出它的缺点是计算量大，结果具有随机性，精度较低。

一般说来精度为

n4/2阶，门增加时精度提高较慢。

但是用随机模拟方法可以计算被积函数非常复

杂的定积分、重积分，并且维数没有限制。

如下介绍多重积分的MALTAB实现

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

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

functions=mtc（f,fai1,fai2,a,b,c,d,n）

ifnargin<8n=10000;end

x=unifrnd（a,b,1,n）;y=unifrnd（c,d,1,n）;

s=0;

fork=1:

n

ifand（y（k）>=feval（fai1,x（k））,y（k）v=feval（fai2,x（k）））

s=s+feval（f,x（k）,y（k））;

end

end

s=s/n*（b-a）*（d-c）;

return

functions=mtcx（f,fai1,fai2,a,b,c,d,n）

ifnargin<8n=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

三重或三重以上积分只须增加若干个输入参数，即新增积分变量的上、下限,

方法不变，对程序稍作修改即可。

~~22

1j丄七2

例2计算-匚亡e2sin（x+y）dydx

一2卄

积分的精确值为0.4119295461763Q在命令窗口调用函数M文件mtc.m和

mtcx.m情况如下：

>>tic;s1=mtc（inline（exp（-x.A2/2）.*sin（x42+y）',

inline（-sqrt（1-x.A2/2）''nline（sqrt（1-x.A2/2）^）,-0.5,1,-1,1,100000）,toc

s1=

0.40901034347183

Elapsedtimeis40.094000seconds.

>>tic;s2=mtcx（inline（exp（-x.A2/2）.*sin（x42+y）',

inline（-Sqrt（1-x.A2/2）',inline（Sqrt（1-x.A2/2）'）,-0.5,1,-1,1,100000）,toc

s2=

0.41183962563018

Elapsedtimeis0.047000seconds.

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

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

functions=mtcxg（f,fai1,fai2,a,b,c,d,n）

ifnargin<6n=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.A2/2）.*sin（x.A2+y）））,inline（‘-sqrt（1-x.A2/2）））,inline（‘sqrt（1-x.A2/2）））,-0.5,1,-1,1,100000）,tocs3=

0.41200550284008

Elapsedtimeis0.922000seconds.

对比常规数值积分和蒙特卡罗方法数值积分,同样数量的n值——也就意味这几乎相同的计算量——常规数值积分结果的精确度要高于蒙特卡罗数值积分的结果。

那么,我们为何还需要用蒙特卡罗来算数值积分呢？

答案的关键在于,常规数值积分的精度直接取决于每个维度上取点数量,维度增加了,但是每个维度上要取的点却不能减少。

在多重积分中，随着被积函数维度增加，需要计算的函数值数量以指数速度递增。

a

例如在一重积分.bf（x）「（x）dx中，只要沿着x轴取n个点；要达到相同大小的精确度，在k重积分11…f（X1X2…xj（XiX2...Xn）dxidx2...dXk中，仍然需要在每个维度上取n个点，k个纬度的坐标相组合，共需要计算nxk个坐标对应的f（x）函数值。

取点越多，会占用计算机大量内存，也需要更长运算时间，最终导致这种计算方法不可行！

蒙特卡罗方法却不同，不管是积分有多少重，取n个点计算的结果精确度都差不多。

因此，即使在一重积分的情形下，蒙特卡罗方法的效率比不过常规数值积分，但随着积分维度增加，常规数值积分的速度呈指数下降，蒙特卡罗方法的效率却基本不变。

经验表明，当积分重数达到4重积分甚至更高时，蒙特卡罗方法将远远优于常规数值积分方法。

4总结

归纳起来本论文主要论述了：

蒙特卡罗算法在数值积分中的运用。

本篇论文给出了蒙特卡罗方法应用在数值积分中的原理，讨论蒙特卡罗算法运用在数值积

分中的两种计算方法，随机投点法和均值估计法，具体蒙特卡罗算法的实现给出了MATLA编程的实例。

最后对蒙特卡罗算法的优缺点进行了讨论。

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

参考文献：

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

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

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

[4]李