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

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

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

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

4、為一個(gè).m文本，被積函數(shù)鍵盤輸入，示性函數(shù)用function 語句構(gòu)造，求不同區(qū)域二重積分，只需改變function 函數(shù)內(nèi)容。編程完整解決用蒙特卡洛方法估計(jì)一重、二重積分值問題。程序代碼及運(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=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; %單次實(shí)驗(yà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時(shí)的樣本均值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; %存放樣本數(shù)字特征%保存為f11.m運(yùn)行結(jié)果：%估計(jì)積分 ，積分真值為1m=f11輸入一元被積函數(shù)如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%估計(jì)積分 真值為0.8862M=f11輸入一元被積函數(shù)如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第一類二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)function I=I2(x,y)if x2+y2=a&x S=

10、quadl(f,0,1)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)V=zeros(1,2);V(1)=input(V1:);%輸入樣本容量V(2)=in

11、put(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; %單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y2*t/10; %十次均值d=0;for i=1:10 d=d+(Y2(i)-EY)2;endd=d/(10-1);EY2(m-V(1)+1)=EY; %樣本容量為10m時(shí)的樣本均值MSE2(m-V(1)

12、+1)=d; %方差endEY2,MSE2outf22=EY2;MSE2; %存放樣本數(shù)字特征%第二類二重積分，程序保存為f22.m運(yùn)行結(jié)果： %估計(jì)積分 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.0759 2.9699 2.8566 2.8269MSE2 = 1.3267 0.0900 0.0060 0.0014m = 3.0759 2.9699

13、 2.8566 2.82691.3267 0.0900 0.0060 0.0014 實(shí)驗(yàn)結(jié)果整理：第一類一重積分：估計(jì)積分 積分真值：1積分估計(jì)值：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估計(jì)積分積分真值：0.8862積分估計(jì)值：0.8871樣本容量：10

14、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第一類二重積分：估計(jì)積分積分真值：5.3981積分估計(jì)值： 5.4041樣本容量：10 100 1000 10000 樣本均值：4.7702 5.1250 5.4317 5.4041絕對誤差： 0.6279 0.2732 0.0335 0.0060相對誤差：0.1163 0.0506 0.0062 0.0011均

15、方誤差：3.8965 0.5564 0.0247 0.0017第二類一重積分：估計(jì)積分 積分估計(jì)值：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第二類二重積分：估計(jì)積分積分估計(jì)值：2.8269樣本容量：10 100 1000 10000 樣本均值：3.0759 2.9699 2.8566 2.8269樣本方差：1.3267 0.0900 0.0060 0.0014實(shí)驗(yàn)結(jié)果分析：從第一類積分看，以估計(jì)積分 為例：

16、積分真值：1積分估計(jì)值：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隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，絕對誤差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，說明估計(jì)具有無偏性；絕對誤差、相對誤差、均方誤差呈減小

17、趨勢，說明增大樣本容量能提高估計(jì)精度；驗(yàn)證了蒙特卡洛方法估計(jì)積分值的可行性，為后續(xù)估計(jì)第二類積分提供了參考。從第二類積分看，以估計(jì)積分 為例：積分估計(jì)值：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由于積分真值未知，無法直接比較估計(jì)值與積分值值；但隨樣本容量增大，樣本方差減小，間接反映了估計(jì)精度的提高。蒙特卡洛方法估計(jì)值1.4590相比用matlab 指令求得的積分結(jié)果1.4627，絕對偏差0.0038，相對偏差0.0025。蒙特卡洛方法估計(jì)值與用matlab 指令求得的積分結(jié)果相互驗(yàn)證?？偨Y(jié)與討論：蒙特卡洛方法是基于隨機(jī)數(shù)的一種統(tǒng)計(jì)方法。蒙特卡洛方法估計(jì)積分值，總的思想是將積分改寫為某個(gè)隨機(jī)變量的數(shù)學(xué)期望，借助相應(yīng)的隨機(jī)數(shù)，利用樣本均值估計(jì)數(shù)學(xué)期望，從而估計(jì)相應(yīng)的積分值。為使方法具有一般性，概率密度函數(shù)一重積分選擇了，二重積分選用。程序設(shè)計(jì)方面，本著使程序具有一般性以及方便以后使用的