用蒙特卡洛方法估计积分方法及matlab编程实现

1、 用蒙特卡洛方法估计积分 方法及matlab编程实现 专业班级：材料43 学生姓名：王宏辉 学 号：60 指导教师：李耀武 完成时间：2016年6月8日 用蒙特卡洛方法估计积分 方法及matlab编程实现 实验内容: 1用蒙特卡洛方法估计积分 xsin xdx , o e-x dx 和ex y dxdy 的值， 0 x2 y2 1 并将估计值与真值进行比较。 2用蒙特卡洛方法估计积分edx和1 dxdy的值，并对误 0 x2y2ijl x4 y4 差进行估计 要求： （1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方 法； （2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复

2、试验，通过计算样本均值以评价估计的无偏性； 通过计算均方误差（针对第1类题）或样本方差（针对第2类题） 以评价估计结果的精度。 目的: （1）能通过MATLAB或其他数学软件了解随机变量的概率密 度、分布函数及其期望、方差、协方差等； （2）熟练使用MATLAB对样本进行基本统计，从而获取数据的 基本信息; (3)能用MATLAB熟练进行样本的一元回归分析 实验原理： 蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变 量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从 而估计相应的积分值。 具体操作如下: 般地，积分 b g(x)dx改写成S bg(x) a f(X) f (

3、x)dx b h(x) f (x)dx 的形 2 式，(其中为f(x) 一随机变量 X的概率密度函数，且f(x)的支持域 x|f(x) 0 (a,b) , h(x) 型)；令 Y=h(X)，则积分 S=E( Y)；利用 f(x) matlab软件，编程产生随机变量X的随机数，在由 y h(x)I(x), I(x) 1 ,x (a,b)，得到随机变量丫的随机数，求出样本均 0 ,x (a,b) 值，以此估计积分值。 积分S g(x, y)dxdy的求法与上述方法类似，在此不赘述。 A 概率密度函数的选取： (a,b)，为使方法普 一重积分，由于要求f(x)的支持域x| f(x) 0 遍适用，考虑

4、到标准正态分布概率密度函数 f(x) 1 x T支持域为 R ,故选用f(x) 类似的，二重积分选用f (x, y) 1e 支持域为R2 估计评价: 进行重复试验，通过计算样本均值以评价估计的无偏性；通过计 算均方误（针对第 1类题，积得出）或样本方差（针对第 2 类题，积 不出）以评价估计结果的精度。 程序设计： 依据问题分四类：第一类一重积分；第一类二重积分；第二类一 重积分，第二类二重积分，相应程序设计成四类。 为了使程序具有一般性以及方便以后使用：一重积分，程序保存 为一个 .m 文本，被积函数，积分区间均采用键盘输入；二重积分， 程序主体保存为一个 .m 文本，被积函数键盘输入，示性

5、函数用 function 语句构造，求不同区域二重积分，只需改变 function 函 数内容。 编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。 程序代码及运行结果： 第一类一重积分程序代码： %构%造示性函数 function I=I1（x,a,b） if x=aD1;RD1;MSE1; % 存放样本数字特征 %保存为 运行结果: %古计积分xsin xdx，积分真值为1 0 m=f11 输入一元被积函数如x.*sin(x):x.*sin(x) g1 = x.*si n(x) 输入积分下界a:0 输入积分上界b:pi/2 积分真值:1 输入样本容量10AV1-10AV2: V1:1 V

6、2:5 n 二 10 100 1000 10000 100000 Real = EY1 = D1 = RD1 = MSE1 = m= 2 %估计积分e-x2dx 真值为 0 M=f11 输入一元被积函数如x.*si n(x):exp(-x42) g1 = exp(-x.A2) 输入积分下界 a:0 输入积分上界 b:+inf 积分真值 :piA2% 输入样本容量 10AV1-10AV2: V1:1 V2:4 10 100 1000 10000 Real = EY1 = D1 = RD1 = MSE1 = 第一类二重积分程序代码： %构%造示性函数，求不同区域上积分只需更改示性函数 functi

