2022年用蒙特卡洛方法估计积分方法及matlab编程实现 .pdf
用蒙特卡洛方法估计积分方法及 matlab 编程实现专业班级：材料 43 学生：王宏辉学号： 2140201060 指导教师：李耀武完成时间：2016 年 6 月 8 日精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 32 页用蒙特卡洛方法估计积分方法及 matlab 编程实现实验内容：1 用蒙特卡洛方法估计积分20sinxxdx，2-0 xedx和22221xyxyedxdy的值，并将估计值与真值进行比较。2 用蒙特卡洛方法估计积分210 xe dx和2244111xydxdyxy的值，并对误差进行估计。要求：1针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；2利用电脑产生所选分布的随机数以估计积分值；3进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差针对第 1 类题 或样本方差针对第 2 类题以评价估计结果的精度。目的：1能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；2 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 32 页3 能用 MATLAB 熟练进行样本的一元回归分析。实验原理：蒙特卡洛方法估计积分值， 总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。具体操作如下：一般地，积分bdxxga)(S改写成bbdxfhdxfgaa)(x)(x)(xf(x)(xS的形式， 其中为)f(x一随机变量X 的概率密度函数，且)f(x的支持域）（bf,a0)(x|x) ，f(x)(x)(xgh; 令 Y=h(X)，则积分 S=E Y ；利用matlab软 件 ， 编 程 产 生 随 机 变 量X的 随 机 数 ， 在 由)b(a,)b(a,01I(x), )(x)(xyxxIh，得到随机变量 Y的随机数，求出样本均值，以此估计积分值。积分AdxdygS)y(x,的求法与上述方法类似，在此不赘述。概率密度函数的选取：一重积分，由于要求)f(x的支持域）（bf,a0)(x|x，为使方法普遍适用，考虑到标准正态分布概率密度函数22e21)(xxf支持域为R，故选用22e21)(xxf。类似的，二重积分选用22221)y(x,yxef，支持域为2R。估计评价：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 32 页进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误针对第 1 类题，积得出或样本方差针对第2 类题，积不出以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存为一个 .m 文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存为一个.m 文本，被积函数键盘输入，示性函数用function 语句构造，求不同区域二重积分，只需改变function 函数内容。编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。程序代码及运行结果：第一类一重积分程序代码：%构造示性函数function I=I1(x,a,b) if x=a&x=b I=1; else I=0; end 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 32 页% % 第一类一重积分 , 程序主体：function outf11=f11() g1=input(输入一元被积函数如x.*sin(x):,s)%输入被积函数g1=inline(g1); a=input(输入积分下界 a:);% 输入积分上下限b=input(输入积分上界 b:); Real=input(积分真值 :);% 输入积分真值fprintf(输入样本容量 10V1-10V2:r) V=zeros(1,2); V(1)=input(V1:);%输入样本容量V(2)=input(V2:); for m=V(1):V(2)% 样本容量 10m1-10m2 n=10m for j=1:10 x=randn(1,n); for i=1:n t1(i)=I1(x(i),a,b);%示性及求和向量end 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 32 页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; for i=1:10 d=d+(Y1(i)-Real)2; end d=d/(10-1); EY1(m-V(1)+1)=EY; %样本容量为 10m时的样本均值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; % 存放样本数字特征精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 32 页运行结果：%估计积分20sinxxdx，积分真值为 1 m=f11 输入一元被积函数如x.*sin(x):x.*sin(x) g1 = x.*sin(x) 输入积分下界 a:0 输入积分上界 b:pi/2 积分真值 :1 输入样本容量10V1-10V2: V1:1 V2:5 n = 10 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 32 页n = 100 n = 1000 n = 10000 n = 100000 Real = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 32 页1 EY1 = D1 = RD1 = MSE1 = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 32 页m= %估计积分2-0 xedxM=f11 输入一元被积函数如x.*sin(x):exp(-x.2) g1 = exp(-x.2) 输入积分下界 a:0 输入积分上界 b:+inf 输入样本容量10V1-10V2: V1:1 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 32 页V2:4 n = 10 n = 100 n = 1000 n = 10000 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 32 页Real = EY1 = D1 = RD1 = MSE1 = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 32 页M = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 32 页第一类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y) if x2+y2=a&x S=quadl(f,0,1) S = 1.4627精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 23 页，共 32 页第二类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y) if x2+y2=1 I=1; else I=0; end % 第二类二重积分函数主体function outf22=f22() g2=input(输入二元被积函数如1./(1+x.4+y.4).0.5:,s)%输入被积函数g2=inline(g2,x,y); fprintf(输入样本容量 10V1*10V1-10V2*10V2:r) V=zeros(1,2); V(1)=input(V1:);%输入样本容量V(2)=input(V2:); for m=V(1):V(2)%样本容量 10m1-10m2 n=10m 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 24 页，共 32 页for j=1:10 x=randn(1,n); y=randn(1,n); for i=1:n 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; %单次实验样本均值end t=ones(1,10); EY=Y2*t/10; %十次均值d=0; for i=1:10 d=d+(Y2(i)-EY)2; end d=d/(10-1); EY2(m-V(1)+1)=EY; %样本容量为 10m时的样本均值MSE2(m-V(1)+1)=d; %方差end EY2,MSE2 outf22=EY2;MSE2; % 存放样本数字特征精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 25 页，共 32 页运行结果：%估计积分2244111xydxdyxym=f22 g2 = 输入样本容量10V1*10V1-10V2*10V2: V1:1 V2:4 n = 10 n = 100 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 26 页，共 32 页n = 1000 n = 10000 EY2 = MSE2 = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 27 页，共 32 页m = 实验结果整理：第一类一重积分：估计积分20sinxxdx积分真值： 1 积分估计值：样本容量： 10 100 1000 10000 100000 样本均值：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 28 页，共 32 页绝对误差：相对误差：均方误差：估计积分2-0 xedx积分真值： 0.8862 积分估计样本容量： 10 100 1000 10000 第一类二重积分：估计积分22221xyxyedxdy积分真值： 5.3981 积分估计值：样本容量： 10 100 1000 10000 样本均值：绝对误差：相对误差：均方误差：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 29 页，共 32 页第二类一重积分：估计积分210 xe dx积分估计样本容量： 10 100 1000 10000 第二类二重积分：估计积分2244111xydxdyxy积分估计样本容量： 10 100 1000 10000 实验结果分析：从第一类积分看，以估计积分20sinxxdx为例：积分真值： 1 积分估计值：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 30 页，共 32 页样本容量： 10 100 1000 10000 100000 样本均值：绝对误差：相对误差：均方误差：随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误差、相对误差、均方误差呈减小趋势；随着样本容量的增大，样本均值有接近积分真值的趋势，说明估计具有无偏性；绝对误差、相对误差、均方误差呈减小趋势， 说明增大样本容量能提高估计精度；验证了蒙特卡洛方法估计积分值的可行性，为后续估计第二类积分提供了参考。从第二类积分看，以估计积分210 xe dx为例：样本容量： 10 100 1000 10000 由于积分真值未知，无法直接比较估计值与积分值值；但随样本容量增大，样本方差减小， 间接反映了估计精度的提高。蒙特卡洛方法估计值 1.4590相比用 matlab 指令求得的积分结果1.4627，绝对偏差 0.0038，相对偏差 0.0025。蒙特卡洛方法估计值与用matlab 指令求得的积分结果相互验证。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 31 页，共 32 页总结与讨论：蒙特卡洛方法是基于随机数的一种统计方法。蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。为 使 方 法 具 有 一 般 性 ， 概 率 密 度 函 数 一 重 积 分 选 择 了22e21)(xxf，二重积分选用22221)y(x,yxef。程序设计方面，本着使程序具有一般性以及方便以后使用的原则，依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类，并存储为.m 文件，用蒙特卡洛方法估计积分值， 一重积分只需调用相应程序即可；二重积分只需依据积分域修改相应示性函数即可调用相应函数求解。极大方便了同类问题求解。实验运行结果说明本方案可操作性良好。遗留问题：本次实验未设计选用不同概率密度函数，估计精度的比较，留有不同条件下选用何种概率密度函数估计效果最正确？如何缩短程序运行时间？如何对程序进行封装？如何更好评价第二类积分估计值无偏性以及精度？等问题。：王宏辉班级：材料 43 学号： 2140201060精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 32 页，共 32 页
