概率论实验报告蒙特卡洛方法估计积分值

概率论实验报告

——蒙特卡洛方法估计积分值

姓名：

学号：

班级：

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

1用蒙特卡洛方法估计积分 20sin x xdx π

?，2-0x e dx +∞?和

22221x y x y e dxdy ++≤??的值，并将估

计值与真值进行比较。

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

0x e dx ?

和

22x y +≤??的值，

并对误差进行估计。

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

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

（3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。 目的：（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；

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

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

实验一、估计2

sin x xdx π

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

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)

运行结果：

b=1.0026

b=1.0061

b=1.0037

b=1.0135

b=0.9932

b=0.9988

b=1.0213

b=1.0310

b=0.9813

b=1.0041

p =

1.0056

z =

1

f=0.000207

>> （运行截图）

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

实验二、估计

2-0x e dx +∞

?的值，并将估计值与真值进行比较。 I=dx e x ?+∞-02=1/2*pi dx e

pi e x x *2***2/1*2/2/22-+∞

∞--? =)(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)

运行结果：

b=0.8779

b=0.8650

b=0.8826

b=0.8551

b=0.8855

b=0.8823

b=0.8771

b=0.8641

b=0.9186

b=0.8740

p =

0.8782

z =

0.8862

f=0.000329

>> （运行截图）结果显示估计结果与真实值的方差为f=0.00329，估计结果与真实值非常接近。

实验三、估计

2222

1x y x y e dxdy ++≤??的值，并将估计值与真值进行比较。

MATLAB 代码：

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

运行结果：

I =

5.4006

d =

0.0012

I =

5.4069

d =

9.5383e-004

>> （运行截图）

实验四、估计2

1

0x e dx ?的值，并对误差进行估计。

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

x e ,)(x f x =1.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,1,n,1);

x=sort(a); y=exp(x(i).^2); s=s+y;

end

b=s./n;

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

h(j)=b;

s=0; m=m+b;

end

p=m./10

for j=1:10

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

end

f=f./9;

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

运行结果：

>> clear

b=1.4546

b=1.4723

b=1.4540

b=1.4584

b=1.4663

b=1.4341

b=1.4695

b=1.4314

b=1.4605

b=1.4938

p =

1.4595

f=0.000331

>> （运行截图）结果显示，误差为0.000331，以平均值作为真实值，均方误差也比较小。

实验五、估计22x y +≤??的值，并对误差进行估计。

MATLAB 代码：

n=1000;m=100;sum=0;S=0;I=0;

x=unifrnd(-2,2,m,n);y=unifrnd(-2,2,m,n); for j=1:m

s=0;

for i=1:n

if x(j,i)^2+y(j,i)^2<=4

s=s+16/sqrt(1+x(j,i)^4+y(j,i)^2); end

end

S(j)=s/n;sum=sum+S(j);

end

I=sum/m; D=0;d=0;

for j=1:m

D=D+(I-S(j))^2;

end

d=D/(m-1)

I

运行结果：

d =

0.0284

I =

7.4566

>> （运行截图）