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

概率论实验报告
——蒙特卡洛方法估计积分值
姓名：
学号：
班级：
实验内容：用蒙特卡洛方法估计积分值
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，以平均值作为真实值，均方误差也比较小。

实验五、估计
22
44 1
1
1
x y
dxdy x 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
>> （运行截图）。