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

1、 用蒙特卡洛方法估计积分方法及matlab编程实现 专业班级： 材料43学生姓名： 王宏辉学 号： 2140201060指导教师：李耀武完成时间： 2016年6月8日 用蒙特卡洛方法估计积分方法及matlab编程实现实验内容：1用蒙特卡洛方法估计积分 ，和的值，并将估计值与真值进行比较。2用蒙特卡洛方法估计积分 和的值，并对误差进行估计。要求： （1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；（2）利用计算机产生所选分布的随机数以估计积分值；（3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针 对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。目的：

2、 （1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数 及其期望、方差、协方差等；（2） 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息；（3） 能用 MATLAB 熟练进行样本的一元回归分析。实验原理： 蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。具体操作如下：一般地，积分改写成的形式，（其中为一随机变量X的概率密度函数，且的支持域)，）;令Y=h(X)，则积分S=E（Y）；利用matlab软件，编程产生随机变量X的随机数，在由，得到随机变量Y的随机数，求出样

3、本均值，以此估计积分值。积分的求法与上述方法类似，在此不赘述。概率密度函数的选取：一重积分，由于要求的支持域，为使方法普遍适用，考虑到标准正态分布概率密度函数支持域为，故选用。类似的，二重积分选用，支持域为。估计评价：进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存

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

5、ne(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=10mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b)

6、;%示性及求和向量endy1=g1(x)*(pi*2)0.5).*exp(x.2/2);Y1(j)=y1*t1'/n; %单次实验样本均值endt=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;endd=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; %方差endReal

7、,EY1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; %存放样本数字特征%保存为f11.m运行结果：%估计积分 ，积分真值为1m=f11输入一元被积函数如x.*sin(x):x.*sin(x)g1 =x.*sin(x)输入积分下界a:0输入积分上界b:pi/2积分真值:1输入样本容量 10V1-10V2:V1:1V2:5n = 10n = 100n = 1000n = 10000n = 100000Real = 1EY1 = 1.2635 1.0088 1.0066 1.0109 1.0018D1 = 0.2635 0.0088 0.0066 0.0109 0.001

8、8RD1 = 0.2635 0.0088 0.0066 0.0109 0.0018MSE1 =0.6439 0.0205 0.0028 0.0006 0.0001m= 1.2635 1.0088 1.0066 1.0109 1.0018 0.2635 0.0088 0.0066 0.0109 0.0018 0.2635 0.0088 0.0066 0.0109 0.00180.6439 0.0205 0.0028 0.0006 0.0001%估计积分 真值为0.8862M=f11输入一元被积函数如x.*sin(x):exp(-x.2)g1 =exp(-x.2)输入积分下界a:0输入积分上界b:

9、+inf积分真值:pi0.5/2%0.8862输入样本容量 10V1-10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000Real = 0.8862EY1 = 0.9333 0.9077 0.8873 0.8871D1 = 0.0470 0.0215 0.0010 0.0009RD1 = 0.0531 0.0243 0.0012 0.0010MSE1 = 0.1927 0.0112 0.0016 0.0000M = 0.9333 0.9077 0.8873 0.8871 0.0470 0.0215 0.0010 0.0009 0.0531 0.0243 0.0

10、012 0.0010 0.1927 0.0112 0.0016 0.0000第一类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if x2+y2<=1 I=1;else I=0;end%保存为I2.m%第一类二重积分程序主体%保存为f12.mfunction outf12=f12()g2=input('输入二元被积函数如exp(x.2+y.2):','s')%输入被积函数g2=inline(g2,'x','y');Real=input('积分真值:');

11、%输入积分真值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=10mfor 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(x.2+y.2)/2);Y2(j)=y2*t2'/n; %单次实验样本均值

12、endt=ones(1,10);EY=Y2*t'/10; %十次均值D=abs(EY-Real); %绝对误差RD=D/Real; %绝对误差d=0;for i=1:10 d=d+(Y2(i)-Real)2;endd=d/(10-1);EY2(m-V(1)+1)=EY; %样本容量为10m时的样本均值D2(m-V(1)+1)=D; %绝对误差RD2(m-V(1)+1)=RD; %绝对误差MSE2(m-V(1)+1)=d; %方差endReal,EY2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; %存放样本数字特征%保存为f12.m运行结果：%估计积分 ，真值为

13、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输入样本容量 10V1*10V1-10V2*10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000Real = 5.3981EY2 = 4.7702 5.1250 5.4317 5.4041D2 = 0.6279 0.2732 0.0335 0.0060RD2 = 0.1163 0.0506 0.0062 0.0011MSE2 =3.8965 0.5564 0.0

