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

1、精品用蒙特卡洛方法估计积分方法及 matlab 编程实现专业班级：材料 43学生姓名：王宏辉学号：2140201060指导教师：李耀武完成时间：2016年6月8日感谢下载载精品用蒙特卡洛方法估计积分方法及 matlab编程实现实验内容：2e- x2ex2 y21 用蒙特卡洛方法估计积分x sin xdx ，dx 和dxdy 的值，00x2 y2 1并将估计值与真值进行比较。1x212 用蒙特卡洛方法估计积分dxdy 的值，并对误e dx 和0x2 y2 1 1 x4y4差进行估计。要求：（ 1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法；（ 2 ）利用计算机产生所选分布的随机数以估计积

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

3、x)dx 改写成Sf (x )dxh(x) f (x ) dx 的形f(x)aaa式，（其中为 f(x ) 一随机变量 X 的概率密度函数，且f(x ) 的支持域 x | f (x ) 0（ a, b）)， h(x )g (x ) ）令Y=h(X)，则积分S=E（ ）；利;Yf(x)用matlab 软 件 ， 编 程 产 生 随 机 变 量X 的随机数，在由yh(x ) I (x ), I(x)1, x(a, b) ，得到随机变量Y 的随机数，求出样本0, x(a, b)均值，以此估计积分值。积分 Sg(x, y)dxdy的求法与上述方法类似，在此不赘述。A概率密度函数的选取：一重积分，由于要

4、求f(x )的支持域f (x )，为使方法普 x |0 （ a, b）1x2遍适用，考虑到标准正态分布概率密度函数f (x)2 支持域为e21x2R ，故选用 f (x )2 。e21x 2y2类似的，二重积分选用 f (x, y)e2，支持域为 R 2 。2感谢下载载精品估计评价：进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误（针对第1 类题，积得出）或样本方差（针对第2 类题，积不出）以评价估计结果的精度。程序设计：依据问题分四类：第一类一重积分；第一类二重积分；第二类一重积分，第二类二重积分，相应程序设计成四类。为了使程序具有一般性以及方便以后使用：一重积分，程序保存为

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

6、函数如x.*sin(x):','s')% 输入被积函数g1=inline(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-10m2n=10

7、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=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;% 样本容量为10m时的样本均值D1(m-V(1)+1)=D;

8、%绝对误差RD1(m-V(1)+1)=RD;% 绝对误差MSE1(m-V(1)+1)=d;%方差endReal,EY1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; %存放样本数字特征%保存为 f11.m感谢下载载精品运行结果：2%估计积分x sin xdx ，积分真值为 10m=f11输入一元被积函数如x.*sin(x):x.*sin(x)g1 =x.*sin(x)输入积分下界 a:0输入积分上界 b:pi/2积分真值 :1输入样本容量10V1-10V2:V1:1V2:5n =10感谢下载载精品n =100n =1000n =10000n =100000Real =感

9、谢下载载精品1EY1 =1.26351.00881.00661.01091.0018D1 =0.26350.00880.00660.01090.0018RD1 =0.26350.00880.00660.01090.0018MSE1 =0.64390.02050.00280.00060.0001感谢下载载精品m=1.26351.00881.00661.01091.00180.26350.00880.00660.01090.00180.26350.00880.00660.01090.00180.64390.02050.00280.00060.00012% 估计积分e- x dx 真值为 0.886

10、20M=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 =10000感谢下载载精品Real =0.8862EY1 =0.93330.90770.88730.8871D1 =0.04700.02150.00100.0009RD1 =0.05310.02430.00120.0010MSE1 =感谢下载载精品0.19270.01120.00160.0000M =0.93330.9

11、0770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.19270.01120.00160.0000感谢下载载精品第一类二重积分程序代码：% 构造示性函数，求不同区域上积分只需更改示性函数function I=I2(x,y)if x2+y2<=1I=1;elseI=0;end%保存为 I2.m%第一类二重积分程序主体%保存为 f12.mfunction outf12=f12()g2=input('输入二元被积函数如exp(x.2+y.2):','s')%输入被积函数g2=inline

