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

1、用蒙特卡洛方法估計(jì)積分方法及matlab編程實(shí)現(xiàn)專業(yè)班級：材料43學(xué)生姓名：王宏輝學(xué)號：2140202160指導(dǎo)教師：李耀武完成時(shí)間：2021年6月8日用蒙特卡洛方法估計(jì)積分方法及matlab編程實(shí)現(xiàn)實(shí)驗(yàn)內(nèi)容：2,一221用家特卡洛方法估計(jì)積分xsinxdx,e-xdx和exydxdy的值,00x2y21并將估計(jì)值與真值進(jìn)行比擬.2用蒙特卡洛方法估計(jì)積分1ex2dx和1dxdy的值,并對誤0x2y211x4y4差進(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ì)的無偏性；通過

2、計(jì)算均方誤差(針對第1類題)或樣本方差(針又t第2類題)以評價(jià)估計(jì)結(jié)果的精度.目的：(1)能通過MATLAB或其他數(shù)學(xué)軟件了解隨機(jī)變量的概率密度、分布函數(shù)及其期望、方差、協(xié)方差等；(2)熟練使用MATLAB對樣本進(jìn)行根本統(tǒng)計(jì),從而獲取數(shù)據(jù)的(3)能用MATLAB熟練進(jìn)行樣本的一元回歸分析.實(shí)驗(yàn)原理：蒙特卡洛方法估計(jì)積分值,總的思想是將積分改寫為某個隨機(jī)變量的數(shù)學(xué)期望,借助相應(yīng)的隨機(jī)數(shù),利用樣本均值估計(jì)數(shù)學(xué)期望,從而估計(jì)相應(yīng)的積分值.具體操作如下：bb()b一一一般地,積分Sg(x)dx改與成Sf(x)dxh(x)f(x)dx的形aaf(x)a式,(其中為f(x)一隨機(jī)變量X的概率密度函數(shù),且f

3、(x)的支持域x|f(x)0(a,b),h(x)幽)；令丫牛兇,那么積分S=E(Y);利用f(x)matlab軟件,編程產(chǎn)生隨機(jī)變量X的隨機(jī)數(shù),在由yh(x)I(x),I(x)1,x(a,b),得到隨機(jī)變量Y的隨機(jī)數(shù),求出樣本均0,x(a,b)值,以此估計(jì)積分值.積分Sg(x,y)dxdy的求法與上述方法類似,在此不贅述.A概率密度函數(shù)的選?。?(a,b),為使方法普一重積分,由于要求f(x)的支持域x|f(x)x2遍適用,考慮到標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)f(x)1e萬支持域?yàn)?1-R,應(yīng)選用f(x)Je2.2支持域?yàn)镽2.八、一13類似的,一重積分選用f(x,y)-e2,估計(jì)評價(jià):進(jìn)行重復(fù)試驗(yàn)

4、,通過計(jì)算樣本均值以評價(jià)估計(jì)的無偏性；通過計(jì)算均方誤針對第1類題,積得出或樣本方差針對第2類題,積不出以評價(jià)估計(jì)結(jié)果的精度.程序設(shè)計(jì)：依據(jù)問題分四類：第一類一重積分；第一類二重積分；第二類一重積分,第二類二重積分,相應(yīng)程序設(shè)計(jì)成四類.為了使程序具有一般性以及方便以后使用：一重積分,程序保存為一個.m文本,被積函數(shù),積分區(qū)間均采用鍵盤輸入；二重積分,程序主體保存為一個.m文本,被積函數(shù)鍵盤輸入,示性函數(shù)用function語句構(gòu)造,求不同區(qū)域二重積分,只需改變function函數(shù)內(nèi)容.編程完整解決用蒙特卡洛方法估計(jì)一重、二重積分值問題.程序代碼及運(yùn)行結(jié)果：第一類一重積分程序代碼：%造示性函數(shù)fun

5、ctionI=I1x,a,bifx>=a&&x<二bI=1;elseI=0;%保存為I1.m%嘴一類一重積分,程序主體：%保存為f11.mfunctionoutf11=f11()g1=input('輸入一元被積函數(shù)如x.*sin(x):','s')%輸入被積函數(shù)g1=inline(g1);a=input('輸入積分下界a:');%輸入積分上下限b=input('輸入積分上界b:');Real=input('積分真值:');輸入積分真值fprintf('輸入樣本容量1071-107

6、2:)V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');form=V(1):V樣本容量10Am1-10Am2n=10Amforj=1:10x=randn(1,n);fori=1:nt1(i)=I1(x(i),a,b);%示性及求和向量y1=g1(x)*(pi*2)A0.5),*exp(x.A2/2);Y1(j)=y1*t1'/n;%單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y1*t'/10;%十次均值D=abs(EY-Real);%絕對誤差RD=D/Real;%絕對誤差d=

