《用蒙特卡洛方法估计积分方法及matlab编程实现.》由会员分享，可在线阅读，更多相关《用蒙特卡洛方法估计积分方法及matlab编程实现.（32页珍藏版）》请在金锄头文库上搜索。
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=inline(g1);a=input(输入积分下界a:);%输入积分上下限b=in
5、put(输入积分上界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);%示性及求和向量endy1=g1(x)*(pi*2)0.5).*exp(x.2/2);Y1(j)=y1*t1/n; %单次实验样本均值endt=ones(1,10);EY=Y1*t/10; %
6、十次均值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,EY1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; %存放样本数字特征%保存为f11.m运行结果：%估计积分 ，积分真值为1m=f11输入一元被积函数如x.*sin(x):x.*sin(
7、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.0018RD1 = 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.
8、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:+inf积分真值:pi0.5/2%0.8862输入样本容量 10V1-10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000Real = 0.8862EY1 = 0.9333 0.9
9、077 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.0012 0.0010 0.1927 0.0112 0.0016 0.0000第一类二重积分程序代码：%构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if x2+y2=1 I=1;
10、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(积分真值:);%输入积分真值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
11、(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: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; %绝对误差MS
12、E2(m-V(1)+1)=d; %方差endReal,EY2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; %存放样本数字特征%保存为f12.m运行结果：%估计积分 ，真值为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.0247 0.0017m = 4.7702 5.1250 5.4317 5.4041 0.6279 0.2732 0.0335 0.006
