matlab蒙特卡洛法估计积分值.docx
《matlab蒙特卡洛法估计积分值.docx》由会员分享，可在线阅读，更多相关《matlab蒙特卡洛法估计积分值.docx（9页珍藏版）》请在冰点文库上搜索。
matlab蒙特卡洛法估计积分值
matlab--蒙特卡洛法估计积分值
西安交通大学实验报告
课程：
概率论与数理统计
实验日期：
报告日期：
专业班级：
姓名：
学号：
实验内容：
用蒙特卡洛方法估计积分值
要求：
（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；
（2）利用计算机产生所选分布的随机数以估计积分值；
（3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。
目的：
（1）能通过MATLAB或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；
（2）熟练使用MATLAB对样本进行基本统计，从而获取数据的基本信息；
（3）能用MATLAB熟练进行样本的一元回归分析。
1用蒙特卡洛方法估计积分
，
和
的值，并将估计值与真值进行比较。
1）
仍是用均匀分布来估计此积分的大小，g（x）=xsinx,
=1/（
）.x>0.分别取10个估计值h（j），求得估计值的均值p，对照积分的真实值求得估计均方误差f。
Matlab程序代码如下：
s=0;m=0;f=0;r=0;n=50;
h（1:
10）=0;
forj=1:
10
fori=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
forj=1:
10
r=（h（j）-z）.^2;f=f+r;
end
f=f./10;
fprintf（'f=%.6f\n',f）
结果显示f=0.000221，表明估计结果与理论值非常接近。
2）
I=
=1/2*
g（x）=e
为标准正态分布的概率密度.分别取10个估计值h（j），求得估计值的均值p，对照积分的真实值求得估计均方误差f。
Matlab程序代码如下：
s=0;m=0;f=0;n=50;r=0;
h（1:
10）=0;
forj=1:
10
fori=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
forj=1:
10
r=（h（j）-z）.^2;
f=f+r;
end
f=f./10;
fprintf（'f=%.6f\n',f）
结果如下：
结果显示估计结果与真实值的方差为f=0.00322，估计结果与真实值非常接近。
3）
m=10000;sum=0;n=50;D=0;
X=unifrnd（-1,1,n,m）;Y=unifrnd（-1,1,n,m）;
fori=1:
n
a=0;
forj=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
fori=1:
n
D=D+（S（i）*4-pi*（exp
（1）-1））^2;
end
d=D/n
2用蒙特卡洛方法估计积分
和
的值，并对误差进行估计。
1）
此积分采用的是均匀分布。
g（x）=
=1.x>0.分别取10个估计值h（j），求得估计值的均值p，对照积分的真实值求得估计均方误差f。
Matlab程序代码如下：
s=0;m=0;f=0;r=0;n=50;
h（1:
10）=0;
forj=1:
10
fori=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
forj=1:
10
r=（h（j）-p）.^2;f=f+r;
end
f=f./9;
fprintf（'f=%.6f\n',f）
结果如下：
结果显示，误差为0.000322，以平均值作为真实值，均方误差也比较小。
2）
n=1000;m=100;sum=0;S=0;I=0;
x=unifrnd（-2,2,m,n）;y=unifrnd（-2,2,m,n）;
forj=1:
m
s=0;
fori=1:
n
ifx（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;
forj=1:
m
D=D+（I-S（j））^2;
end
d=D/（m-1）
升级会员 