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

1、用蒙特卡洛方法估计积分方法及matlab编程实现专业班级：材料43学生姓名：王宏辉学 号：2140201060指导教师：李耀武 完成时间：2016年6月8日用蒙特卡洛方法估计积分方法及matlab 编程实现实验内容：1用蒙特卡洛方法估计积分Jxsinxdx, Je-xdx和JJ e*“dxdy的00x2：；y2<1值，并将估计值与真值进行比较。 12用家特卡洛方法估计积分(e dx和1 , 4 4 dxdy的值,0x2 72 j . 1 x y并对误差进行估计。要求：(1)针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；(2 )利用计算机产生所选分布的随机数以估计积分值；(3)进行重

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

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

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

5、同区域二重积分，只需改变function函数内容编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。程序代码及运行结果：第一类一重积分程序代码：%构造示性函数function I=I1(x,a,b)if x>=a&&x<二bI=1;elseI=0;end%保存为I1.m%第一类一重积分，程序主体：%保存为fll.mfunction outf11=f11()g1=input('输入一元被积函数如 x.*sin(x):','s')%输入被积函数g1=inline(g1);a=input('输入积分下界a:');%输入积分上

6、下限b=input('输入积分上界 b:');Real=input('积分真值：');输入积分真值fprintf('输入样本容量10AV1-10AV2:r')V=zeros(1,2);V(1)=input('V1:');%输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量end y1=g1(x)*(pi*2)A0.5).*ex

7、p(x.A2/2);Y1(j)=y1*t1'/n; % 单次实验样本均值 end t=ones(1,10);EY=Y1*t'/10;%十次均值D=abs(EY-Real); % 绝对误差 RD=D/Real;%绝对误差d=0;foh=1:10d=d+(Y1(i)-Real)A2;%样本容量为 10Am%绝对误差%绝对误差时的end d=d/(10-1); EY1(m-V(1)+1)=EY;样本均值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;

8、RD1;MSE1; % 存放样本数字特征 %保存为fll.m运行结果:%估计积分fxsinxdx,积分真值为1m=f11输入一元被积函数如x.*sin(x):x.*sin(x) g1 = x.*sin(x)输入积分下界a:0输入积分上界b:pi/2积分真值:1输入样本容量10AV1-10AV2:V1:1V2:510100 n =1000n =10000n =100000Real =11.0018EY1 =1.26351.00881.00661.0109D1 =0.26350.00880.00660.01090.00180.00180.00011.00180.00180.00180.00010.

9、64390.02050.00280.0006m=1.26351.00881.00661.01090.26350.00880.00660.01090.26350.00880.00660.01090.64390.02050.00280.0006RD1 =0.26350.00880.00660.0109MSE1 =%估计积分 ；e-x2dx真值为0.8862 0M=f11输入一元被积函数如x.*sin(x):exp(-x.A2)g1 =exp(-xd2)输入积分下界a:0输入积分上界b:+inf积分真值：piA0.5/2%0.8862输入样本容量 10AV1-10AV2:V1:1V2:4101001

10、00010000Real =0.8862EY1 =0.9333D1 =0.0470RD1 =0.0531MSE1 =0.1927M =0.93330.04700.05310.19270.90770.88730.88710.02150.00100.00090.02430.00120.00100.01120.00160.00000.90770.88730.88710.02150.00100.00090.02430.00120.00100.01120.00160.0000第一类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if xA2+yA2&

11、lt;=1I=1;elseI=0;end%保存为I2.m%第一类二重积分程序主体%保存为f12.mfunction outf12=f12()g2=input(' 输入二元被积函数如exp(x.A2+yA2)：','s')%输入被积函数g2=inline(g2,'x','y');Real=input('积分真值：');输入积分真值fprintf(' 输 入 样 本 容 量10AVi*i0Avi-i0AV2*i0AV2：r')V=zeros(1,2);V(1)=input('V1:');

12、%输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);y=randn(1,n);for i=1:nt2(i)=I2(x(i),y(i);%示性及求和向量endy2=g2(x,y)*(2*pi).*exp(xA2+yA2)/2);Y2(j)=y2*t2'/n; %单次实验样本均值 endt=ones(1,10);EY=Y2*t'/10;D=abs(EY-Real); % RD=D/Real;d=0;foh=1:10d=d+(Y2(i)-Real)A

13、2; end d=d/(10-1);EY2(m-V(1)+1)=EY;样本均值D2(m-V(1)+1)=D;RD2(m-V(1)+1)=RD;MSE2(m-V(1)+1)=d;%十次均值绝对误差%绝对误差%样本容量为 10Am 时的%绝对误差%绝对误差%方差end存放样本数字特征Real,EY2,D2,RD2,MSE2 outf12=EY2;D2;RD2;MSE2; % %保存为f12.m运行结果:%估计积分 ex*dxdy,真值为 pi*(exp(1)-1)%5.3981 22x y也m=f12输入二元被积函数如 exp(x.A2+y.A2):exp(x.A2+y.A2)g2 =exp(x.

14、A2+yd2)积分真值：pi*(exp(1)-1)%5.3981输入样本容量 10AV1*10AV1-10AV2*10AV2：V1:1V2:4 n = 101001000n =10000Real =5.3981EY2 =4.77025.12505.43175.4041D2 =0.62790.27320.03350.0060RD2 =0.11630.05060.00620.0011MSE2 =3.89650.55640.02470.00174.77025.12505.43175.40410.62790.27320.03350.00600.11630.05060.00620.00113.89650

