用蒙特卡洛方法估计积分方法及matlab编程实现专业班级：材料43学生姓名：王宏辉学号：2140201060指导教师：李耀武完成时间：2016年6月8日用蒙特卡洛方法估计积分 方法及matlab 编程实现实验内容：1用蒙特卡洛方法估计积分 20sin x xdx π⎰，2-0x e dx +∞⎰和22221xy x y e dxdy ++≤⎰⎰的值，并将估计值与真值进行比较。

2用蒙特卡洛方法估计积分 210x e dx ⎰和22x y +≤⎰⎰的值，并对误差进行估计。

要求：（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；（2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。

目的：（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；（2） 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息；（3） 能用 MATLAB 熟练进行样本的一元回归分析。

实验原理：蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。

具体操作如下：一般地，积分⎰=bdx x g a )(S 改写成⎰⎰==bb dx f h dx f g aa )(x )(x )(x f(x))(x S 的形式，（其中为)f(x 一随机变量X 的概率密度函数，且)f(x 的支持域）（｝｛b f ,a 0)(x |x ⊇>)，f(x ))(x )(x g h =）;令Y=h(X)，则积分S=E （Y ）；利用matlab 软件，编程产生随机变量X 的随机数，在由⎩⎨⎧∉∈==)b (a,,)b (a,,01I(x) ,)(x )(x y x x I h ，得到随机变量Y 的随机数，求出样本均值，以此估计积分值。

积分⎰⎰=Adxdy g S )y (x,的求法与上述方法类似，在此不赘述。

概率密度函数的选取：一重积分，由于要求)f(x 的支持域）（｝｛b f ,a 0)(x |x ⊇>，为使方法普遍适用，考虑到标准正态分布概率密度函数22e 21)(x xf -=π支持域为R ，故选用22e 21)(x x f -=π。

类似的，二重积分选用22221)y (x,y x e f +-=π，支持域为2R 。

估计评价：进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度。

程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。

为了使程序具有一般性以及方便以后使用：一重积分，程序保存为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存为一个.m文本，被积函数键盘输入，示性函数用function 语句构造，求不同区域二重积分，只需改变function 函数内容。

编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。

程序代码及运行结果：第一类一重积分程序代码：%%%构造示性函数function I=I1(x,a,b)if x>=a&&x<=bI=1;elseI=0;%保存为I1.m%%%%%%%%%%%%%%%%%%第一类一重积分,程序主体：%保存为f11.mfunction 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:');for m=V(1):V(2)%样本容量10^m1--10^m2n=10^mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量y1=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:10d=d+(Y1(i)-Real)^2;endd=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; %方差endReal,EY1,D1,RD1,MSE1outf11=[EY1;D1;RD1;MSE1]; %存放样本数字特征%保存为f11.m运行结果：%估计积分2sinx xdxπ⎰，积分真值为1 m=f11输入一元被积函数如x.*sin(x):x.*sin(x) g1 =x.*sin(x)输入积分下界a:0输入积分上界b:pi/2积分真值:1输入样本容量10^V1--10^V2:V1:1V2:5n =10n =100n =1000 n =10000 n =100000 Real =1EY1 =1.2635 1.0088 1.0066 1.0109 1.0018 D1 =0.2635 0.0088 0.0066 0.0109 0.0018 RD1 =0.2635 0.0088 0.0066 0.0109 0.0018 MSE1 =0.6439 0.0205 0.0028 0.0006 0.0001m=1.2635 1.0088 1.0066 1.0109 1.00180.2635 0.0088 0.0066 0.0109 0.00180.2635 0.0088 0.0066 0.0109 0.00180.6439 0.0205 0.0028 0.0006 0.0001%估计积分 2-0x edx +∞⎰真值为0.8862M=f11输入一元被积函数如x.*sin(x):exp(-x.^2) g1 =exp(-x.^2)输入积分下界a:0 输入积分上界b:+inf 积分真值:pi^0.5/2%0.8862 输入样本容量 10^V1--10^V2:V1:1V2:4n =10n =100n =1000 n =10000Real =0.8862EY1 =0.9333 0.9077 0.8873 0.8871 D1 =0.0470 0.0215 0.0010 0.0009 RD1 =0.0531 0.0243 0.0012 0.0010 MSE1 =0.1927 0.0112 0.0016 0.0000 M =0.9333 0.9077 0.8873 0.88710.0470 0.0215 0.0010 0.00090.0531 0.0243 0.0012 0.00100.1927 0.0112 0.0016 0.0000第一类二重积分程序代码：%%%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if x^2+y^2<=1I=1;elseI=0;end%保存为I2.m%第一类二重积分程序主体%保存为f12.mfunction outf12=f12()g2=input('输入二元被积函数如exp(x.^2+y.^2):','s')%输入被积函数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('V2:');for m=V(1):V(2)%样本容量10^m1--10^m2n=10^mfor 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=abs(EY-Real); %绝对误差RD=D/Real; %绝对误差d=0;for i=1:10d=d+(Y2(i)-Real)^2;endd=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; %方差 endReal,EY2,D2,RD2,MSE2outf12=[EY2;D2;RD2;MSE2]; %存放样本数字特征 %保存为f12.m 运行结果： %估计积分 22221xy x y e dxdy ++≤⎰⎰，真值为pi*(exp(1)-1)%5.3981m=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:1 V2:4 n =10n =100n =1000 n =10000 Real =5.3981EY2 =4.77025.1250 5.4317 5.4041 D2 =0.6279 0.2732 0.0335 0.0060RD2 =0.1163 0.0506 0.0062 0.0011 MSE2 =3.8965 0.5564 0.0247 0.0017m =4.77025.1250 5.4317 5.40410.6279 0.2732 0.0335 0.00600.1163 0.0506 0.0062 0.00113.8965 0.5564 0.0247 0.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.^2):','s')%输入被积函数g1=inline(g1);a=input('输入积分下界a:');%输入积分上下限b=input('输入积分上界b:');fprintf('输入样本容量10^V1--10^V2:\r')V=zeros(1,2);V(1)=input('V1:');%输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量10^m1--10^m2n=10^mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量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=0;for i=1:10d=d+(Y1(i)-EY)^2;endd=d/(10-1);EY1(m-V(1)+1)=EY; %样本容量为10^m时的样本均值MSE1(m-V(1)+1)=d; %方差endEY1,MSE1outf21=[EY1;MSE1]; %存放样本数字特征%%%%程序保存为f21.m 运行结果：%估计积分21xe dxm=f21输入一元被积函数如exp(x.^2):exp(x.^2) g1 =exp(x.^2)输入积分下界a:0输入积分上界b:1输入样本容量10^V1--10^V2:V1:1V2:4n =10n =100n =1000n =10000EY1 =2.0782 1.6583 1.5029 1.4590 MSE1 =0.4315 0.0889 0.0057 0.0008m =2.0782 1.6583 1.5029 1.45900.4315 0.0889 0.0057 0.0008%用matlab 指令求积分21xe dxf=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)if x^2+y^2<=1I=1;elseI=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('输入样本容量10^V1*10^V1--10^V2*10^V2:\r')V=zeros(1,2);V(1)=input('V1:');%输入样本容量V(2)=input('V2:');for m=V(1):V(2)%样本容量10^m1--10^m2n=10^mfor 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:10d=d+(Y2(i)-EY)^2;endd=d/(10-1);EY2(m-V(1)+1)=EY; %样本容量为10^m时的样本均值MSE2(m-V(1)+1)=d; %方差endEY2,MSE2outf22=[EY2;MSE2]; %存放样本数字特征%第二类二重积分，程序保存为f22.m运行结果：%估计积分⎰⎰22+≤x ym=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输入样本容量10^V1*10^V1--10^V2*10^V2:V1:1V2:4n =10n =100n =1000n =10000EY2 =3.0759 2.9699 2.8566 2.8269 MSE2 =1.3267 0.0900 0.0060 0.0014m =3.0759 2.9699 2.8566 2.82691.3267 0.0900 0.0060 0.0014实验结果整理：第一类一重积分：估计积分2sinx xdxπ⎰积分真值：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估计积分2-0x e dx +∞⎰积分真值：0.8862 积分估计值：0.8871 样本容量：101001000 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第一类二重积分： 估计积分22221xy x y e dxdy ++≤⎰⎰积分真值：5.3981 积分估计值： 5.4041 样本容量：101001000 10000样本均值：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第二类一重积分： 估计积分 21x e dx ⎰积分估计值：1.4590 样本容量：101001000 10000样本均值：2.0782 1.6583 1.5029 1.4590 样本方差：0.4315 0.0889 0.0057 0.0008 用matlab 指令求得积分结果1.4627 第二类二重积分：估计积分22x y +≤⎰⎰积分估计值：2.8269 样本容量：101001000 10000样本均值：3.0759 2.9699 2.8566 2.8269 样本方差：1.3267 0.0900 0.0060 0.0014实验结果分析：从第一类积分看，以估计积分2sinx xdxπ⎰为例：积分真值：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 随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误差、相对误差、均方误差呈减小趋势；随着样本容量的增大，样本均值有接近积分真值的趋势，说明估计具有无偏性；绝对误差、相对误差、均方误差呈减小趋势，说明增大样本容量能提高估计精度；验证了蒙特卡洛方法估计积分值的可行性，为后续估计第二类积分提供了参考。