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

用蒙特卡洛方法估计积分方法及matlab编程实现专业班级：材料43学生姓名：王宏辉学号：2140201060指导教师：李耀武完成时间：2016年6月8日用蒙特卡洛方法估计积分方法及matlab编程实现实验内容：1用蒙特卡洛方法估计积分20sinxxdx，2-0xedx和22221xyxyedxdy的值，并将估计值与真值进行比较。2用蒙特卡洛方法估计积分210xedx和2244111xydxdyxy的值，并对误差进行估计。要求：（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；（2）利用计算机产生所选分布的随机数以估计积分值；（3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。目的：（1）能通过MATLAB或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等；（2）熟练使用MATLAB对样本进行基本统计，从而获取数据的基本信息；（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。估计评价：进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积不出）以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分，程序主体保存为一个.m文本，被积函数键盘输入，示性函数用function语句构造，求不同区域二重积分，只需改变function函数内容。编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。程序代码及运行结果：第一类一重积分程序代码：%%%构造示性函数functionI=I1(x,a,b)ifx=a&&x=bI=1;elseI=0;end%保存为I1.m%%%%%%%%%%%%%%%%%%第一类一重积分,程序主体：%保存为f11.mfunctionoutf11=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:');form=V(1):V(2)%样本容量10^m1--10^m2n=10^mforj=1:10x=randn(1,n);fori=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=abs(EY-Real);%绝对误差RD=D/Real;%绝对误差d=0;fori=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运行结果：%估计积分20sinxxdx，积分真值为1m=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=1000n=10000n=100000Real=1EY1=1.26351.00881.00661.01091.0018D1=0.26350.00880.00660.01090.0018RD1=0.26350.00880.00660.01090.0018MSE1=0.64390.02050.00280.00060.0001m=1.26351.00881.00661.01091.00180.26350.00880.00660.01090.00180.26350.00880.00660.01090.00180.64390.02050.00280.00060.0001%估计积分2-0xedx真值为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=1000n=10000Real=0.8862EY1=0.93330.90770.88730.8871D1=0.04700.02150.00100.0009RD1=0.05310.02430.00120.0010MSE1=0.19270.01120.00160.0000M=0.93330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.19270.01120.00160.0000第一类二重积分程序代码：%%%构造示性函数，求不同区域上积分只需更改示性函数functionI=I2(x,y)ifx^2+y^2=1I=1;elseI=0;end%保存为I2.m%第一类二重积分程序主体%保存为f12.mfunctionoutf12=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:');form=V(1):V(2)%样本容量10^m1--10^m2n=10^mforj=1:10x=randn(1,n);y=randn(1,n);fori=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;fori=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运行结果：%估计积分22221xyxyedxdy，真值为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:1V2:4n=10n=100n=1000n=10000Real=5.3981EY2=4.77025.12505.43175.4041D2=0.62790.27320.03350.0060RD2=0.11630.05060.00620.0011MSE2=3.89650.55640.02470.0017m=4.77025.12505.43175.40410.62790.27320.03350.00600.11630.05060.00620.00113.89650.55640.02470.0017第二类一重积分程序代码：%%%构造示性函数functionI=I1(x,a,b)ifx=a&&x=bI=1;elseI=0;end%保存为I1.m%第二类一重积分程序主体%程序保存为f21.mfunctionoutf21=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:');form=V(1):V(2)%样本容量10^m1--10^m2n=10^mforj=1:10x=randn(1,n);fori=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;fori=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运行结果：%估计积分210xedxm=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.07821.65831.50291.4590MSE1=0.43150.08890.00570.0008m=2.07821.65831.50291.45900.43150.08890.00570.0008%用matlab指令求积分210xedxf=inline('exp(x.^2)')f=Inlinefunction:f(x)=exp(x.^2)S=quadl(f,0,1