用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx
- 文档编号:22081480
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:20
- 大小:30.63KB
用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx
《用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx》由会员分享,可在线阅读,更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现文档格式.docx(20页珍藏版)》请在冰豆网上搜索。
ifx>
=a&
&x<=b
I=1;
else
I=0;
end
%保存为I1、m
%%%%%%%%%%%%%%%%
%%第一类一重积分,程序主体:
%保存为f11、m
function outf11=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);
%示性及求与向量
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:
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;
%方差
Real,EY1,D1,RD1,MSE1
outf11=[EY1;
D1;
RD1;
MSE1];
%存放样本数字特征
%保存为f11、m
运行结果:
%估计积分,积分真值为1
m=f11
输入一元被积函数如x、*sin(x):
x、*sin(x)
g1=
0
输入积分上界b:
pi/2
积分真值:
1
输入样本容量 10^V1--10^V2:
V1:
V2:
5
n=
10
100
1000
n =
10000
100000
Real =
1
EY1=
1、2635 1、0088 1、0066 1、01091、0018
D1=
0、2635 0、00880、00660、0109 0、0018
RD1=
0、2635 0、00880、0066 0、0109 0、0018
MSE1=
0、64390、02050、00280、0006 0、0001
m=
1、2635 1、0088 1、0066 1、0109 1、0018
0、26350、00880、00660、01090、0018
0、2635 0、0088 0、0066 0、0109 0、0018
0、6439 0、0205 0、0028 0、00060、0001ﻩ
%估计积分真值为0、8862
M=f11
输入一元被积函数如x、*sin(x):
exp(-x、^2)
exp(-x、^2)
输入积分下界a:
+inf
pi^0、5/2%0、8862
输入样本容量10^V1--10^V2:
V1:
V2:
4
10
n=
n=
1000
10000
Real =
0、8862
EY1 =
0、93330、90770、8873 0、8871
D1 =
0、04700、0215 0、00100、0009
RD1=
0、05310、0243 0、0012 0、0010
0、1927 0、01120、00160、0000
M =
0、9333 0、90770、8873 0、8871
0、0470 0、02150、0010 0、0009
0、0531 0、0243 0、0012 0、0010
0、19270、01120、0016 0、0000
第一类二重积分程序代码:
%%%构造示性函数,求不同区域上积分只需更改示性函数
functionI=I2(x,y)
ifx^2+y^2<
=1
I=1;
else
I=0;
%保存为I2、m
%第一类二重积分程序主体
%保存为f12、m
functionoutf12=f12()
g2=input('
输入二元被积函数如exp(x、^2+y、^2):
','
)%输入被积函数
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('
for m=V
(1):
V(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
x=randn(1,n);
y=randn(1,n);
fori=1:
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;
%单次实验样本均值
t=ones(1,10);
EY=Y2*t'
/10;
%十次均值
D=abs(EY-Real);
RD=D/Real;
%绝对误差
for i=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;
Real,EY2,D2,RD2,MSE2
outf12=[EY2;
D2;
RD2;
MSE2];
%保存为f12、m
运行结果:
%估计积分,真值为pi*(exp
(1)-1)%5、3981
m=f12
输入二元被积函数如exp(x、^2+y、^2):
exp(x、^2+y、^2)
g2=
exp(x、^2+y、^2)
pi*(exp(1)-1)%5、3981
输入样本容量10^V1*10^V1--10^V2*10^V2:
V1:
V2:
10
100
1000
10000
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、00620、0011
MSE2 =
3、89650、55640、02470、0017
m=
4、7702 5、12505、43175、4041
0、6279 0、2732 0、0335 0、0060
0、11630、05060、00620、0011
3、8965 0、5564 0、02470、0017
第二类一重积分程序代码:
%%%构造示性函数
functionI=I1(x,a,b)
if x>
=a&
&
x<
=b
I=1;
I=0;
%保存为I1、m
%第二类一重积分程序主体
%程序保存为f21、m
functionoutf21=f21()
g1=input('输入一元被积函数如exp(x、^2):
','
g1=inline(g1);
a=input('输入积分下界a:
%输入积分上下限
b=input('
输入积分上界b:
fprintf('输入样本容量10^V1--10^V2:
\r'
V=zeros(1,2);
V
(1)=input('
V
(2)=input('
form=V
(1):
V
(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
x=randn(1,n);
for i=1:
t1(i)=I1(x(i),a,b);
y1=g1(x)*((pi*2)^0、5)、*exp(x、^2/2);
Y1(j)=y1*t1'
t=ones(1,10);
EY=Y1*t'/10;
%十次均值
d=0;
fori=1:
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;
%方差
EY1,MSE1
outf21=[EY1;
MSE1];
%%%%程序保存为f21、m
%估计积分
m=f21
输入一元被积函数如exp(x、^2):
exp(x、^2)
g1 =
exp(x、^2)
输入积分下界a:
1
输入样本容量10^V1--10^V2:
4
n=
10
100
1000
10000
EY1 =
2、07821、6583 1、5029 1、4590
MSE1=
0、4315 0、0889 0、00570、0008
m=
2、0782 1、65831、5029 1、4590
0、4315 0、0889 0、0057 0、0008
%用matlab指令求积分
f=inline('
exp(x、^2)'
f=
Inline function:
f(x)=exp(x、^2)
>
S=quadl(f,0,1)
S=
1、4627ﻬ第二类二重积分程序代码:
%%%构造示性函数,求不同区域上积分只需更改示性函数
function I=I2(x,y)
ifx^2+y^2<
=1
I=1;
%保存为I2、m
%第二类二重积分函数主体
%,程序保存为f22、m
function outf22=f22()
g2=input('
输入二元被积函数如1、/(1+x、^4+y、^4)、^0、5:
','s')%输入被积函数
g2=inline(g2,'
x'
'
fprintf('
输入样本容量10^V1*10^V1--10^V2*10^V2:
\r'
V=zeros(1,2);
V
(1)=input('
V
(2)=input('V2:
form=V
(1):
V(2)%样本容量10^m1--10^m2
n=10^m
forj=1:
x=randn(1,n);
y=randn(1,n);
t2(i)=I2(x(i),y(i));
y2=g2(x,y)*(2*pi)、*exp((x、^2+y、^2)/2);
Y2(j)=y2*t2'
/n;
%单次实验样本均值
t=ones(1,10);
EY=Y2*t'
/10;
%十次均值
d=0;
fori=1:
d=d+(Y2(i)-EY)^2;
d=d/(10-1);
EY2(m-V
(1)+1)=EY;
%样本容量为10^m时得样本均值
MSE2(m-V
(1)+1)=d;
%方差
EY2,MSE2
outf22=[EY2;
MSE2];
%第二类二重积分,程序保存为f22、m
%估计积分
m=f22
输入二元被积函数如1、/(1+x、^4+y、^4)、^0、5:
1、/(1+x、^4+y、^4)、^0、5
g2=
1、/(1+x、^4+y、^4)、^0、5
输入样本容量10^V1*10^V1--10^V2*10^V2:
n =
10
100
1000
10000
EY2=
3、07592、9699 2、8566 2、8269
MSE2=
1、3267 0、09000、0060 0、0014
3、0759 2、9699 2、85662、8269
1、3267 0、09000、0060 0、0014
实验结果整理:
第一类一重积分:
估计积分
1ﻩ积分估计值:
1、0018
样本容量:
10ﻩ 100ﻩ1000 10000100000
样本均值:
1、2635 1、0088 1、0066 1、0109 1、0018
绝对误差:
0、26350、00880、00660、01090、0018
相对误差:
0、2635 0、0088 0、0066 0、0109 0、0018
均方误差:
0、6439 0、0205 0、0028 0、0006 0、0001
估计积分
0、8862ﻩ积分估计值:
0、8871
10 100ﻩ 1000 10000
样本均值:
0、93330、9077 0、88730、8871
0、04700、02150、0010 0、0009
相对误差:
0、05310、02430、0012 0、0010
0、1927 0、0112 0、0016 0、0000
第一类二重积分:
5、3981积分估计值:
5、4041
10 100ﻩ1000 10000
4、77025、1250 5、43175、4041
0、62790、2732 0、0335 0、0060
0、1163 0、0506 0、00620、0011
3、8965 0、55640、0247 0、0017
第二类一重积分:
积分估计值:
1、4590
10 100ﻩ100010000
2、07821、65831、5029 1、4590
样本方差:
0、4315 0、08890、00570、0008
用matlab指令求得积分结果1、4627
第二类二重积分:
积分估计值:
2、8269
样本容量:
10100ﻩ1000 10000
3、0759 2、9699 2、85662、8269
1、3267 0、09000、0060 0、0014
实验结果分析:
从第一类积分瞧,以估计积分为例:
1积分估计值:
1、0018
10ﻩ 100ﻩ 1000 10000100000
1、26351、00881、0066 1、01091、0018
0、2635 0、0088 0、00660、0109 0、0018
0、26350、0088 0、00660、0109 0、0018
0、6439 0、02050、0028 0、0006 0、0001
随着样本容量得增大,样本均值有接近积分真值得趋势,绝对误差、相对误差、均方误差呈减小趋势;
随着样本容量得增大,样本均值有接近积分真值得趋势,说明估计具有无偏性;
绝对误差、相对误差、均方误差呈减小趋势,说明增大样本容量能提高估计精度;
验证了蒙特卡洛方法估计积分值得可行性,为后续估计第二类积分提供了参考。
从第二类积分瞧,以估计积分为例:
1、4590
10 1001000 10000
2、07821、65831、5029 1、4590
0、4315 0、08890、00570、0008
用matlab 指令求得积分结果1、4627
由于积分真值未知,无法直接比较估计值与积分值值;
但随样本容量增大,样本方差减小,间接反映了估计精度得提高。
蒙特卡洛方法估计值1、4590相比用matlab指令求得得积分结果1、4627,绝对偏差0、0038,相对偏差0、0025。
蒙特卡洛方法估计值与用matlab 指令求得得积分结果相互验证。
总结与讨论:
蒙特卡洛方法就是基于随机数得一种统计方法。
蒙特卡洛方法估计积分值,总得思想就是将积分改写为某个随机变量得数学期望,借助相应得随机数,利用样本均值估计数学期望,从而估计相应得积分值。
为使方法具有一般性,概率密度函数一重积分选择了,二重积分选用。
程序设计方面,本着使程序具有一般性以及方便以后使用得原则,依据问题分四类:
第二类一重积分,第二类二重积分,相应程序设计成四类,并存储为、m文件,用蒙特卡洛方法估计积分值,一重积分只需调用相应程序即可;
二重积分只需依据积分域修改相应示性函数即可调用相应函数求解。
极大方便了同类问题求解。
实验运行结果表明本方案可操作性良好。
遗留问题:
本次实验未设计选用不同概率密度函数,估计精度得比较,留有不同条件下选用何种概率密度函数估计效果最佳?
如何缩短程序运行时间?
如何对程序进行封装?
如何更好评价第二类积分估计值无偏性以及精度?
等问题。
姓名:
王宏辉
班级:
材料43
学号:
2140201060
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 用蒙特卡洛 方法 估计 积分 matlab 编程 实现