matlab--蒙特卡洛法估计积分值

matlab--蒙特卡洛法估计积分值

西安交通大学实验报告

课程：概率论与数理统计

实验日期：

报告日期：

专业班级：

姓名：学号：

实验内容：用蒙特卡洛方法估计积分值

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

（2）利用计算机产生所选分布的随机数以估计积分值；

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

目的：（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；

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

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

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

2

sin x xdx π

⎰，

2

x e

dx +∞

⎰和

22

221

x y x y e

dxdy

++≤⎰⎰

的值，并将估计值与真值进行比较。

1)dx x x ⎰20

sin π

仍是用均匀分布来估计此积分的大小，g(x)=xsinx,)(x f x

=1/(2

pi ).x>0.分别取10个估计值h(j)，求得估计值的均值p ，对照积分的真实值求得估计均方误差f 。 Matlab 程序代码如下： s=0;m=0;f=0;r=0;n=50; h(1:10)=0; for j=1:10 for i=1:n

a=unifrnd(0,pi/2,n,1); x=sort(a); y=pi/2*mean(x.*sin(x));

s=s+y; end b=s./n;

fprintf('b=%.4f\n',b);

h(j)=b;

s=0;

m=m+b;

end

p=m./10

z=1

for j=1:10

r=(h(j)-z).^2; f=f+r; end

f=f./10;

fprintf('f=%.6f\n',f)

结果显示f=0.000221，表明估计结果与理论值非常接近。

2)

dx

e x ⎰+∞-

2

I=

dx

e x

⎰+∞-

2

=1/2*

pi

dx

e

pi

e x

x*2

*

*

*2

/1*2/

2/2

2-

+∞

∞

-

-

⎰

=)(x

f x

2/2

*

*2

/1x e

pi- g(x)=e pi

x*

2

*2/2-

)(x

f

x为标准正态分布的概率密度.分别取10个估计值h(j)，求得估计值的均值p，对照积分的真实值求得估计均方误差f。

Matlab程序代码如下：

s=0;m=0;f=0;n=50;r=0;

h(1:10)=0;

for j=1:10

for i=1:n

a=normrnd(0,1,1,n);

x=sort(a);

z=(sqrt(2.*pi)).*exp(-x(i).^2./2); s=s+z;

end

b=(s./n)./2; fprintf('b=%.4f\n',b);

h(j)=b; s=0; m=m+b;

end

p=m./10

z=sqrt(pi)./2

for j=1:10

r=(h(j)-z).^2;

f=f+r;

end

f=f./10;

fprintf('f=%.6f\n',f)

结果如下：

结果显示估计结果与真实值的方差为f=0.00322，估计结果与真实值非常接近。

3)

22

221

x y

x y

e dxdy

+

+≤

⎰⎰

m=10000;sum=0;n=50;D=0;

X=unifrnd(-1,1,n,m);Y=unifrnd(-1,1,n,m); for i=1:n a=0; for j=1:m

if(X(i,j)^2+Y(i,j)^2<=1)

Z(i,j)=exp(X(i,j)^2+Y(i,j)^2); a=a+Z(i,j); end end

S(i)=a/m;sum=sum+S(i); end

I=sum/n*4 for i=1:n

D=D+(S(i)*4-pi*(exp(1)-1))^2; end

d=D/n

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

1

0x e

dx

⎰和224

4

1x y x y

+≤++⎰⎰

的值，并对误差进行估计。 1)dx e x ⎰1

02

此积分采用的是均匀分布。g(x)=2

x

e ,)(x

f x

=1.x>0.