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)