12、(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-10m2n=10m感谢下载载精品for j=1:10x=randn(1,n);y=randn(1,n);for i=1:nt2(i)=I2(x(i),y(i);%示性及求和向量e

13、ndy2=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;% 样本容量为10m时的样本均值D2(m-V(1)+1)=D;%绝对误差RD2(m-V(1)+1)=RD;% 绝对误差感谢下载载精品MSE2(m-V(1)+1)=d;%方差endReal,EY

14、2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; %存放样本数字特征%保存为 f12.m运行结果：%估计积分ex2 y2dxdy ，真值为 pi*(exp(1)-1)%5.3981x2y2 1m=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.3981感谢下载载精品EY2 =4.77025.12505.43175.

15、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感谢下载载精品第二类一重积分程序代码：% 构造示性函数function I=I1(x,a,b)if x>=a&&x<=bI=1;elseI=0;end%保存为 I1.m%第二类一重积分程序主体

16、%程序保存为 f21.mfunction outf21=f21()g1=input('输入一元被积函数如exp(x.2):','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)

17、%样本容量 10m1-10m2n=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;%十次均值d=0;for i=1:10d=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;

18、MSE1; %存放样本数字特征感谢下载载精品% 程序保存为 f21.m运行结果：12%估计积分ex dx0m=f21输入一元被积函数如exp(x.2):exp(x.2)g1 =exp(x.2)输入积分下界 a:0输入积分上界 b:1输入样本容量10V1-10V2:V1:1V2:4n =10n =感谢下载载精品100n =1000n =10000EY1 =2.07821.65831.50291.4590MSE1 =0.43150.08890.00570.0008感谢下载载精品m =2.07821.65831.50291.45900.43150.08890.00570.00081%用 matlab

19、 指令求积分ex2dx0f=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 x2+y2<=1I=1;elseI=0;end%保存为 I2.m%第二类二重积分函数主体%，程序保存为 f22.mfunction outf22=f22()g2=input('输入二元被积函数如1./(1+x.4+y.4).0.5:',

20、's')%输入被积函数g2=inline(g2,'x','y');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-10m2n=10m感谢下载载精品for j=1:10x=randn(1,n);y=randn(1,n);for i=1:nt2(i)=I2(x(i),y(i);%示性及求和向量endy2=g2(x,

21、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;%样本容量为 10m时的样本均值MSE2(m-V(1)+1)=d;%方差endEY2,MSE2outf22=EY2;MSE2; %存放样本数字特征%第二类二重积分，程序保存为f22.m感谢下载载精品运行结果：1%估计积分dxdyx2y2 1 1x4y 4m=f22输入二元被积函数如1./(1

22、+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.07592.96992.85662.8269MSE2 =1.32670.09000.00600.0014感谢下载载精品m =3.07592.96992.85662.82691.32670.09000.00600.0014实验结果整理：第一类一重积分：2估计积分x sin xdx0积分真值： 1积分估计值： 1.0018样本容量： 101001

23、00010000100000感谢下载载精品样本均值： 1.26351.00881.00661.01091.0018绝对误差： 0.26350.00880.00660.01090.0018相对误差： 0.26350.00880.00660.01090.0018均方误差： 0.64390.02050.00280.00060.0001- x2估计积分 e dx0积分真值： 0.8862积分估计值： 0.8871样本容量： 10100100010000样本均值： 0.93330.90770.88730.8871绝对误差： 0.04700.02150.00100.0009相对误差： 0.05310.02

24、430.00120.0010均方误差： 0.19270.01120.00160.0000第一类二重积分：估计积分x2y2edxdyx2y2 1积分真值：5.3981积分估计值： 5.4041样本容量：10100100010000样本均值：4.77025.12505.43175.4041绝对误差： 0.62790.27320.03350.0060相对误差：0.11630.05060.00620.0011感谢下载载精品均方误差： 3.89650.55640.02470.0017第二类一重积分：1x2估计积分 e dx0积分估计值： 1.4590样本容量： 10100100010000样本均值： 2

25、.07821.65831.50291.4590样本方差： 0.43150.08890.00570.0008用 matlab 指令求得积分结果 1.4627第二类二重积分：估计积分1dxdyx2 y2 1 1 x4y4积分估计值： 2.8269样本容量： 10100100010000样本均值： 3.07592.96992.85662.8269样本方差： 1.32670.09000.00600.0014感谢下载载精品实验结果分析：2从第一类积分看，以估计积分x sin xdx 为例：0积分真值： 1积分估计值： 1.0018样本容量： 10100100010000100000样本均值： 1.26351.00881.00661.01091.0018绝对误差： 0.26350.00880.00660.01090.0018相对误差： 0.26350.00880.00660.01090.0018均方误差： 0.64390.02050.00280.00060.0001随着样本容量的增大，样本均值有接近积分真值的趋势，绝对误差、相对误差、均方误差呈减小趋势；随着样本容量的增大，样本均值有接近积分真值的趋势，说明估计具有无偏性；绝对误差、相对误差、均方误差呈减小趋势， 说明