用蒙特卡洛方法估计积分方法及matlab编程实现.docx
《用蒙特卡洛方法估计积分方法及matlab编程实现.docx》由会员分享,可在线阅读,更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现.docx(21页珍藏版)》请在冰豆网上搜索。
用蒙特卡洛方法估计积分方法及matlab编程实现
用蒙特卡洛方法估计积分
方法及matlab编程实现
专业班级:
材料43
学生姓名:
王宏辉
学号:
2140201060
指导教师:
李耀武
完成时间:
2016年6月8日
用蒙特卡洛方法估计积分
方法及matlab编程实现
实验内容:
1用蒙特卡洛方法估计积分
,
和
的值,并将估计值与真值进行比拟。
2用蒙特卡洛方法估计积分
和
的值,并对误差进行估计。
要求:
〔1〕针对要估计的积分选择适当的概率分布设计蒙特卡洛方法;
〔2〕利用计算机产生所选分布的随机数以估计积分值;
〔3〕进行重复试验,通过计算样本均值以评价估计的无偏性;通过计算均方误差〔针对第1类题〕或样本方差〔针对第2类题〕以评价估计结果的精度。
目的:
〔1〕能通过MATLAB或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等;
〔2〕熟练使用MATLAB对样本进行根本统计,从而获取数据的根本信息;
〔3〕能用MATLAB熟练进行样本的一元回归分析。
实验原理:
蒙特卡洛方法估计积分值,总的思想是将积分改写为某个随机变量的数学期望,借助相应的随机数,利用样本均值估计数学期望,从而估计相应的积分值。
具体操作如下:
一般地,积分
改写成
的形式,〔其中为
一随机变量X的概率密度函数,且
的支持域
),
〕;令Y=h(X),那么积分S=E〔Y〕;利用matlab软件,编程产生随机变量X的随机数,在由
,得到随机变量Y的随机数,求出样本均值,以此估计积分值。
积分
的求法与上述方法类似,在此不赘述。
概率密度函数的选取:
一重积分,由于要求
的支持域
,为使方法普遍适用,考虑到标准正态分布概率密度函数
支持域为
,应选用
。
类似的,二重积分选用
,支持域为
。
估计评价:
进行重复试验,通过计算样本均值以评价估计的无偏性;通过计算均方误〔针对第1类题,积得出〕或样本方差〔针对第2类题,积不出〕以评价估计结果的精度。
程序设计:
依据问题分四类:
第一类一重积分;第一类二重积分;第二类一重积分,第二类二重积分,相应程序设计成四类。
为了使程序具有一般性以及方便以后使用:
一重积分,程序保存为一个.m文本,被积函数,积分区间均采用键盘输入;二重积分,程序主体保存为一个.m文本,被积函数键盘输入,示性函数用function语句构造,求不同区域二重积分,只需改变function函数内容。
编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。
程序代码及运行结果:
第一类一重积分程序代码:
%%%构造示性函数
functionI=I1(x,a,b)
ifx>=a&&x<=b
I=1;
else
I=0;
end
%%%%%%%%%%%%%%%%
%%第一类一重积分,程序主体:
functionoutf11=f11()
g1=input('输入一元被积函数如x.*sin(x):
','s')%输入被积函数
g1=inline(g1);
a=input('输入积分下界a:
');%输入积分上下限
b=input('输入积分上界b:
');
Real=input('积分真值:
');%输入积分真值
fprintf('输入样本容量10^V1--10^V2:
\r')
V=zeros(1,2);
V
(1)=input('V1:
');%输入样本容量
V
(2)=input('V2:
');
form=V
(1):
V
(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
10
x=randn(1,n);
fori=1:
n
t1(i)=I1(x(i),a,b);%示性及求和向量
end
y1=g1(x)*((pi*2)^0.5).*exp(x.^2/2);
Y1(j)=y1*t1'/n;%单次实验样本均值
end
t=ones(1,10);
EY=Y1*t'/10;%十次均值
D=abs(EY-Real);%绝对误差
RD=D/Real;%绝对误差
d=0;
fori=1:
10
d=d+(Y1(i)-Real)^2;
end
d=d/(10-1);
EY1(m-V
(1)+1)=EY;%样本容量为10^m时的样本均值
D1(m-V
(1)+1)=D;%绝对误差
RD1(m-V
(1)+1)=RD;%绝对误差
MSE1(m-V
(1)+1)=d;%方差
end
Real,EY1,D1,RD1,MSE1
outf11=[EY1;D1;RD1;MSE1];%存放样本数字特征
运行结果:
%估计积分
,积分真值为1
m=f11
输入一元被积函数如x.*sin(x):
x.*sin(x)
g1=
x.*sin(x)
输入积分下界a:
0
输入积分上界b:
pi/2
积分真值:
1
输入样本容量10^V1--10^V2:
V1:
1
V2:
5
n=
10
n=
100
n=
1000
n=
10000
n=
100000
Real=
1
EY1=
D1=
RD1=
MSE1=
m=
%估计积分
M=f11
输入一元被积函数如x.*sin(x):
exp(-x.^2)
g1=
exp(-x.^2)
输入积分下界a:
0
输入积分上界b:
+inf
输入样本容量10^V1--10^V2:
V1:
1
V2:
4
n=
10
n=
100
n=
1000
n=
10000
Real=
EY1=
D1=
RD1=
MSE1=
M=
第一类二重积分程序代码:
%%%构造示性函数,求不同区域上积分只需更改示性函数
functionI=I2(x,y)
ifx^2+y^2<=1
I=1;
else
I=0;
end
%第一类二重积分程序主体
functionoutf12=f12()
g2=input('输入二元被积函数如exp(x.^2+y.^2):
','s')%输入被积函数
g2=inline(g2,'x','y');
Real=input('积分真值:
');%输入积分真值
fprintf('输入样本容量10^V1*10^V1--10^V2*10^V2:
\r')
V=zeros(1,2);
V
(1)=input('V1:
');%输入样本容量
V
(2)=input('V2:
');
form=V
(1):
V
(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
10
x=randn(1,n);
y=randn(1,n);
fori=1:
n
t2(i)=I2(x(i),y(i));%示性及求和向量
end
y2=g2(x,y)*(2*pi).*exp((x.^2+y.^2)/2);
Y2(j)=y2*t2'/n;%单次实验样本均值
end
t=ones(1,10);
EY=Y2*t'/10;%十次均值
D=abs(EY-Real);%绝对误差
RD=D/Real;%绝对误差
d=0;
fori=1:
10
d=d+(Y2(i)-Real)^2;
end
d=d/(10-1);
EY2(m-V
(1)+1)=EY;%样本容量为10^m时的样本均值
D2(m-V
(1)+1)=D;%绝对误差
RD2(m-V
(1)+1)=RD;%绝对误差
MSE2(m-V
(1)+1)=d;%方差
end
Real,EY2,D2,RD2,MSE2
outf12=[EY2;D2;RD2;MSE2];%存放样本数字特征
运行结果:
%估计积分
,真值为
m=f12
输入二元被积函数如exp(x.^2+y.^2):
exp(x.^2+y.^2)
g2=
exp(x.^2+y.^2)
输入样本容量10^V1*10^V1--10^V2*10^V2:
V1:
1
V2:
4
n=
10
n=
100
n=
1000
n=
10000
Real=
EY2=
D2=
RD2=
MSE2=
m=
第二类一重积分程序代码:
%%%构造示性函数
functionI=I1(x,a,b)
ifx>=a&&x<=b
I=1;
else
I=0;
end
%第二类一重积分程序主体
%
functionoutf21=f21()
g1=input('输入一元被积函数如exp(x.^2):
','s')%输入被积函数
g1=inline(g1);
a=input('输入积分下界a:
');%输入积分上下限
b=input('输入积分上界b:
');
fprintf('输入样本容量10^V1--10^V2:
\r')
V=zeros(1,2);
V
(1)=input('V1:
');%输入样本容量
V
(2)=input('V2:
');
form=V
(1):
V
(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
10
x=randn(1,n);
fori=1:
n
t1(i)=I1(x(i),a,b);%示性及求和向量
end
y1=g1(x)*((pi*2)^0.5).*exp(x.^2/2);
Y1(j)=y1*t1'/n;%单次实验样本均值
end
t=ones(1,10);
EY=Y1*t'/10;%十次均值
d=0;
fori=1:
10
d=d+(Y1(i)-EY)^2;
end
d=d/(10-1);
EY1(m-V
(1)+1)=EY;%样本容量为10^m时的样本均值
MSE1(m-V
(1)+1)=d;%方差
end
EY1,MSE1
outf21=[EY1;MSE1];%存放样本数字特征
%%%%
运行结果:
%估计积分
m=f21
输入一元被积函数如exp(x.^2):
exp(x.^2)
g1=
exp(x.^2)
输入积分下界a:
0
输入积分上界b:
1
输入样本容量10^V1--10^V2:
V1:
1
V2:
4
n=
10
n=
100
n=
1000
n=
10000
EY1=
MSE1=
m=
%用matlab指令求积分
f=inline('exp(x.^2)')
f=
Inlinefunction:
f(x)=exp(x.^2)
>>S=quadl(f,0,1)
S=
第二类二重积分程序代码:
%%%构造示性函数,求不同区域上积分只需更改示性函数
functionI=I2(x,y)
ifx^2+y^2<=1
I=1;
else
I=0;
end
%第二类二重积分函数主体
functionoutf22=f22()
g2=input('输入二元被积函数如1./(1+x.^4+y.^4).^0.5:
','s')%输入被积函数
g2=inline(g2,'x','y');
fprintf('输入样本容量10^V1*10^V1--10^V2*10^V2:
\r')
V=zeros(1,2);
V
(1)=input('V1:
');%输入样本容量
V
(2)=input('V2:
');
form=V
(1):
V
(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
10
x=randn(1,n);
y=randn(1,n);
fori=1:
n
t2(i)=I2(x(i),y(i));%示性及求和向量
end
y2=g2(x,y)*(2*pi).*exp((x.^2+y.^2)/2);
Y2(j)=y2*t2'/n;%单次实验样本均值
end
t=ones(1,10);
EY=Y2*t'/10;%十次均值
d=0;
fori=1:
10
d=d+(Y2(i)-EY)^2;
end
d=d/(10-1);
EY2(m-V
(1)+1)=EY;%样本容量为10^m时的样本均值
MSE2(m-V
(1)+1)=d;%方差
end
EY2,MSE2
outf22=[EY2;MSE2];%存放样本数字特征
运行结果:
%估计积分
m=f22
g2=
输入样本容量10^V1*10^V1--10^V2*10^V2:
V1:
1
V2:
4
n=
10
n=
100
n=
1000
n=
10000
EY2=
MSE2=
m=
实验结果整理:
第一类一重积分:
估计积分
积分真值:
1积分估计值:
样本容量:
10100100010000100000
样本均值:
绝对误差:
相对误差:
均方误差:
估计积分
积分估计
样本容量:
10100100010000
第一类二重积分:
估计积分
积分估计值:
样本容量:
10100100010000
样本均值:
绝对误差:
相对误差:
均方误差:
第二类一重积分:
估计积分
积分估计
样本容量:
10100100010000
第二类二重积分:
估计积分
积分估计
样本容量:
10100100010000
实验结果分析:
从第一类积分看,以估计积分
为例:
积分真值:
1积分估计值:
样本容量:
10100100010000100000
样本均值:
绝对误差:
相对误差:
均方误差:
随着样本容量的增大,样本均值有接近积分真值的趋势,绝对误差、相对误差、均方误差呈减小趋势;随着样本容量的增大,样本均值有接近积分真值的趋势,说明估计具有无偏性;绝对误差、相对误差、均方误差呈减小趋势,说明增大样本容量能提高估计精度;验证了蒙特卡洛方法估计积分值的可行性,为后续估计第二类积分提供了参考。
从第二类积分看,以估计积分
为例:
样本容量:
10100100010000
由于积分真值未知,无法直接比拟估计值与积分值值;但随样本容量增大,样本方差减小,间接反映了估计精度的提高。
蒙特卡洛方法估计值1.4590相比用matlab指令求得的积分结果1.4627,绝对偏差0.0038,相对偏差0.0025。
蒙特卡洛方法估计值与用matlab指令求得的积分结果相互验证。
总结与讨论:
蒙特卡洛方法是基于随机数的一种统计方法。
蒙特卡洛方法估计积分值,总的思想是将积分改写为某个随机变量的数学期望,借助相应的随机数,利用样本均值估计数学期望,从而估计相应的积分值。
为使方法具有一般性,概率密度函数一重积分选择了
,二重积分选用
。
程序设计方面,本着使程序具有一般性以及方便以后使用的原那么,依据问题分四类:
第一类一重积分;第一类二重积分;第二类一重积分,第二类二重积分,相应程序设计成四类,并存储为.m文件,用蒙特卡洛方法估计积分值,一重积分只需调用相应程序即可;二重积分只需依据积分域修改相应示性函数即可调用相应函数求解。
极大方便了同类问题求解。
实验运行结果说明本方案可操作性良好。
遗留问题:
本次实验未设计选用不同概率密度函数,估计精度的比拟,留有不同条件下选用何种概率密度函数估计效果最正确?
如何缩短程序运行时间?
如何对程序进行封装?
如何更好评价第二类积分估计值无偏性以及精度?
等问题。
姓名:
王宏辉
班级:
材料43
学号:
2140201060