用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)

1、 用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn) 專業(yè)班級： 材料43學(xué)生姓名： 王宏輝學(xué) 號： 2140201060指導(dǎo)教師：李耀武完成時間： 2016年6月8日 用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)實驗內(nèi)容：1用蒙特卡洛方法估計積分 ，和的值，并將估計值與真值進(jìn)行比較。2用蒙特卡洛方法估計積分 和的值，并對誤差進(jìn)行估計。要求： （1）針對要估計的積分選擇適當(dāng)?shù)母怕史植荚O(shè)計蒙特卡洛方法；（2）利用計算機(jī)產(chǎn)生所選分布的隨機(jī)數(shù)以估計積分值；（3）進(jìn)行重復(fù)試驗，通過計算樣本均值以評價估計的無偏性；通過計算均方誤差（針 對第1類題）或樣本方差（針對第2類題）以評價估計結(jié)果的精度。目的：

2、 （1）能通過 MATLAB 或其他數(shù)學(xué)軟件了解隨機(jī)變量的概率密度、分布函數(shù) 及其期望、方差、協(xié)方差等；（2） 熟練使用 MATLAB 對樣本進(jìn)行基本統(tǒng)計，從而獲取數(shù)據(jù)的基本信息；（3） 能用 MATLAB 熟練進(jìn)行樣本的一元回歸分析。實驗原理： 蒙特卡洛方法估計積分值，總的思想是將積分改寫為某個隨機(jī)變量的數(shù)學(xué)期望，借助相應(yīng)的隨機(jī)數(shù)，利用樣本均值估計數(shù)學(xué)期望，從而估計相應(yīng)的積分值。具體操作如下：一般地，積分改寫成的形式，（其中為一隨機(jī)變量X的概率密度函數(shù)，且的支持域)，）;令Y=h(X)，則積分S=E（Y）；利用matlab軟件，編程產(chǎn)生隨機(jī)變量X的隨機(jī)數(shù)，在由，得到隨機(jī)變量Y的隨機(jī)數(shù)，求出樣

3、本均值，以此估計積分值。積分的求法與上述方法類似，在此不贅述。概率密度函數(shù)的選?。阂恢胤e分，由于要求的支持域，為使方法普遍適用，考慮到標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)支持域為，故選用。類似的，二重積分選用，支持域為。估計評價：進(jìn)行重復(fù)試驗，通過計算樣本均值以評價估計的無偏性；通過計算均方誤（針對第1類題，積得出）或樣本方差（針對第2類題，積不出）以評價估計結(jié)果的精度。程序設(shè)計：依據(jù)問題分四類：第一類一重積分；第一類二重積分；第二類一重積分，第二類二重積分，相應(yīng)程序設(shè)計成四類。為了使程序具有一般性以及方便以后使用：一重積分，程序保存為一個.m文本，被積函數(shù)，積分區(qū)間均采用鍵盤輸入；二重積分，程序主體保存

4、為一個.m文本，被積函數(shù)鍵盤輸入，示性函數(shù)用function 語句構(gòu)造，求不同區(qū)域二重積分，只需改變function 函數(shù)內(nèi)容。編程完整解決用蒙特卡洛方法估計一重、二重積分值問題。程序代碼及運(yùn)行結(jié)果：第一類一重積分程序代碼：%構(gòu)造示性函數(shù)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('輸入一元被積函數(shù)如x.*sin(x):','s')%輸入被積函數(shù)g1=inli

5、ne(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-10m2 n=10mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b)

6、;%示性及求和向量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: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

7、,EY1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1; %存放樣本數(shù)字特征%保存為f11.m運(yùn)行結(jié)果：%估計積分 ，積分真值為1m=f11輸入一元被積函數(shù)如x.*sin(x):x.*sin(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.001

8、8RD1 = 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.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輸入一元被積函數(shù)如x.*sin(x):exp(-x.2)g1 =exp(-x.2)輸入積分下界a:0輸入積分上界b:

9、+inf積分真值:pi0.5/2%0.8862輸入樣本容量 10V1-10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000Real = 0.8862EY1 = 0.9333 0.9077 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.0

10、012 0.0010 0.1927 0.0112 0.0016 0.0000第一類二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)function I=I2(x,y)if x2+y2<=1 I=1;else I=0;end%保存為I2.m%第一類二重積分程序主體%保存為f12.mfunction outf12=f12()g2=input('輸入二元被積函數(shù)如exp(x.2+y.2):','s')%輸入被積函數(shù)g2=inline(g2,'x','y');Real=input('積分真值:');

