用Monte Carlo方法计算定积分.docVIP

用MonteCarlo方法计算定积分(随即投点法)理论思路

设，求在区间[0，1]上的积分值:fx()0()1,,fx

1

Jfxdx,(),0

设(X,Y)服从正方形上的均匀分布，则可知X服从[0，1]上01,01,,,,xy,,

的均匀分布，Y也服从[0,1]上的均匀分布,且X与Y独立，又记事件

AYfx,,(),,

则A的概率为

1()1fxpPYfxdydxfxdxJ,,,,,(())(),,,000

即定积分的值J就是事件A的概率p。由伯努利大数定律，我们可以用重复试验中A出现的频率作为p的估计值。这种求定积分的方法也就是MonteCarlo随即投点法，即将(X,Y)看成是正方形内的随机投点，用随即01,01,,,,xy,,

点落在区域

中的频率作为定积分的近似值。{()}Yfx,

(1)先用计算机产生(0,1)上均匀分布的2n个随机数:。xyin,,1,2,...,,ii

(2)对n对数据(xyin,),1,2,...,,，记录满足如下不等式ii

yfx,()ii

uunnu的次数，这就是事件A发生的概率。由此可得事件A发生的频率,则。J,nnn用MonteCarlo方法计算定积分(平均值法)理论思路

为计算定积分:

1Jfxdx,(),0

设随即变量X服从(0,1)上的均匀分布，则YfX,()的数学期望为

1EfXfxdxJ,,()，，，，,0

所以估计J的值就是估计f(X)的数学期望的值。有辛钦大数定律，可以用f(X)的观察值的平均去估计f(X)的数学期望的值。

具体做法如下:

先用计算机产生n个(0,1)上均匀分布的随机数:，然后对xin,1,2,...,,i

n1每个计算，最后的J的估计值为。xJfx,fx()(),iiin,1i

一、MonteCarlo在一重积分上的应用

考虑一个简单的一重定积分的应用:

b1,,fxdx()，，,a

不少统计问题，如计算概率、各阶矩等，最后都归结为定积分的近似计算问题，我们首先介绍两种求的简单的MonteCarlo方法，并给出求积分的几个实例。随机投点法:

考虑(1)的积分，简单起见，设a、b有限，y=f(x)在[a,b]上连续非负，有

max其中。实际上，如果对于函数f(x)，不是在区间[a,b]0(),,fxMMfx,()axb,,

上都有，则进行计算时可以利用恒等式:fx()0,

bb=，其中h,0是适当选择的正数，使对一切fxdx()[(())()]fxhdxhba，,,,,aa

，都有。fxh()0，,xxb:a,,

图3

,M令，}，并设(X，Y)是在上均匀分布的二,,,,,,{(,):,0}xyaxbyM

b1,维随机变量，联合密度函数为则易见是中曲,,fxdx()I(,0)axbyM,,,,,aMba,()

,线下方的面积(如上图所示)。随机投点法的思想是:向中进行随机yfx,()

投点。若点落在下方称为中的，否则成为不中，则点中的概率为:yfx,()

,p,，若进行了N次投点，其中n次中的，则得到的一个估计,Mba(),

n。,,,Mba()1N

其实施步骤为:

1)独立产生2N个U(0,1)的随机数uviN,,1,2,...,;,ii

2)计算和;fx()xaubayMv,，,,(),iiiii

3)统计个数n;fxy(),ii

n4)计算θ估计值。,,,Mba()1N

,sinx例[1]:用MonteCarlo方法求积分值dx,0x

sinxmax解:用公式(1),其中,去N=20，利用Matlab编程,,,,,abM0,,1axb,,x相应程序如下:

s=0;

a=unifrnd(0,pi,1,20);

x=sort(a)

fori=1:19

y=x(i+1)-x(i)

z=sin(x(i))/x(i)

s=s+y*z;

end

s

Matlab编程可得到结果:

sinxifx(),ixii

10.7225

20.3345

30.2423

40.2805

50.9052

60.8992

70.7201

80.403

90.3961

100.9431

110.8474

120.4438

130.4385

140.6692

150.9327

160.2334

170.2148

180.9569

190.4379

200.5157

y则满足条件yfx,()的个数n=12,利用公式(1)得到，所求积分的估计值为iii

n12,,,,,,Mba()1.885N20

sinx由数学分析的知识可知，例2中的被积函数原函数是不能用初等函数表x,sinx