用蒙特卡洛方法估计积分方法及matlab编程实现

用蒙特卡洛方法估计积分方法及matlab编程实现

专业班级：材料43

学生姓名：王宏辉

学号：2140201060

指导教师：李耀武

完成时间：2016年6月8日

用蒙特卡洛方法估计积分 方法及matlab 编程实现

实验内容：

1用蒙特卡洛方法估计积分 2

0sin x xdx π?，2

-0x e dx +∞

?和

2

2

221

x

y x y e dxdy ++≤??

的值，

并将估计值与真值进行比较。 2用蒙特卡洛方法估计积分 21

0x e dx ?

和

22x y +≤??

的值，

并对误差进行估计。

要求：

（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；

（2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针

对第1类题）或样本方差（针对第2类

题）以评价估计结果的精度。 目的：

（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数

及其期望、方差、协方差等；

（2） 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息；

（3） 能用 MATLAB 熟练进行样本的一元回归分析。

实验原理：

蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。 具体操作如下：

一般地，积分?=b

dx x g a )(S 改写成??==b

b dx f h dx f g a

a )(x )(x )(x f(x))

(x S 的形

式，（其中为)f(x 一随机变量X 的概率密度函数，且)f(x 的支持域

）（｝｛b f ,a 0)(x |x ?>)，f(x )

)

(x )(x g h =

）;令Y=h(X)，则积分S=E （Y ）；利用matlab 软件，编程产生随机变量X 的随机数，在由

??

??∈==)

b (a,,)

b (a,,01I(x) ,)(x )(x y x x I h ，得到随机变量Y 的随机数，求出样本均值，以此估计积分值。

积分??=A

dxdy g S )y (x,的求法与上述方法类似，在此不赘述。

概率密度函数的选取：

一重积分，由于要求)f(x 的支持域）（｝｛b f ,a 0)(x |x ?>，为使方法普

遍适用，考虑到标准正态分布概率密度函数2

2

e 21)(x x

f -=π

支持域为

R ，故选用2

2e 21)(x x f -

=

π

。

类似的，二重积分选用2

2

221)y (x,y x e f +-=π

，支持域为2R 。

估计评价：

进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度。

程序设计：

依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。

为了使程序具有一般性以及方便以后使用：一重积分，程序保存为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存为一个.m文本，被积函数键盘输入，示性函数用function 语句构造，求不同区域二重积分，只需改变function 函数内容。

编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。

程序代码及运行结果：

第一类一重积分程序代码：

%%%构造示性函数

function I=I1(x,a,b)

if x>=a&&x<=b

I=1;

else

I=0;

%保存为I1.m

%%%%%%%%%%%%%%%%

%%第一类一重积分,程序主体：

%保存为f11.m

function outf11=f11()

g1=input('输入一元被积函数如x.*sin(x):','s')%输入被积函数

g1=inline(g1);

a=input('输入积分下界a:');%输入积分上下限

b=input('输入积分上界b:');

Real=input('积分真值:');%输入积分真值

fprintf('输入样本容量10^V1--10^V2:\r')

V=zeros(1,2);

V(1)=input('V1:');%输入样本容量

V(2)=input('V2:');

for m=V(1):V(2)%样本容量10^m1--10^m2

n=10^m

for j=1:10

x=randn(1,n);

for i=1:n

t1(i)=I1(x(i),a,b);%示性及求和向量

y1=g1(x)*((pi*2)^0.5).*exp(x.^2/2);

Y1(j)=y1*t1'/n; %单次实验样本均值

end

t=ones(1,10);

EY=Y1*t'/10; %十次均值

D=abs(EY-Real); %绝对误差

RD=D/Real; %绝对误差

d=0;

for i=1:10

d=d+(Y1(i)-Real)^2;

end

d=d/(10-1);

EY1(m-V(1)+1)=EY; %样本容量为10^m时的样本均值

D1(m-V(1)+1)=D; %绝对误差

RD1(m-V(1)+1)=RD; %绝对误差

MSE1(m-V(1)+1)=d; %方差

end

Real,EY1,D1,RD1,MSE1

outf11=[EY1;D1;RD1;MSE1]; %存放样本数字特征

%保存为f11.m