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

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

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

3、，且f(x)的支持域x|f(x)O=( a,b) , h(x-gx)；令 Y=h(X)，则积分 S=E( Y；利用 f(x)matlab软件，编程产生随机变量X的随机数，在由 y = h(x)l(x), l(x) =公匸，得到随机变量丫的随机数，求出样本均0 X(a,b)值，以此估计积分值。积分S = . g(x, y)dxdy的求法与上述方法类似，在此不赘述。A概率密度函数的选取：一重积分，由于要求f(x)的支持域x|f(x)0二(a,b)，为使方法普1故选用f(xF彳 廿类似的，二重积分选用f(x, y) e 2，支持域为R22兀估计评价:进行重复试验，通过计算样本均值以评价估计的无偏性；

4、通过计 算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积 不出）以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一 重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存 为一个.m文本，被积函数，积分区间均采用键盘输入；二重积分， 程序主体保存为一个.m文本，被积函数键盘输入，示性函数用 fun ction 语句构造，求不同区域二重积分，只需改变 fun ction 函数内容。编程完整解决用蒙特卡洛方法估计一重、二重积分值问题。程序代码及运行结果：第一类一重积分程序代码：%构造示性函数fun ct

5、io n I=l1（x,a,b）if x=a&x=bl=1;elsel=0;end%保存为I1.m%第一类一重积分，程序主体：%保存为f11.mfun ctio n outf11=f11()输入被积g1=input(输入一元被积函数如 x.*sin(x):,s)% 函数g1=i nlin e(g1);a=input(输入积分下界a:);%输入积分上下限b=input(输入积分上界b:);Real二input(积分真值:);%输入积分真值 fprintf( 输入样本容量 10AV1-10AV2:r) V=zeros(1,2);V(1)=i nput(V1:);% 输入样本容量V(2)=i npu

6、t(V2:);for m=V(1):V(2)% 样本容量 10八口1-10八口2n=10八口for j=1:10x=ra ndn (1, n);for i=1: nt1(i)=l1(x(i),a,b);% 示性及求和向量endy仁 g1(x)*(pi*2)A0.5).*exp(x.A2/2);Y 1(j)=y1*t1/n; % 单次实验样本均值end t=o nes(1,10);EY二Y 1*t/10; % 十次均值D=abs(EY-Real); % 绝对误差RD二D/Real; % 绝对误差d=0;for i=1:10d=d+(Y 1(i)-Real)A2;endd=d/(10-1);EY

7、1(m-V(1)+1)=E Y; % 样本容量为10八m时的样本均值D1(m-V(1)+1)=D; % 绝对误差RD1(m-V(1)+1)=RD; % 绝对误差MSE1(m-V(1)+1)=d; % 方差endReal,E Y1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; % 存放样本数字特征%保存为f11.m运行结果：2%古计积分 xs in xdx，积分真值为10m=f11输入一元被积函数如x.*sin(x):x.*sin(x)g1 = x.*si n(x)输入积分下界a:0输入积分上界b:pi/2积分真值:1输入样本容量10AV1-10AV2:V1:1V2:5n

8、 二10n 二100n 二1000n 二10000n 二100000Real =1EY1 =1.2635 1.0088 1.0066 1.0109 1.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.000

9、60.0001%估计积分 e-x2dx真值为0.88620M=f11输入一元被积函数如x.*si n(x):exp(-x42)g1 =exp(-x.A2)输入积分下界a:0输入积分上界b:+inf积分真值：piA0.5/2%0.8862输入样本容量10AV1-10AV2:V1:1V2:4n 二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.93

10、330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.1927 0.0112 0.00160.0000第一类二重积分程序代码:%造示性函数，求不同区域上积分只需更改示性函数fun cti on I=l2(x,y)if xA2+yA2=1l=1;elsel=0;end%保存为I2.m%第一类二重积分程序主体%保存为f12.mfun ctio n outf12二f12()g2=in put(输入二元被积函数如 exp(x.A2+y.A2):,s)%入被积函数g2=i nlin e(g2,x,y);Real二inp

11、ut(积分真值:);%输入积分真值fprintf( 输入样本容量 10AV1*10AV1-10AV2*10AV2:r)V=zeros(1,2);V(1)=i nput(V1:);% 输入样本容量V(2)=i nput(V2:);for m=V(1):V(2)% 样本容量 10八口1-10八口2n=10八口for j=1:10x=ra ndn (1, n);y=ra ndn (1, n);for i=1: nt2(i)=l2(x(i),y(i);% 示性及求和向量endy2=g2(x,y)*(2*pi).*exp(x.A2+y.A2)/2);Y 2(j)=y2*t2/n; % 单次实验样本均值e

12、ndt=o nes(1,10);EY二Y 2*t/10; % 十次均值D=abs(EY-Real); % 绝对误差RD二D/Real; % 绝对误差d=0;for i=1:10d=d+( Y2(i)-RealF2;endd=d/(10-1);EY 2(m-V(1)+1)=E Y;%样本容量为10a m时的样本均值D2(m-V(1)+1)=D;%绝对误差RD2(m-V(1)+1)=RD;%绝对误差MSE2(m-V(1)+1)=d;%方差Real,E Y2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; % 存放样本数字特征%保存为f12.m运行结果：%古计积分 ex2 y2dxdy，真值为 pi*(exp(1)-1)%5.3981xy2=a&x=bI=1;elseI=0;%保存为I1.m%第二类一重积分程序主体%程序保存为f21.mfun