7、0;fori=1:10d=d+(Y1(i)-Real)A2;endd=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10Am時(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ù)字特征運(yùn)行結(jié)果:%估計(jì)積分xsinxdx,積分真值為10m=f11輸入一元被積函數(shù)如x.*sin(x):x.*sin(x)gi=x.*sin(x)輸入積分下界a:0輸入積分上界b:pi/2積分真值：1輸入樣本容量10A

8、V1-10AV2:V1:1V2:5n=100100010000100000Real=EY1=1.26351.00881.00661.01091.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.00060.0001

9、一一一2%估計(jì)積分e-xdx真值為0.88620M=f11輸入一元被積函數(shù)如x.*sin(x):exp(-x.A2)gi=exp(-x.A2)輸入積分下界a:0輸入積分上界b:+inf積分真值：piA0.5/2%0.8862輸入樣本容量10AV1-10AV2:10100100010000Real=0.8862EY1=0.93330.90770.88730.8871D1=0.04700.02150.00100.0009RD1=0.05310.02430.00120.00100.19270.01120.00160.00000.93330.90770.04700.02150.05310.02430.

10、19270.01120.88730.88710.00100.00090.00120.00100.00160.0000第一類二重積分程序代碼:%造示性函數(shù),求不同區(qū)域上積分只需更改示性函數(shù)functionI=I2(x,y)ifxA2+yA2<=1I=1;elseI=0;end%保存為I2.m施一類二重積分程序主體%保存為f12.mfunctionoutf12=f12()g2=input('輸入二元被積函數(shù)如exp(x.A2+y.A2):','s')%輸入被積函數(shù)g2=inline(g2,'x','y');Real=input(

11、'積分真值:');輸入積分真值fprintf('輸入樣本容量10AV1*10AV1-10AV2*10AV2:r')V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');form=V(1):V樣本容量10Am1-10Am2forj=1:10x=randn(1,n);y=randn(1,n);fori=1:nt2(i)=I2(x(i),y(i);%示性及求和向量endy2=g2(x,y)*(2*pi).*exp(x.A2+y.A2)/2);Y2(j)=y2*t2'/n;

12、%單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y2*t'/10;%十次均值D=abs(EY-Real);%絕對誤差RD=D/Real;%絕對誤差d=0;fori=1:10d=d+(Y2(i)-Real)A2;endd=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為10Am時(shí)的樣本均值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é)果：%估計(jì)積分ex2y

13、2dxdy,真值為pi*(exp(1)-1)%5.3981x2y21m=f12輸入二元被積函數(shù)如exp(x,A2+y.A2):exp(x,A2+y.A2)g2=exp(x.A2+yJ2)積分真值:pi*(exp(1)-1)%5.3981輸入樣本容量10AV1*10AV1-10AV2*10AV2:V1:1V2:4n=100100010000Real=EY2=4.77025.12505.4317D2=0.62790.27320.0335RD2=0.11630.05060.0062MSE2=3.89650.55640.0247m=4.77025.12505.43170.62790.27320.033

14、50.11630.05060.00623.89650.55640.02475.40410.00600.00110.00175.40410.00600.00110.0017第二類一重積分程序代碼：%造示性函數(shù)functionI=I1(x,a,b)ifx>=a&&x<二bI=1;elseI=0;end%保存為I1.m施二類一重積分程序主體%程序保存為f21.mfunctionoutf21=f21()輸入被積函g1=input('輸入一元被積函數(shù)如exp(x.A2):','s')%數(shù)g1=inline(g1);a=input(,輸入積分下界

15、a:');%輸入積分上下限b=input('輸入積分上界b:');fprintf('輸入樣本容量10AV1-10AV2:丁)V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');form=V(1):V(2)%樣本容量10Am1-10Am2n=10Amforj=1:10x=randn(1,n);fori=1:nt1(i)=I1(x(i),a,b);%示性及求和向量endy1=g1(x)*(pi*2)A0.5),*exp(x.A2/2);Y1(j)=y1*t1'/n;%

16、單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y1*t'/10;%十次均值d=0;fori=1:10d=d+(Y1(i)-EY)A2;endd=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10Am時(shí)的樣本均值MSE1(m-V(1)+1)=d;%方差endEY1,MSE1outf21=EY1;MSE1;%存放樣本數(shù)字特征%程序保存為f21.m運(yùn)行結(jié)果：1、一2%估計(jì)積分exdx0m=f21輸入一元被積函數(shù)如exp(x.A2):exp(x.A2)gi=exp(x.A2)輸入積分下界a:0輸入積分上界b:1輸入樣本容量10AV1-10AV2:V1:1V2:4n=1