11、%輸入積分真值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(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; %單次實驗樣本均值

12、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; %絕對誤差MSE2(m-V(1)+1)=d; %方差endReal,EY2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; %存放樣本數(shù)字特征%保存為f12.m運(yùn)行結(jié)果：%估計積分 ，真值為

13、pi*(exp(1)-1)%5.3981 m=f12輸入二元被積函數(shù)如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.0

14、247 0.0017m = 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第二類一重積分程序代碼：%構(gòu)造示性函數(shù)function I=I1(x,a,b)if x>=a&&x<=b I=1;else I=0;end%保存為I1.m%第二類一重積分程序主體%程序保存為f21.mfunction outf21=f21()g1=input('輸入一元被積函數(shù)如exp(x.2):',&#

15、39;s')%輸入被積函數(shù)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)%樣本容量10m1-10m2 n=10mfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量en

16、dy1=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:10 d=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;MSE1; %存放樣本數(shù)字特征%程序保存為f21.m運(yùn)行結(jié)果：%估計積分 m=f21輸入一元被積函數(shù)如exp(x.2):exp(x.2)g1 =exp(x.2

17、)輸入積分下界a:0輸入積分上界b:1輸入樣本容量 10V1-10V2:V1:1V2:4n = 10n = 100n = 1000n = 10000EY1 = 2.0782 1.6583 1.5029 1.4590MSE1 = 0.4315 0.0889 0.0057 0.0008m = 2.0782 1.6583 1.5029 1.4590 0.4315 0.0889 0.0057 0.0008%用matlab 指令求積分 f=inline('exp(x.2)')f =Inline function: f(x) = exp(x.2)>> S=quadl(f,0,1

18、)S =1.4627第二類二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)function I=I2(x,y)if x2+y2<=1 I=1;else I=0;end%保存為I2.m%第二類二重積分函數(shù)主體%，程序保存為f22.mfunction outf22=f22()g2=input('輸入二元被積函數(shù)如1./(1+x.4+y.4).0.5:','s')%輸入被積函數(shù)g2=inline(g2,'x','y');fprintf('輸入樣本容量 10V1*10V1-10V2*10V2:r')

19、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(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:10 d

20、=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; %存放樣本數(shù)字特征%第二類二重積分，程序保存為f22.m運(yùn)行結(jié)果： %估計積分 m=f22輸入二元被積函數(shù)如1./(1+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.075

21、9 2.9699 2.8566 2.8269MSE2 = 1.3267 0.0900 0.0060 0.0014m = 3.0759 2.9699 2.8566 2.82691.3267 0.0900 0.0060 0.0014 實驗結(jié)果整理：第一類一重積分：估計積分 積分真值：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均

22、方誤差：0.6439 0.0205 0.0028 0.0006 0.0001估計積分積分真值：0.8862積分估計值：0.8871樣本容量：10 100 1000 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第一類二重積分：估計積分積分真值：5.3981積分估計值： 5.4041樣本容量：10 100 1000 10000 樣本均值：4.7702 5.1250 5.4317

23、 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第二類一重積分：估計積分 積分估計值：1.4590樣本容量：10 100 1000 10000 樣本均值：2.0782 1.6583 1.5029 1.4590樣本方差：0.4315 0.0889 0.0057 0.0008用matlab 指令求得積分結(jié)果1.4627第二類二重積分：估計積分積分估計值：2.8269樣本容量：10 100 1000 10000 樣本均值：3.0759 2.9

24、699 2.8566 2.8269樣本方差：1.3267 0.0900 0.0060 0.0014實驗結(jié)果分析：從第一類積分看，以估計積分 為例：積分真值：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隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，絕對誤差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大，樣本均值有接近