用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx
- 文档编号:13553523
- 上传时间:2022-10-11
- 格式:DOCX
- 页数:19
- 大小:43.89KB
用蒙特卡洛方法估计积分方法及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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 用蒙特卡洛 方法 估计 积分 matlab 编程 实现