7、on I=I2(x,y) if xA2+yA2 A . 方差 end Real,EY2,D2,RD2,MSE2 outf12=EY2;D2;RD2;MSE2; % 存放样本数字特征 %保存为 运行结果： %估计积分ex2 y2 dxdy ，真值为 pi*(exp(1)-1)% x2 y2 1 m=f12 输入二元被积函数如 exp(x.A2+y.A2):exp(x.A2+y.A2) g2 = exp(x.A2+y.A2) 积分真值 :pi*(exp(1)-1)% 输入样本容量 10AV1*10AV1-10AV2*10AV2: V1:1 V2:4 n = 10 100 1000 10000 Re

8、al = EY2 = D2 = RD2 = MSE2 = 第二类一重积分程序代码： %构%造示性函数 function I=I1(x,a,b) if x=a else I=0; end %保存为 %第二类二重积分函数主体 %，程序保存为 function outf22=f22() 输入 g2=input( 输入二元被积函数如 1./(1+x.A4+y.A4).A:,s)% 被积函数 g2=inline(g2,x,y); fprintf( 输入样本容量 10AV1*10AV1-10AV2*10AV2:r) V=zeros(1,2); V(1)=input(V1:);% 输入样本容量 V(2)=i

9、nput(V2:); for m=V(1):V(2)% 样本容量 10Am1-10Am2 n=10Am 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.A2+y.A2)/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)-E 丫)八2; end d=d/(10-1); EY 2(m-V(1)+1

10、)=E Y;%样本容量为10Am时的样本均值 MSE2(m-V(1)+1)=d;% 方差 end EY2,MSE2 outf22=EY2;MSE2; % 存放样本数字特征 %第二类二重积分，程序保存为 运行结果: %古计积分 dxdy x2 y2 J x4 y4 m=f22 输入二元被积函数如 1./(1+x.A4+y.A4).A:1./(1+x.A4+y.A4).A g2 = 1./(1+x.A4+y.A4).A 输入样本容量 10AV1*10AV1-10AV2*10AV2: V1:1 V2:4 n 二 10 100 1000 10000 EY2 = MSE2 = 实验结果整理： 第一类一重

11、积分： 2 估计积分 xsin xdx o 积分真值：1积分估计值： 100000 样本容量：10100100010000 样本均值: 绝对误差： 相对误差： 均方误差： 估计积分 e- x dx 0 积分真值： 样本容量： 10 样本均值： 绝对误差： 相对误差： 均方误差： 积分估计值： 100 1000 10000 第一类二重积分： 估计积分ex y dxdy x2 y2 1 积分真值： 积分估计值： 样本容量： 10 100 1000 样本均值： 绝对误差： 相对误差： 10000 均方误差： 10000 10000 xdx为例: 第二类一重积分: 1 估计积分 ex dx 0 积分估

12、计值： 样本容量：101001000 样本均值： 样本方差： 用matlab指令求得积分结果 第二类二重积分： 估计积分f 1 dxdy x2 y2 1J1 x4 y4 积分估计值： 样本容量：101001000 样本均值： 样本方差： 实验结果分析： 2 从第一类积分看，以估计积分xsin 0 积分真值：1积分估计值: 样本容量： 10 100 1000 10000 100000 样本均值： 绝对误差： 相对误差： 均方误差： 随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误 差、相对误差、均方误差呈减小趋势；随着样本容量的增大，样本均 值有接近积分真值的趋势，说明估计具有无偏性；绝

13、对误差、相对误 差、均方误差呈减小趋势， 说明增大样本容量能提高估计精度；验证 了蒙特卡洛方法估计积分值的可行性， 为后续估计第二类积分提供了 1 从第二类积分看，以估计积分 ex dx 为例： 0 积分估计值： 样本容量： 10 100 1000 10000 样本均值： 样本方差： 用 matlab 指令求得积分结果 由于积分真值未知，无法直接比较估计值与积分值值；但随样本 容量增大，样本方差减小， 间接反映了估计精度的提高。蒙特卡洛方 法估计值相比用 matlab 指令求得的积分结果， 绝对偏差，相对偏差 蒙特卡洛方法估计值与用 matlab 指令求得的积分结果相互验证。 总结与讨论: 蒙特卡洛方法是基于随机数的一种统计方法。 蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变 量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从 而估计相应的积分值。 为使方法具有一般性，概率密度函数一重积分选择了 1x21 x2 y2 f(x) =e 2，二重积分选用 f(x, y) e 2 。 J22 程序设计方面，本着使程序具有一般性以及方便以后使用的原则， 依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积 分，第二类二重积分，相应程序设计成四类，并存储为.m