用蒙特卡洛方法估计积分方法及matlab编程实现.docx
《用蒙特卡洛方法估计积分方法及matlab编程实现.docx》由会员分享，可在线阅读，更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现.docx（19页珍藏版）》请在冰豆网上搜索。
用蒙特卡洛方法估计积分方法及matlab编程实现
用蒙特卡洛方法估计积分
方法及matlab编程实现
专业班级：
材料43
学生姓名：
王宏辉
学号：
指导教师：
李耀武
完成时间：
2016年6月8日
用蒙特卡洛方法估计积分
方法及matlab编程实现
实验内容:
31
2乂
1用蒙特卡洛方法估计积分xsinxdx,e-xdx和iiexydxdy的值，
00x2：
:
;y2d
并将估计值与真值进行比较。
2用蒙特卡洛方法估计积分exdx和dxdy的值，并对误
0x2#gJl+x4+y4
差进行估计。
要求：
（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方
法；
（2）利用计算机产生所选分布的随机数以估计积分值；
（3）进行重复试验，通过计算样本均值以评价估计的无偏性；
通过计算均方误差（针对第1类题）或样本方差（针对第2类题）
以评价估计结果的精度。
目的：
（1）能通过MATLAB或其他数学软件了解随机变量的概率密
度、分布函数及其期望、方差、协方差等；
（2）熟练使用MATLAB对样本进行基本统计，从而获取数据的
基本信息；
（3）能用MATLAB熟练进行样本的一元回归分析
实验原理：
蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。
具体操作如下:
式，（其中为f（x）一随机变量X的概率密度函数，且f（x）的支持域
{x|f（x）〉O}=（a,b））,h（x^-g^x））；令Y=h（X），则积分S=E（Y；利用f（x）
matlab软件，编程产生随机变量X的随机数，在由y=h（x）l（x）,l（x）=」公匸⑻①，得到随机变量丫的随机数，求出样本均
0X（a,b）
值，以此估计积分值。
积分S=..g（x,y）dxdy的求法与上述方法类似，在此不赘述。
A
概率密度函数的选取：
一重积分，由于要求f（x）的支持域｛x|f（x）・0｝二（a,b），为使方法普
1
故选用f（xF
彳¥廿
类似的，二重积分选用f（x,y）—e2，支持域为R2
2兀
估计评价:
进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度。
程序设计：
依据问题分四类：
第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。
为了使程序具有一般性以及方便以后使用：
一重积分，程序保存为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存为一个.m文本，被积函数键盘输入，示性函数用function语句构造，求不同区域二重积分，只需改变function函
数内容。
编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。
程序代码及运行结果：
第一类一重积分程序代码：
%%构造示性函数
functionI=l1（x,a,b）
ifx>=a&&x<=b
l=1;
else
l=0;
end
%保存为I1.m
%%%%%%%%%%%%%%%%
%第一类一重积分，程序主体：
%保存为f11.m
functionoutf11=f11（）
输入被积
g1=input（'输入一元被积函数如x.*sin（x）:
','s'）%函数
g1=inline（g1）;
a=input（'输入积分下界a:
'）;%输入积分上下限
b=input（'输入积分上界b:
'）;
Real二input（'积分真值:
'）;%输入积分真值fprintf（'输入样本容量10AV1--10AV2:
\r'）V=zeros（1,2）;
V
（1）=input（'V1:
'）;%输入样本容量
V
（2）=input（'V2:
'）;
form=V
（1）:
V
（2）%样本容量10八口1--10八口2
n=10八口
forj=1:
10
x=randn（1,n）;
fori=1:
n
t1（i）=l1（x（i）,a,b）;%示性及求和向量
end
y仁g1（x）*（（pi*2）A0.5）.*exp（x.A2/2）;
Y1（j）=y1*t1'/n;%单次实验样本均值
endt=ones（1,10）;
EY二Y1*t'/10;%十次均值
D=abs（EY-Real）;%绝对误差
RD二D/Real;%绝对误差
d=0;
fori=1:
10
d=d+（Y1（i）-Real）A2;
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];%存放样本数字特征
%保存为f11.m
运行结果：
2
%古计积分xsinxdx，积分真值为1
0
m=f11
输入一元被积函数如x.*sin（x）:
x.*sin（x）
g1=x.*sin（x）
输入积分下界a:
0
输入积分上界b:
pi/2
积分真值:
1
输入样本容量10AV1--10AV2:
V1:
1
V2:
5
n二
10
n二
100
n二
1000
n二
10000
n二
100000
Real=
1
EY1=
1.26351.00881.00661.01091.0018
D1=
0.2635
0.0088
0.0066
0.0109
0.0018
RD1=
0.2635
0.0088
0.0066
0.0109
0.0018
MSE1=
0.6439
0.0205
0.0028
0.0006
0.0001
m=
1.2635
1.0088
1.0066
1.0109
1.0018
0.2635
0.0088
0.0066
0.0109
0.0018
0.2635
0.0088
0.0066
0.0109
0.0018
0.6439
0.0205
0.0028
0.0006
0.0001
%估计积分e-x2dx真值为0.8862
0
M=f11
输入一元被积函数如x.*sin（x）:
exp（-x42）
g1=
exp（-x.A2）
输入积分下界a:
0
输入积分上界b:
+inf
积分真值：
piA0.5/2%0.8862
输入样本容量10AV1--10AV2:
V1:
1
V2:
4
n二
100
n二
1000
n二
10000
Real=
0.8862
EY1=
0.9333
0.9077
0.8873
0.8871
D1=
0.0470
0.0215
0.0010
0.0009
RD1=
0.0531
0.0243
0.0012
0.0010
MSE1=
0.1927
0.0112
0.0016
0.0000
M=
0.9333
0.9077
0.8873
0.8871
0.0470
0.0215
0.0010
0.0009
0.0531
0.0243
0.0012
0.0010
0.19270.01120.0016
0.0000
第一类二重积分程序代码:
%%%造示性函数，求不同区域上积分只需更改示性函数
functionI=l2（x,y）
ifxA2+yA2<=1
l=1;
else
l=0;
end
%保存为I2.m
%第一类二重积分程序主体
%保存为f12.m
functionoutf12二f12（）
g2=input（'输入二元被积函数如exp（x.A2+y.A2）:
','s'）%
入被积函数
g2=inline（g2,'x','y'）;
Real二input（'积分真值:
'）;%输入积分真值
fprintf（'输入样本容量10AV1*10AV1--10AV2*10AV2:
\r'）
V=zeros（1,2）;
V
（1）=input（'V1:
'）;%输入样本容量
V
（2）=input（'V2:
'）;
form=V
（1）:
V
（2）%样本容量10八口1--10八口2
n=10八口
forj=1:
10
x=randn（1,n）;
y=randn（1,n）;
fori=1:
n
t2（i）=l2（x（i）,y（i））;%示性及求和向量
end
y2=g2（x,y）*（2*pi）.*exp（（x.A2+y.A2）/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）-RealF2;
end
d=d/（10-1）;
EY2（m-V
（1）+1）=EY;
%
样本容量为10am时的样本均值
D2（m-V
（1）+1）=D;
%
绝对误差
RD2（m-V
（1）+1）=RD;
%
绝对误差
MSE2（m-V
（1）+1）=d;
%
方差
Real,EY2,D2,RD2,MSE2
outf12=[EY2;D2;RD2;MSE2];%存放样本数字特征
%保存为f12.m
运行结果：
%古计积分ex2y2dxdy，真值为pi*（exp
（1）-1）%5.3981
x^y2<
m=f12
输入二元被积函数如exp（x.A2+y.A2）:
exp（x.A2+y.A2）
g2=
exp（x.A2+y.A2）
积分真值:
pi*（exp
（1）-1）%5.3981
输入样本容量10AV1*10AV1--10AV2*10AV2:
V1:
1
V2:
4
n二
10
n二
100
n二
10000
Real=
5.3981
EY2=
4.7702
5.1250
5.4317
5.4041
D2=
0.6279
0.2732
0.0335
0.0060
RD2=
0.1163
0.0506
0.0062
0.0011
MSE2=
3.8965
0.5564
0.0247
0.0017
m=
4.7702
5.1250
5.4317
5.4041
0.6279
0.2732
0.0335
0.0060
0.1163
0.0506
0.0062
0.0011
3.8965
0.5564
0.0247
0.0017
第二类一重积分程序代码:
%%构造示性函数
functionI=l1（x,a,b）
ifx>=a&&x<=b
I=1;
else
I=0;
%保存为I1.m
%第二类一重积分程序主体
%程序保存为f21.m
fun
升级会员 