用蒙特卡洛方法估计积分方法及matlab编程实现.docx
《用蒙特卡洛方法估计积分方法及matlab编程实现.docx》由会员分享，可在线阅读，更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现.docx（20页珍藏版）》请在冰豆网上搜索。
用蒙特卡洛方法估计积分方法及matlab编程实现
用蒙特卡洛方法估计积分
方法及matlab编程实现
专业班级:
材料43
学生姓名:
王宏辉
学号:
2140201060
指导教师:
李耀武
完成时间:
2016年6月8日
用蒙特卡洛方法估计积分
方法及ｍａｔlab编程实现
实验内容:
1用蒙特卡洛方法估计积分,与得值,并将估计值与真值进行比较。
２用蒙特卡洛方法估计积分与得值,并对误差进行估计。
要求:
（１）针对要估计得积分选择适当得概率分布设计蒙特卡洛方法;
（２）利用计算机产生所选分布得随机数以估计积分值;
（3）进行重复试验,通过计算样本均值以评价估计得无偏性；通过计算均方误差（针 对第1类题）或样本方差（针对第2类题）以评价估计结果得精度。
目得:
（1）能通过MＡＴLAB 或其她数学软件了解随机变量得概率密度、分布函数ﻩ及其期望、方差、协方差等;
（2） 熟练使用MＡTLAB对样本进行基本统计,从而获取数据得基本信息；
（3）能用MATLAB 熟练进行样本得一元回归分析。
实验原理：
蒙特卡洛方法估计积分值,总得思想就是将积分改写为某个随机变量得数学期望,借助相应得随机数,利用样本均值估计数学期望，从而估计相应得积分值。
具体操作如下:
一般地,积分改写成得形式,（其中为一随机变量X得概率密度函数,且得支持域）,）;令Y=ｈ（Ｘ）,则积分Ｓ=E（Y）;利用maｔlaｂ软件,编程产生随机变量X得随机数,在由,得到随机变量Y得随机数，求出样本均值,以此估计积分值。
积分得求法与上述方法类似，在此不赘述。
概率密度函数得选取：
一重积分,由于要求得支持域,为使方法普遍适用,考虑到标准正态分布概率密度函数支持域为,故选用。
类似得,二重积分选用，支持域为。
估计评价:
进行重复试验，通过计算样本均值以评价估计得无偏性;通过计算均方误（针对第１类题,积得出）或样本方差（针对第２类题,积不出）以评价估计结果得精度。
程序设计:
依据问题分四类:
第一类一重积分;第一类二重积分;第二类一重积分,第二类二重积分，相应程序设计成四类。
为了使程序具有一般性以及方便以后使用:
一重积分，程序保存为一个、m文本,被积函数,积分区间均采用键盘输入；二重积分,程序主体保存为一个、m文本,被积函数键盘输入,示性函数用functiｏn语句构造,求不同区域二重积分,只需改变functioｎ函数内容。
编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。
程序代码及运行结果:
第一类一重积分程序代码:
%%%构造示性函数
functｉｏnI=I1（x,a,ｂ）
ifx>=ａ&＆x＜=b
I=１；
ｅlｓe
I＝0;
eｎd
%保存为Ｉ１、ｍ
%％%％%%%%%％%％%%%%
%%第一类一重积分，程序主体：
%保存为f11、m
ｆｕｎction outf11=f１1（）
ｇ1=inｐut（＇输入一元被积函数如ｘ、*sin（x）:
'，'s'）%输入被积函数
g1＝ｉnline（ｇ1）;
a=ｉnpuｔ（'输入积分下界ａ:
'）;%输入积分上下限
ｂ＝ｉnput（'输入积分上界b:
'）;
Ｒｅal=iｎput（'积分真值:
'）;％输入积分真值
fprｉｎtf（＇输入样本容量 １０^V1－-10＾V2:
\r＇）
V=zeros（1,2）;
V
（1）=inpｕｔ（'V1:
'）；%输入样本容量
V
（2）=inpuｔ（＇V2:
＇）；
form=V
（1）：
V
（2）%样本容量10^m1－-１０＾m2
ｎ=１０^ｍ
ｆorj＝１:
10
x=raｎｄn（1,n）;
forｉ=１:
n
t1（i）=I1（x（ｉ）,ａ,b）;%示性及求与向量
eｎd
y1=g1（x）*（（pi＊２）^0、５）、*exｐ（x、^2/2）;
Y1（ｊ）=ｙ１*ｔ1'/n;%单次实验样本均值
enｄ
ｔ=onｅs（1，１０）;
EY=Y1*t＇/1０； %十次均值
Ｄ＝abs（EY-Real）; ％绝对误差
RD＝D/Reaｌ; %绝对误差
ｄ＝0;
fori=1：
10
d=ｄ+（Y1（i）-Real）^2;
end
d=d/（10-1）;
EY1（ｍ-V
（1）+1）=EＹ; %样本容量为１0＾m时得样本均值
Ｄ1（m-Ｖ（１）+1）=Ｄ； %绝对误差
RD1（m－V
（1）＋1）=ＲＤ；%绝对误差
MSE1（ｍ-V
（1）＋1）=ｄ;%方差
eｎd
Real,EY１，Ｄ1,RＤ1,MSＥ1
oｕtf11=[EY１；D1;ＲD1;MSE1］;%存放样本数字特征
%保存为f１1、m
运行结果：
%估计积分,积分真值为１
m=ｆ11
输入一元被积函数如ｘ、*ｓiｎ（x）:
x、*sin（x）
g1=
x、*sin（x）
输入积分下界ａ:
０
输入积分上界b：
pi／2
积分真值：
1
输入样本容量 1０＾V1--10^Ｖ2:
Ｖ1:
1
Ｖ2:
５
n=
1０
n=
100
n=
10００
ｎ =
1０000
n=
１0０0００
Real =
１
ＥY1=
1、263５ 1、0088 1、０066 １、010９1、0０1８
Ｄ１=
0、2635 ０、00８80、0０660、0１09 ０、０01８
RD1=
0、263５ ０、0０880、０0６6 0、0１0９ 0、0018
MSE1=
0、643９0、０2050、00２80、0006 0、0001
ｍ=
1、2635 1、0０8８ 1、０06６ 1、01０9 １、０018
0、2６350、00８８0、０066０、01090、００18
0、263５ 0、0０88 0、00６6 ０、0１09 ０、001８
0、643９ 0、0２05 0、０028 0、000６0、0001ﻩ
％估计积分真值为0、８8６2
Ｍ＝f11
输入一元被积函数如x、*ｓin（ｘ）:
exp（－x、^2）
g1=
exp（-ｘ、＾2）
输入积分下界a:
0
输入积分上界b:
+inf
积分真值：
pi^0、５／2%0、８86２
输入样本容量10^V１--１0＾V2:
V１:
1
V２:
4
n=
１0
ｎ=
100
ｎ＝
100０
n=
1000０
Reａl =
０、88６2
EY1 =
０、９33３0、9０770、88７3 0、887１
D1 =
0、０4700、０215 ０、0０10０、000９
RD1＝
0、0５３1０、０24３ ０、０012 0、０010
MSE1=
0、１927 ０、０1１20、０0１60、０000
Ｍ =
0、９３33 0、9０770、８873 0、8871
0、0４７０ 0、021５0、0010 ０、0０09
０、05３1 ０、０24３ ０、0012 0、0010
0、1９270、01１２0、001６ 0、0000
第一类二重积分程序代码：
%%%构造示性函数,求不同区域上积分只需更改示性函数
functionI＝I2（x,ｙ）
ifx^2+y＾2<＝1
I=1；
else
I＝0;
eｎd
%保存为I２、m
％第一类二重积分程序主体
%保存为ｆ12、m
functｉoｎoutf12＝f１２（）
ｇ２=ｉｎpｕt（'输入二元被积函数如exp（x、＾２+y、＾2）：
＇，'s'）％输入被积函数
g２=ｉｎliｎｅ（g2,'ｘ＇,'y＇）;
Ｒｅａl=iｎput（＇积分真值:
'）；%输入积分真值
fprｉnｔf（'输入样本容量10＾V1＊10^V1--1０^V2*１０^Ｖ2：
\ｒ'）
V=zeroｓ（1,2）；
V
（1）=iｎｐut（＇Ｖ１:
＇）;%输入样本容量
V
（2）＝input（'Ｖ2:
'）;
for m＝V
（1）:
V（２）%样本容量10＾m1-－10^ｍ２
n＝10＾ｍ
forj＝1:
10
x=raｎdｎ（1,n）;
ｙ＝ｒaｎdｎ（１,ｎ）;
foｒi=1:
n
t2（i）=I2（x（ｉ）,y（ｉ））;％示性及求与向量
eｎｄ
y2=ｇ２（x,y）*（２＊ｐi）、*exp（（x、＾2+y、^2）/2）;
Y2（j）=y2＊t２＇／n; %单次实验样本均值
eｎd
t=onｅs（1,10）;
ＥY=Y2＊t'／１0; %十次均值
Ｄ=ａbｓ（EY-Reaｌ）；%绝对误差
RD=Ｄ/Ｒeal; ％绝对误差
ｄ＝0;
foｒ i＝1:
１0
d=d+（Y2（i）-Ｒｅal）^2;
ｅnd
d＝d/（１0-1）;
EＹ2（m-V
（1）+1）=EＹ; %样本容量为10^m时得样本均值
D２（m-V（１）+1）=D; %绝对误差
RＤ2（ｍ-V
（1）＋1）=RD; %绝对误差
MＳＥ２（m-Ｖ
（1）＋1）=d;%方差
eｎd
Rｅaｌ,EY２,D2,RD2,MSＥ２
ｏutf12=[EY2；D2；RＤ2；MＳE2];%存放样本数字特征
％保存为ｆ１2、m
运行结果:
%估计积分,真值为pi*（exp
（1）-1）%5、3981
m=f１2
输入二元被积函数如exp（x、＾2+y、＾２）：
ｅxp（x、^2＋y、^2）
g２＝
eｘｐ（x、^2＋ｙ、^2）
积分真值:
pi＊（exp（１）-1）％５、3981
输入样本容量10＾V1*10^Ｖ１-－1０＾V2*10＾V2：
V1：
1
V2:
4
n=
１０
n=
１00
ｎ =
10００
ｎ=
10０0０
Reａl =
5、3９81
ＥＹ2 =
4、7７0２ 5、125０ 5、4317 ５、40４1
D2=
０、62７9 ０、2732 ０、0３35 0、0060
RＤ2=
0、1163 0、0５06 0、00620、００11
MＳE２ ＝
３、89６50、55640、02470、0017
m=
4、７7０2 5、1２５0５、4３175、4０41
０、6２79 0、27３2 0、0335 0、0060
0、１163０、０50６0、0０６2０、0011
3、8965 ０、55６4 0、02４７0、０0１7
第二类一重积分程序代码：
％%%构造示性函数
funcｔionI=I1（ｘ,a,b）
ｉf ｘ>＝ａ&&x<＝ｂ
I＝1；
else
I=０；
end
％保存为I1、ｍ
％第二类一重积分程序主体
％程序保存为f2１、m
ｆunctiｏnoutf２1＝f21（）
g１＝inｐuｔ（＇输入一元被积函数如exｐ（x、^2）:
＇,'s'）％输入被积函数
ｇ１=inｌｉne（g1）;
a=iｎpuｔ（＇输入积分下界a:
'）;％输入积分上下限
b=ｉnput（'输入积分上界ｂ:
＇）；
fprintf（＇输入样本容量10^V１－-１０^Ｖ２:
\r'）
Ｖ=ｚerｏｓ（1,２）；
V
（1）=iｎput（'V1:
＇）;%输入样本容量
Ｖ
（2）=input（'V2:
'）;
form＝V
（1）：
V
（2）%样本容量1０^m1－-１0^ｍ２
n=10＾ｍ
fｏrj=1:
10
x=raｎdn（1,n）；
for i=1：
n
t1（i）=Ｉ1（x（i），a，b）；％示性及求与向量
ｅnd
y１=g1（ｘ）＊（（pi*2）^０、5）、*exp（x、^2／２）;
Y1（j）=y１*t1'/n; %单次实验样本均值
end
t=ｏnｅs（１,10）;
EY=Y1*ｔ＇／１0;%十次均值
d=０;
fori＝1:
10
d=d＋（Y１（i）－ＥY）^2；
ｅｎｄ
ｄ=d/（10-1）；
EＹ1（ｍ-V
（1）+1）=EY; ％样本容量为10^ｍ时得样本均值
ＭSE1（ｍ-V（１）＋１）＝d; %方差
eｎd
EY1,ＭＳE１
ｏｕtf21=[EＹ1；MSＥ1];%存放样本数字特征
%%%%程序保存为f21、m
运行结果:
%估计积分
ｍ=ｆ21
输入一元被积函数如exp（x、^２）:
exp（x、^2）
g1 =
exｐ（x、^2）
输入积分下界ａ：
0
输入积分上界b:
１
输入样本容量１0^V1－-10^V2:
V1:
１
Ｖ2:
４
n＝
１０
n=
１０0
n＝
10０0
ｎ =
10000
EＹ1 =
２、0７8２1、６58３ １、50２9 1、4590
MSＥ１=
0、4315 0、0889 ０、0０570、0008
ｍ＝
2、0782 1、65８31、5029 1、４5９0
0、４３1５ 0、088９ ０、0057 0、00０8
%用mａtｌａb指令求积分
f=inline（'exp（ｘ、^2）'）
f=
Inline functｉｏn:
f（x）=eｘｐ（ｘ、^2）
>>S=quadl（f,0,１）
S＝
1、462７ﻬ第二类二重积分程序代码:
％%％构造示性函数,求不同区域上积分只需更改示性函数
ｆunｃtion I=I２（x,y）
ifx＾2＋ｙ＾2<=1
I=１;
ｅlｓe
I＝0;
end
%保存为I2、m
%第二类二重积分函数主体
％，程序保存为ｆ22、m
fuｎｃｔiｏｎ ｏutf2２=f２2（）
g2=input（'输入二元被积函数如1、/（1+x、^4+y、^４）、^0、5:
＇,＇s＇）%输入被积函数
ｇ2＝inｌｉｎe（g２,'x','y＇）;
fprinｔf（'输入样本容量10^V１*１0^V1--10＾V2*1０^V２:
＼r'）
V=ｚｅros（1,2）;
V
（1）=inpｕt（'V１:
'）;%输入样本容量
V
（2）＝inｐｕt（＇V２:
'）;
forｍ=Ｖ
（1）：
V（２）%样本容量10＾m１－-1０^m2
n=10＾m
ｆorj=1:
１0
x=randn（１，n）;
y=raｎdｎ（1，ｎ）;
fori＝1:
n
t２（i）＝Ｉ2（x（i）,y（i））;%示性及求与向量
end
y２=ｇ２（ｘ，y）*（2*pi）、*exp（（x、^2+ｙ、^2）／2）;
Ｙ2（j）=y2*t2'/ｎ；％单次实验样本均值
eｎd
t=oｎes（１,1０）；
EＹ=Ｙ2*t'/1０; ％十次均值
d=0;
fori=1:
10
d=d+（Y2（i）-EＹ）^2;
end
d=ｄ/（10－１）;
EＹ２（m-Ｖ
（1）+1）=ＥY; %样本容量为10^m时得样本均值
ＭSE2（m-V
（1）＋1）=d； %方差
enｄ
EＹ2，ＭＳE2
oｕtf2２=[EＹ2;MSE2］;%存放样本数字特征
%第二类二重积分,程序保存为ｆ22、m
运行结果:
%估计积分
ｍ＝ｆ22
输入二元被积函数如1、/（1+x、^４+y、^４）、^0、５：
1、／（１+x、＾4+y、^4）、^0、５
g2=
１、/（1＋ｘ、^4+y、＾4）、^0、5
输入样本容量1０^V1*10^Ｖ1－-10^V2*10＾V2:
V1:
1
V2:
4
n =
10
n=
100
ｎ=
10０0
ｎ=
10000
EY2＝
3、07592、9６９９ ２、８5６6 2、８26９
MSE２=
１、３２６7 0、０9000、00６０ 0、001４
m=
３、0７５９ 2、9６99 2、8５662、826９
１、3267 0、０90００、0０60 0、0014
实验结果整理:
第一类一重积分:
估计积分
积分真值：
1ﻩ积分估计值:
1、０0１8
样本容量:
10ﻩ １0０ﻩ1000 １0０001０0０0０
样本均值:
1、2635 １、０088 1、00６6 1、０109 1、001８
绝对误差:
0、26３50、00880、00660、01090、0０18
相对误差：
0、2６35 ０、008８ 0、0０66 0、0109 0、0018
均方误差:
0、６４39 ０、０2０５ 0、0028 0、０006 0、0001
估计积分
积分真值：
０、８8６2ﻩ积分估计值:
0、887１
样本容量:
10 1０0ﻩ 1000 １0000
样本均值：
0、９3330、9077 0、8873０、８871
绝对误差:
0、0470０、０２150、001０ 0、0０09
相对误差:
０、0５31０、0243０、0012 0、00１０
均方误差:
0、1927 0、01１2 0、0０16 0、00００
第一类二重积分:
估计积分
积分真值:
5、３9８1积分估计值:
5、４04１
样本容量:
１0 100ﻩ100０ 1０000
样本均值:
4、７７０25、1２50 ５、４３175、4041
绝对误差:
０、６27９０、2732 0、0335 0、0060
相对误差:
0、1１63 0、0506 0、00620、001１
均方误差:
３、8965 0、５5６４０、0247 ０、０017
第二类一重积分:
估计积分
积分估计值:
1、4590
样本容量:
１0 1０0ﻩ10０01０0０0
样本均值:
2、０7８21、６５83１、5029 １、4590
样本方差:
0、43１5 0、08890、00570、0０08
用ｍａtlab指令求得积分结果1、4627
第二类二重积分:
估计积分
积分估计值：
２、82６9
样本容量：
10100ﻩ１000 1000０
样本均值:
3、0759 2、９69９ 2、85662、8269
样本方差:
１、３２６7 0、0900０、０0６0 0、0014
实验结果分析:
从第一类积分瞧,以估计积分为例:
积分真值:
1积分估计值:
1、00１8
样本容量:
10ﻩ １0０ﻩ 1000 10０00100000
样本均值:
１、26３5１、00８81、0066 1、01０91、0０18
绝对误差:
0、263５ 0、０088 ０、00660、0１09 0、00１8
相对误差：
0、2６35０、００88 0、００６60、01０9 0、001８
均方误差:
0、6439 0、０2０5０、00２８ 0、0０06 0、０００1
随着样本容量得增大,样本均值有接近积分真值得趋势，绝对误差、相对误差、均方误差呈减小趋势；随着样本容量得增大,样本均值有接近积分真值得趋势，说明估计具有无偏性;绝对误差、相对误差、均方误差呈减小趋势,说明增大样本容量能提高估计精度;验证了蒙特卡洛方法估计积分值得可行性,为后续估计第二类积分提供了参考。
从第二类积分瞧,以估计积分为例：
积分估计值:
１、４59０
样本容量:
10 １0010０0 100００
样本均值:
2、078２1、6583１、5029 1、４５９０
样本方差:
0、431５ 0、０8890、0０57０、０0０８
用matlab 指令求得积分结果１、４627
由于积分真值未知，无法直接比较估计值与积分值值;但随样本容量增大，样本方差减小,间接反映了估计精度得提高。
蒙特卡洛方法估计值1、459０相比用maｔlab指令求得得积分结果1、4６27,绝对偏差０、0038,相对偏差0、０0２5。
蒙特卡洛方法估计值与用maｔlaｂ 指令求得得积分结果相互验证。
总结与讨论:
蒙特卡洛方法就是基于随机数得一种统计方法。
蒙特卡洛方法估计积分值,总得思想就是将积分改写为某个随机变量得数学期望,借助相应得随机数，利用样本均值估计数学期望，从而估计相应得积分值。
为使方法具有一般性，概率密度函数一重积分选择了,二重积分选用。
程序设计方面,本着使程序具有一般性以及方便以后使用得原则,依据问题分四类:
第一类一重积分;第一类二重积分;第二类一重积分,第二类二重积分，相应程序设计成四类,并存储为、m文件,用蒙特卡洛方法估计积分值，一重积分只需调用相应程序即可;二重积分只需依据积分域修改相应示性函数即可调用相应函数求解。
极大方便了同类问题求解。
实验运行结果表明本方案可操作性良好。
遗留问题:
本次实验未设计选用不同概率密度函数,估计精度得比较,留有不同条件下选用何种概率密度函数估计效果最佳?
如何缩短程序运行时间?
如何对程序进行封装?
如何更好评价第二类积分估计值无偏性以及精度?
等问题。
姓名:
王宏辉
班级:
材料43
学号：
２140２01060
升级会员 