15、.55640.02470.0017第二类一重积分程序代码：%构造示性函数function I=I1(x,a,b)if x>=a&&x<二bI=1;elseI=0;end%保存为I1.m%第二类一重积分程序主体%程序保存为f21.mfunction outf21=f21()输入被g1=input('输入一元被积函数如exp(x.A2):','s')%积函数g1=inline(g1);a=input('输入积分下界a:');%输入积分上下限 b=input('输入积分上界 b:');fprintf('

16、;输入样本容量10AV1-10AV2:r')V=zeros(1,2);V(1)=input('V1:');% 输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量 10Am1-10Am2n=10Amfor j=1:10 x=randn(1,n); for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量endy1=g1(x)*(pi*2)A0.5).*exp(x.A2/2);Y1(j)=y1*t1'/n; % 单次实验样本均值 end t=ones(1,10);EY=Y1*t/10；%十次均值d

17、=0; fo门=1:10d=d+(Y1(i)-EY)A2;end d=d/(10-1);EY1(m-V(1)+1)=EY;%样本容量为10Am 时的样本均值MSE1(m-V(1)+1)=d;% 方差end EY1,MSE1outf21=EY1;MSE1; % 存放样本数字特征 %程序保存为f21.m运行结果：1 %估计积分exdx0m=f21输入一元被积函数如 exp(x.A2):exp(x.A2)g1 =exp(xd2)输入积分下界a:0 输入积分上界b:1 输入样本容量 10AV1-10AV2:V1:1V2:4 n =10100n = 1000n = 10000EY1 =1.45900.0

18、0082.07821.65831.5029MSE1 =0.43150.08890.0057m =2.07821.65831.50291.45900.43150.08890.00570.000812%用matlab指令求积分fe dx0f=inline('exp(x.A2)')f =Inline function:f(x) = exp(x-2)>> S=quadl(f,0,1)S =1.4627第二类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存

19、为I2.m%第二类二重积分函数主体%,程序保存为f22.mfunction outf22=f22()g2=input(' 输 入 二 元 被 积 函 数 如1./(1+xA4+yA4)A0.5：','s')%输入被积函数g2=inline(g2,'x','y');fprintf('输入样本容量10AVl*l0AVi-i0AV2*i0AV2：r')V=zeros(1,2);V(1)=input('V1:');%输入样本容量V(2)=input('V2:');for m=V(1):V(2

20、)% 样本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);y=randn(1,n);for i=1:nt2(i)=I2(x(i),y(i);% 示性及求和向量 end y2=g2(x,y)*(2*pi).*exp(xA2+yA2)/2);Y2(j)=y2*t2'/n; % 单次实验样本均值 end t=ones(1,10);EY=Y2*t'/10;%十次均值d=0; fo门=1:10d=d+(Y2(i)-EY)A2;end d=d/(10-1);EY2(m-V(1)+1)=EY;%样本容量为 10Am 时的样本均值 MSE2(m-V(1)+

21、1)=d;% 方差end EY2,MSE2 outf22=EY2;MSE2; % 存放样本数字特征 %第二类二重积分，程序保存为f22.m运行结果：%估计积分一1-dxdyx>yZ:1 x4 -y4m=f22输 入 二 元 被 积 函 数 如1./(1+x.A4+yA4)A0.5：1./(1+xA4+yA4)A0.5g2 =1./(1+xA4+yA4)A0.5输入样本容量 10AV1*10AV1-10AV2*10AV2:V1:1V2:4 n =10n =100 n =100010000EY2 =3.07592.96992.85662.8269MSE2 =1.32670.09000.006

22、00.0014m =3.07592.96992.85661.32670.09000.00602.82690.0014实验结果整理:第一类一重积分:估计积分积分真值： 样本容量： 样本均值 1.0018 绝对误差 0.0018 相对误差 0.0018均方误差 0.0001n2 xsin xdx01积分估计值:10100:1.26350.26350.26350.64391000100001000001.00881.00661.01090.00880.00660.01090.00880.00660.01090.02050.00280.00061.0018-be-x2 .e dx样本容量：101001

23、00010000样本均值：0.93330.90770.88730.8871绝对误差：0.04700.02150.00100.0009相对误差：0.05310.02430.00120.0010均方误差：0.19270.01120.00160.0000估计积分积分真值积分估计值:0.88620.88710第一类二重积分:22估计积分ex ydxdy22 dX招色积分真值：5.3981样本容量：10100积分估计值： 5.4041100010000样本均值:绝对误差:相对误差:均方误差:4.77020.62790.11633.89655.12500.27320.05060.55645.43170.0

24、3350.00620.02475.40410.00600.00110.0017第二类一重积分：12估计积分exdx0积分估计值：1.4590样本容量：10100100010000样本均值：2.07821.65831.50291.4590样本方差：0.43150.08890.00570.0008用matlab指令求得积分结果1.4627第二类二重积分：估计积分14 dxdyx2y2j-.1 X4 y4积分估计值：2.8269样本容量：101001000100002.85662.82690.00600.00141000100001.00881.00660.00880.00660.00880.00660.02050.00281000001.01090.01090.01090.0006样本均值：3.07592.9699样本方差：1.32670.0900实验结果分析:工2从第一类积分看，以估计积分 fxsinxdx为例:0积分真值：1积分估计值：1.0018样本容量：10100样本均值：1.26351.0018绝对误差：0.26350.0018相对误差：0.26350.0018均方误差：0.64390.0001随着样本容量的增大，样本均值有接近积分真值的趋势, 绝对误差、相对误差、均方误差呈减小趋势；随着样本容量 的增大，样本均值有接近积分真值的趋势，说明估计具有