用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx
《用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx》由会员分享，可在线阅读，更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx（19页珍藏版）》请在冰豆网上搜索。
（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）;
%示性及求和向量
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;
d=d+（Y1（i）-Real）A2;
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;
%方差
Real,EY1,D1,RD1,MSE1
outf11=[EY1;
D1;
RD1;
MSE1];
%存放样本数字特征
运行结果：
2
%古计积分xsinxdx，积分真值为1
m=f11
输入一元被积函数如x.*sin（x）:
x.*sin（x）
g1=x.*sin（x）
pi/2
输入样本容量10AV1--10AV2:
5
n二
100
1000
10000
100000
Real=
EY1=
1.26351.00881.00661.01091.0018
D1=
0.2635
0.0088
0.0066
0.0109
0.0018
RD1=
MSE1=
0.6439
0.0205
0.0028
0.0006
0.0001
m=
1.2635
1.0088
1.0066
1.0109
1.0018
%估计积分e-x2dx真值为0.8862
M=f11
输入一元被积函数如x.*sin（x）:
exp（-x42）
g1=
exp（-x.A2）
+inf
积分真值：
piA0.5/2%0.8862
4
0.8862
0.9333
0.9077
0.8873
0.8871
0.0470
0.0215
0.0010
0.0009
0.0531
0.0243
0.0012
0.1927
0.0112
0.0016
0.0000
M=
0.19270.01120.0016
第一类二重积分程序代码:
%%%造示性函数，求不同区域上积分只需更改示性函数
functionI=l2（x,y）
ifxA2+yA2<
=1
%保存为I2.m
%第一类二重积分程序主体
%保存为f12.m
functionoutf12二f12（）
g2=input（'
输入二元被积函数如exp（x.A2+y.A2）:
）%
入被积函数
g2=inline（g2,'
x'
y'
%输入积分真值
fprintf（'
输入样本容量10AV1*10AV1--10AV2*10AV2:
）
V=zeros（1,2）;
y=randn（1,n）;
t2（i）=l2（x（i）,y（i））;
y2=g2（x,y）*（2*pi）.*exp（（x.A2+y.A2）/2）;
Y2（j）=y2*t2'
t=ones（1,10）;
EY二Y2*t'
d=d+（Y2（i）-RealF2;
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];
%古计积分ex2y2dxdy，真值为pi*（exp
（1）-1）%5.3981
x^y2<
m=f12
exp（x.A2+y.A2）
g2=
pi*（exp
（1）-1）%5.3981
输入样本容量10AV1*10AV1--10AV2*10AV2:
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=
第二类一重积分程序代码:
I=1;
I=0;
%第二类一重积分程序主体
%程序保存为f21.m
fun
升级会员 