14、247 0.0017m = 4.7702 5.1250 5.4317 5.4041 0.6279 0.2732 0.0335 0.0060 0.1163 0.0506 0.0062 0.0011 3.8965 0.5564 0.0247 0.0017第二类一重积分程序代码：%构造示性函数function I=I1(x,a,b)if x>=a&&x<=b I=1;else I=0;end%保存为I1.m%第二类一重积分程序主体%程序保存为f21.mfunction outf21=f21()g1=input('输入一元被积函数如exp(x.2):',&#

15、39;s')%输入被积函数g1=inline(g1);a=input('输入积分下界a:');%输入积分上下限b=input('输入积分上界b:');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=10mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量en

16、dy1=g1(x)*(pi*2)0.5).*exp(x.2/2);Y1(j)=y1*t1'/n; %单次实验样本均值endt=ones(1,10);EY=Y1*t'/10; %十次均值d=0;for i=1:10 d=d+(Y1(i)-EY)2;endd=d/(10-1);EY1(m-V(1)+1)=EY; %样本容量为10m时的样本均值MSE1(m-V(1)+1)=d; %方差endEY1,MSE1outf21=EY1;MSE1; %存放样本数字特征%程序保存为f21.m运行结果：%估计积分 m=f21输入一元被积函数如exp(x.2):exp(x.2)g1 =exp(x.2

17、)输入积分下界a:0输入积分上界b:1输入样本容量 10V1-10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000EY1 = 2.0782 1.6583 1.5029 1.4590MSE1 = 0.4315 0.0889 0.0057 0.0008m = 2.0782 1.6583 1.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

18、)S =1.4627第二类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if x2+y2<=1 I=1;else I=0;end%保存为I2.m%第二类二重积分函数主体%，程序保存为f22.mfunction outf22=f22()g2=input('输入二元被积函数如1./(1+x.4+y.4).0.5:','s')%输入被积函数g2=inline(g2,'x','y');fprintf('输入样本容量 10V1*10V1-10V2*10V2:r')

19、V=zeros(1,2);V(1)=input('V1:');%输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量10m1-10m2 n=10mfor 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(x.2+y.2)/2);Y2(j)=y2*t2'/n; %单次实验样本均值endt=ones(1,10);EY=Y2*t'/10; %十次均值d=0;for i=1:10 d

20、=d+(Y2(i)-EY)2;endd=d/(10-1);EY2(m-V(1)+1)=EY; %样本容量为10m时的样本均值MSE2(m-V(1)+1)=d; %方差endEY2,MSE2outf22=EY2;MSE2; %存放样本数字特征%第二类二重积分，程序保存为f22.m运行结果： %估计积分 m=f22输入二元被积函数如1./(1+x.4+y.4).0.5:1./(1+x.4+y.4).0.5g2 =1./(1+x.4+y.4).0.5输入样本容量 10V1*10V1-10V2*10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000EY2 = 3.075

21、9 2.9699 2.8566 2.8269MSE2 = 1.3267 0.0900 0.0060 0.0014m = 3.0759 2.9699 2.8566 2.82691.3267 0.0900 0.0060 0.0014 实验结果整理：第一类一重积分：估计积分 积分真值：1积分估计值：1.0018样本容量：10 100 1000 10000 100000样本均值：1.2635 1.0088 1.0066 1.0109 1.0018绝对误差：0.2635 0.0088 0.0066 0.0109 0.0018相对误差：0.2635 0.0088 0.0066 0.0109 0.0018均

22、方误差：0.6439 0.0205 0.0028 0.0006 0.0001估计积分积分真值：0.8862积分估计值：0.8871样本容量：10 100 1000 10000 样本均值：0.9333 0.9077 0.8873 0.8871绝对误差：0.0470 0.0215 0.0010 0.0009相对误差：0.0531 0.0243 0.0012 0.0010均方误差：0.1927 0.0112 0.0016 0.0000第一类二重积分：估计积分积分真值：5.3981积分估计值： 5.4041样本容量：10 100 1000 10000 样本均值：4.7702 5.1250 5.4317

23、 5.4041绝对误差： 0.6279 0.2732 0.0335 0.0060相对误差：0.1163 0.0506 0.0062 0.0011均方误差：3.8965 0.5564 0.0247 0.0017第二类一重积分：估计积分 积分估计值：1.4590样本容量：10 100 1000 10000 样本均值：2.0782 1.6583 1.5029 1.4590样本方差：0.4315 0.0889 0.0057 0.0008用matlab 指令求得积分结果1.4627第二类二重积分：估计积分积分估计值：2.8269样本容量：10 100 1000 10000 样本均值：3.0759 2.9

24、699 2.8566 2.8269样本方差：1.3267 0.0900 0.0060 0.0014实验结果分析：从第一类积分看，以估计积分 为例：积分真值：1积分估计值：1.0018样本容量：10 100 1000 10000 100000样本均值：1.2635 1.0088 1.0066 1.0109 1.0018绝对误差：0.2635 0.0088 0.0066 0.0109 0.0018相对误差：0.2635 0.0088 0.0066 0.0109 0.0018均方误差：0.6439 0.0205 0.0028 0.0006 0.0001随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误差、相对误差、均方误差呈减小趋势；随着样本容量的增大，样本均值有接近