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

1、精品用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)專業(yè)班級：材料43學生姓名：王宏輝學號：2140201060指導教師：李耀武完成時間：2016年6月8日感謝下載載用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)實驗內容：i用蒙特卡洛方法估計積分2xsin xdxo2 一22e-xdx 和ex ydxdy 的值,0x2 y2 1并將估計值與真值進行比較2用蒙特卡洛方法估計積分12exdx和0差進行估計。要求：(1)針對要估計的積分選擇適當?shù)母怕史植荚O計蒙特卡洛方法；(2)利用計算機產(chǎn)生所選分布的隨機數(shù)以估計積分值；(3)進行重復試驗，通過計算樣本均值以評價估計的無偏性；通過計算均方誤差(針對

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

3、x|f(x)0(a,b),h(x)嗎);4Y=h(X),則積分S=E(Y);利f(x)用matlab軟件，編程產(chǎn)生隨機變量X的隨機數(shù)，在由yh(x)I(x),I(x)1,x(；：)，得到隨機變量Y的隨機數(shù)，求出樣本均值，以此估計積分值。積分Sg(x,y)dxdy的求法與上述方法類似，在此不贅述。A概率密度函數(shù)的選取:0(a,b),為使方法普一重積分，由于要求f(x)的支持域x|f(x)遍適用，考慮到標準正態(tài)分布概率密度函數(shù)1與f (x) 2je e 支持域為1x2R，故選用f(x)7re2。22支持域為R2oxy類似的，二重積分選用f(x,y)2e2精品估計評價：進行重復試驗，通過計算樣本均值

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

5、nI=I1（x,a,b）ifx=a&x二bI=1;elseend%保存為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-102:)V=zeros(1,2);V(1)=input(V1:);%輸入樣本容量V(2)=input(V2:);form!=丫(1)出樣本容量10Am1-

6、10Am2n=10Amforj=1:10x=randn(1,n);fori=1:nt1(i)=I1(x(i),a,b);%示性及求和向量endy1=g1(x)*(pi*2)八0.5).*exp(x.A2/2);Y1(j)=y1*t1/n;%單次實驗樣本均值endt=ones(1,10);EY=Y1*t/10;%十次均值D=abs(EY-Real);%絕對誤差RD=D/Real;%絕對誤差d=0;fo門=1:10d=d+(Y1(i)-Real)A2;endd=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10Am時的樣本均D1(m-V(1)+1)=D;%絕對誤差RD1(m-V(

7、1)+1)=RD;%絕對誤差MSE1(m-V(1)+1)=d;%方差endReal,EY1,D1,RD1,MSE1outf11=EY1;D1;RD1;MSE1;%存放樣本數(shù)字特征%保存為f11.m運行結果：2%估計積分xsinxdx,積分真值為10m=f11輸入一元被積函數(shù)如x.*sin(x):x.*sin(x)g1=x.*sin(x)輸入積分下界a:0輸入積分上界b:pi/2積分真值：1輸入樣本容量10AV1-10AV2:V1:1V2:5n=10感謝下載載100100010000100000Real=EY1=1.26351.00881.00661.0109D1=0.26350.00880.0

8、0660.0109RD1=0.26350.00880.00660.0109MSE1=0.64390.02050.00280.00061.00180.00180.00180.0001m=1.01091.00181.26351.00881.00660.26350.00880.00660.01090.00180.26350.00880.00660.01090.00180.64390.02050.00280.00060.00012%估計積分e-xdx真值為0.88620M=f11輸入一元被積函數(shù)如x.*sin(x):exp(-xJ2)g1=exp(-x.A2)輸入積分下界a:0輸入積分上界b:+inf

9、積分真值鐘0.5/2%0.8862輸入樣本容量10AV1-10AV2:V1:1V2:4n=10n=100n=1000n=10000精品MSE1 =感謝下載載Real=0.8862EY1=0.93330.90770.88730.8871D1=0.04700.02150.00100.0009RD1=0.05310.02430.00120.0010精品0.19270.01120.00160.0000M=0.93330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.19270.01120.00160.0000感謝下載