17、00100010000EY1=2.07821.65831.50291.4590MSE1=0.43150.08890.00570.00081.45900.00082.07821.65831.50290.43150.08890.005712%用matlab指令求積分exdx0f=inline('exp(x.A2)')f=Inlinefunction:f(x)=exp(x.A2)>>S=quadl(f,0,1)1.4627第二類二重積分程序代碼：%造示性函數(shù),求不同區(qū)域上積分只需更改示性函數(shù)functionI=I2(x,y)ifxA2+yA2<=1I=1;elseI

18、=0;end%保存為I2.m施二類二重積分函數(shù)主體%程序彳存為f22.mfunctionoutf22=f22()g2=input('輸入二元被積函數(shù)如1./(1+乂八4+丫八4)八0.5:',&')%輸入被積函數(shù)g2=inline(g2,'x','y');fprintf('輸入樣本容量10AV1*10AV1-10AV2*10AV2:丁)V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');form=V(1):V(2)%樣本容量10Am

19、1-10Am2n=10Amforj=1:10x=randn(1,n);y=randn(1,n);fori=1:nt2(i)=I2(x(i),y(i);%示性及求和向量endy2=g2(x,y)*(2*pi).*exp(x.A2+y.A2)/2);Y2(j)=y2*t2'/n;%單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y2*t'/10;%十次均值d=0;fori=1:10d=d+(Y2(i)-EY)A2;endd=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為10Am時(shí)的樣本均值MSE2(m-V(1)+1)=d;%方差endEY2,MSE2outf

20、22=EY2;MSE2;%存放樣本數(shù)字特征施二類二重積分,程序保存為f22.m運(yùn)行結(jié)果:%估計(jì)積分1dxdyx2y21Jix,/m=f22輸入二元被積函數(shù)如1./(1+xA4+yA4)A0.5:1./(1+xA4+yA4)A0.5g2=1./(1+xA4+yA4)A0.5輸入樣本容量10AV1*10AV1-10AV2*10AV2:V1:1V2:4n=10n=100100010000EY2=3.07592.96992.85662.8269MSE2=3.07592.96991.32670.09002.85662.82690.00600.0014實(shí)驗(yàn)結(jié)果整理：第一類一重積分：2估計(jì)積分xsinxdx

21、0積分真值：1積分估計(jì)值：1.0018樣本容量：10100100010000100000樣本均值：1.26351.00881.00661.01091.0018絕對誤差：0.26350.00880.00660.01090.0018相對誤差：0.2635均方誤差：0.6439、一一2估計(jì)積分e-xdx0積分真值：0.8862樣本容量：10樣本均值：0.9333絕對誤差：0.0470相對誤差：0.0531均方誤差：0.19270.00880.00660.02050.00280.01090.00180.00060.0001積分估計(jì)值：10010000.90770.88730.02150.00100.0

22、2430.00120.01120.00160.8871100000.88710.00090.00100.0000第一類二重積分:估計(jì)積分exydxdyx2y21積分真值：5.3981積分估計(jì)值：5.4041樣本容量：10100100010000樣本均值：4.77025.12505.43175.4041絕對誤差：0.62790.27320.03350.0060相對誤差：0.11630.05060.00620.0011土勻方誤差：3.89650.55640.02470.0017第二類一重積分:12估計(jì)積分exdx0積分估計(jì)值：1.4590樣本容量：10100100010000樣本均值：2.0782

23、1.65831.50291.4590樣本方差：0.43150.08890.00570.0008用matlab指令求得積分結(jié)果1.4627第二類二重積分:估計(jì)積分f1dxdy51x4y4積分估計(jì)值：2.8269樣本容量：10100100010000樣本均值：3.07592.96992.85662.8269樣本方差：1.32670.09000.00600.0014實(shí)驗(yàn)結(jié)果分析:2從第一類積分看,以估計(jì)積分xsinxdx為例:0樣本容量：10樣本均值：1.2635絕對誤差：0.2635相對誤差：0.2635均方誤差：0.643910010001.00881.00660.00880.00660.00880.00660.02050.0028100001000001.01091.00180.01090.00180.01090.00180.00060.0001隨著樣本容量的增大,樣本均值有接近積分真值的趨勢,絕對誤差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大,樣本均值有接近積分真值的趨勢,說明估計(jì)具有無偏性；絕對誤差、相對誤差、均方誤差呈