10、載第一類二重積分程序代碼：%構造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)functionI=I2(x,y)ifxA2+yA2=a&xS=quadl(f,0,1)S=1.4627感謝下載載第二類二重積分程序代碼：%構造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)functionI=I2(x,y)ifxA2+yA2=1I=1;elseI=0;end%保存為I2.m%第二類二重積分函數(shù)主體%,程序保存為f22.mfunctionoutf22=f22()g2=input(輸入二元被積函數(shù)如1./(1+乂八4+丫八4)八0.5:,&)%輸入被積函數(shù)g2=inline(g2,x,y);fprintf(輸入

11、樣本容量10A丫1*10八丫1-10八丫2*10八丫2：)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);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;%單次實驗樣本均值endt=ones(1,10);EY=Y2*t/10;%十次均值d=0;fo門=1:10d=d+(Y

12、2(i)-EY)A2;endd=d/(10-1);時的樣本均值EY2(m-V(1)+1)=EY;%樣本容量為10AmMSE2(m-V(1)+1)=d;%方差endEY2,MSE2outf22=EY2;MSE2;%存放樣本數(shù)字特征%第二類二重積分，程序保存為f22.m運行結果：%估計積分,1=dxdyx2y2i,1x4y4m=f22輸入二元被積函數(shù)如1./(1+x.A4+y.A4).A0.5:1./(1+x.A4+y.A4).A0.5g2=1./(1+x.A4+yA4)A0.5輸入樣本容量10A丫1*10八丫1-10八丫2*10八丫2:V1:1V2:4n=10n=精品100感謝下載載100010

13、000EY2=3.07592.96992.85662.8269MSE2=1.32670.09000.00600.0014精品3.07592.96991.32670.09002.85662.82690.00600.0014感謝下載載實驗結果整理：第一類一重積分：2估計積分xsinxdx0積分真值：1積分估計值：1.0018樣本容量：101001000100001000001.26351.00881.00661.01090.26350.00880.00660.01090.26350.00880.00660.01090.64390.02050.00280.00062e-xdx0.8862積分估計值：

14、0.8871101001000100000.93330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.19270.01120.00160.0000樣本均值:絕對誤差:相對誤差:均方誤差:估計積分0積分真值:樣本容量:樣本均值:絕對誤差:相對誤差:均方誤差:1.00180.00180.00180.0001第一類二重積分:估計積分x222exydxdyy21積分真值:5.3981積分估計值:5.4041樣本容量:10100100010000樣本均值:4.77025.12505.43175.4041絕對誤差:0.6

15、2790.27320.03350.0060相對誤差:0.11630.05060.00620.0011均方誤差：3.89650.55640.02470.0017100001.45900.0008100002.82690.0014第二類一重積分：1-2估計積分exdx0積分估計值：1.4590樣本容量：101001000樣本均值：2.07821.65831.5029樣本方差：0.43150.08890.0057用matlab指令求得積分結果1.4627第二類二重積分:估計積分1-dxdyx2y21,1x4y4積分估計值：2.8269樣本容量：101001000樣本均值：3.07592.96992.

16、8566樣本方差：1.32670.09000.0060實驗結果分析:2從第一類積分看，以估計積分xsinxdx為例:0積分真值：1樣本容量：10樣本均值：1.2635絕對誤差：0.2635相對誤差：0.2635均方誤差：0.6439積分估計值：1.001810010001.00881.00660.00880.00660.00880.00660.02050.0028100001000001.01091.00180.01090.00180.01090.00180.00060.0001隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，絕對誤差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大，樣本均

17、值有接近積分真值的趨勢，說明估計具有無偏性；絕對誤差、相對誤差、均方誤差呈減小趨勢，說明增大樣本容量能提高估計精度；驗證了蒙特卡洛方法估計積分值的可行性，為后續(xù)估計第二類積分提供了1一、一一一2-從第二類積分看，以估計積分exdx為例:0積分估計值：1.4590樣本容量:樣本均值:樣本方差:102.07820.43151001.65830.088910001.50290.0057100001.45900.0008用matlab指令求得積分結果1.4627由于積分真值未知，無法直接比較估計值與積分值值；但隨樣本容量增大，樣本方差減小，間接反映了估計精度的提高。蒙特卡洛方法估計值1.4590相比用matlab指令求得的積分結果1.4627,絕對偏差0.0038,相對偏差0.0025。蒙特卡洛方法估計值與用matlab指令求得的積分結果相互驗證。總結與討論：蒙特卡洛方法是基于隨機數(shù)的一種統(tǒng)計方法。蒙特卡洛方法估計積分值，總的思想是將積分改寫為某個隨機變量的數(shù)學期望，借助相應的隨機數(shù)，利用樣本均值估計數(shù)學期望，從而估計相應的積分值。為使方法具有一般性，概率密度函數(shù)一重積分選擇了1x21xLif(x)-y=e2,二重積分選用f(x,y)2e2。程序設計方面，本著使程序具有一般性以及方便以后使用的原則，依據(jù)問題分四類：第一類一重積分；第一類二重積分；第二類一重積分，