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

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

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

3、，(其中為f(x) 一隨機(jī)變量 X的概率密度函數(shù)，且f(x)的支持域x|f(x)>0；( a,b) ,)；令 Y=h(X)，則積分 S=E( Y；利用f(x)matlab軟件，編程產(chǎn)生隨機(jī)變量X的隨機(jī)數(shù)，在由 y = h(x)l(x), l(x)公匸b)，得到隨機(jī)變量丫的隨機(jī)數(shù)，求出樣本均0 ,x(a,b)值，以此估計積分值。積分S = . g(x, y)dxdy的求法與上述方法類似，在此不贅述。A概率密度函數(shù)的選?。阂恢胤e分，由于要求f(x)的支持域x|f(x)0二(a,b)，為使方法普遍適用，考慮到標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)f(x)= 1 J2支持域為常2兀R，故選用 f (x) =

4、= e2x2O類似的，二重積分選用f(x,y)=e，支持域為R2。估計評價:進(jìn)行重復(fù)試驗，通過計算樣本均值以評價估計的無偏性；通過計 算均方誤（針對第1類題，積得出）或樣本方差（針對第2類題，積 不出）以評價估計結(jié)果的精度。程序設(shè)計：依據(jù)問題分四類：第一類一重積分；第一類二重積分；第二類一 重積分，第二類二重積分，相應(yīng)程序設(shè)計成四類。為了使程序具有一般性以及方便以后使用：一重積分，程序保存 為一個.m文本，被積函數(shù)，積分區(qū)間均采用鍵盤輸入；二重積分， 程序主體保存為一個.m文本，被積函數(shù)鍵盤輸入，示性函數(shù)用 fun ction 語句構(gòu)造，求不同區(qū)域二重積分，只需改變fun ction 函數(shù)內(nèi)容

5、。編程完整解決用蒙特卡洛方法估計一重、二重積分值問題。程序代碼及運(yùn)行結(jié)果：第一類一重積分程序代碼：%構(gòu)造示性函數(shù)fun ctio n I=l1（x,a,b）if x>=a&&x<=bl=1;elsel=0;end%保存為I1.m%第一類一重積分，程序主體：%保存為f11.mfun ctio n outf11=f11()輸入被積g1=input('輸入一元被積函數(shù)如 x.*sin(x):','s')%函數(shù)g1=i nlin e(g1);a=input('輸入積分下界a:');%輸入積分上下限b=input('輸入

6、積分上界b:');Real二input('積分真值:');%輸入積分真值 fprintf(' 輸入樣本容量 10AV1-10AV2:r') V=zeros(1,2);V(1)=i nput('V1:');% 輸入樣本容量V(2)=i nput('V2:');for m=V(1):V(2)% 樣本容量 10八口1-10八口2n=10八口for j=1:10x=ra ndn (1, n);for i=1: nt1(i)=I1(x(i),a,b);%示性及求和向量endy仁g1(x)*(pi*2F0.5).*exp(x.A2/2

7、);Y1(j)=y1*t17 n; %單次實(shí)驗樣本均值endt=o nes(1,10);EY二Y 1*t'/10;% 十次均值D=abs(EY-Real); % 絕對誤差RD二D/Real; % 絕對誤差d=0;for i=1:10d=d+( Y1(i)-RealF2;endd=d/(10-1);EY 1(m-V(1)+1)=E Y;% 樣本容量為10八m時的樣本均值D1(m-V(1)+1)=D;%絕對誤差RD1(m-V(1)+1)=RD;% 絕對誤差MSE1(m-V(1)+1)=d;% 方差endReal,E Y1,D1,RD1,MSE1outfl仁EY 1;D1;RD1;MSE1;

8、 %存放樣本數(shù)字特征%保存為f11.m運(yùn)行結(jié)果:2%古計積分xsinxdx,積分真值為10m=f11輸入一元被積函數(shù)如x.*sin(x):x.*sin(x)g1 = x.*si n(x)輸入積分下界a:0輸入積分上界b:pi/2積分真值:1輸入樣本容量10AV1-10AV2:V1:1V2:5100100010000100000Real =EY1 =1.26351.00881.00661.01091.0018D1 =0.26350.00880.00660.01090.0018RD1 =0.26350.00880.00660.01090.0018MSE1 =0.64390.02050.00280.

9、00060.0001m=1.26351.00881.00661.01091.00180.26350.00880.00660.01090.00180.26350.00880.00660.01090.00180.64390.02050.00280.00060.0001-be2%估計積分.e-x dx真值為0.88620M=f11輸入一元被積函數(shù)如x.*si n(x):exp(-x42)g1 =exp(-x.A2)輸入積分下界a:0輸入積分上界b:+inf積分真值:piA0.5/2%0.8862輸入樣本容量10AV1-10AV2:V1:1V2:410100100010000Real =0.8862E

10、Y1 =0.93330.90770.88730.8871D1 =0.04700.02150.00100.0009RD1 =0.05310.02430.00120.00100.19270.01120.00160.0000M =0.93330.90770.88730.88710.04700.02150.00100.00090.05310.02430.00120.00100.19270.01120.00160.0000第一類二重積分程序代碼：%造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)fun cti on I=l2(x,y)if xA2+yA2<=1l=1;elsel=0;end%保存為I2

11、.m%第一類二重積分程序主體%保存為f12.mfun ctio n outf12二f12()g2=in put('輸入二元被積函數(shù)如exp(x.A2+y.A2):','s')%入被積函數(shù)g2=i nlin e(g2,'x','y');Real二input('積分真值:');%輸入積分真值fprintf(' 輸入樣本容量 10AV1*10AV1-10AV2*10AV2:r')V=zeros(1,2);V(1)=i nput('V1:');%輸入樣本容量V(2)=i nput('

12、V2:');for m=V(1):V(2)% 樣本容量 10八口1-10八口2for j=1:10x=ra ndn (1, n);y=ra ndn (1, n);for i=1: nt2(i)=l2(x(i),y(i);%示性及求和向量endy2=g2(x,y)*(2*pi).*exp(x.A2+y.A2)/2);Y 2(j)=y2*t2'/n; %單次實(shí)驗樣本均值endt=o nes(1,10);EY二Y 2*t'/10; %十次均值D=abs(EY-Real); % 絕對誤差RD二D/Real;%絕對誤差d=0;for i=1:10d=d+( Y2(i)-RealF

13、2;endd=d/(10-1);EY 2(m-V(1)+1)=E Y;%樣本容量為10a m時的樣本均值D2(m-V(1)+1)=D;%絕對誤差RD2(m-V(1)+1)=RD;%絕對誤差endReal,E Y2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; % 存放樣本數(shù)字特征%保存為f12.m運(yùn)行結(jié)果：%古計積分ex2 y2dxdy，真值為 pi*(exp(1)-1)%5.3981x2 -ydm=f12輸入二元被積函數(shù)如 exp(x.A2+y.A2):exp(x.A2+y.A2)g2 =exp(x.A2+y.A2)積分真值:pi*(exp(1)-1)%5.3981輸

14、入樣本容量 10AV1*10AV1-10AV2*10AV2:V1:1V2:4100100010000Real =5.3981EY2 =4.77025.12505.43175.4041D2 =0.62790.27320.03350.0060RD2 =0.1163MSE2 =0.05060.00620.00113.89650.55640.02470.0017m =4.77025.12505.43175.40410.62790.27320.03350.00600.11630.05060.00620.00113.89650.55640.02470.0017第二類一重積分程序代碼：%構(gòu)造示性函數(shù)fun

15、ctio n I=l1(x,a,b)if x>=a&&x<=bl=1;elsel=0;end%保存為I1.m%第二類一重積分程序主體%5序保存為f21.mfun ctio n outf21二f21()輸入被積函g仁input('輸入一元被積函數(shù)如exp(x.A2):','s')%數(shù)g1= in li ne(g1);a=input('輸入積分下界a:');%輸入積分上下限b=input('輸入積分上界b:');fprintf('輸入樣本容量 10AV1-10AV2:r')V=zeros(1

16、,2);V(1)=i nput('V1:');% 輸入樣本容量V(2)=i nput('V2:');for m=V(1):V(2)% 樣本容量 10Am1-10Am2n=10八口for j=1:10x=ra ndn (1, n);for i=1: nt1(i)=l1(x(i),a,b);%示性及求和向量endy1=g1(x)*(pi*2)A0.5).*exp(x.A2/2);Y1(j)=y1*t1'/n; %單次實(shí)驗樣本均值endt=o nes(1,10);EY二Y 1*t'/10; %十次均值d=0;for i=1:10d=d+( Y1(i)-

17、E 丫)八2;endd=d/(10-1);EY 1(m-V(1)+1)=E Y;% 樣本容量為10八m時的樣本均值MSE1(m-V(1)+1)=d;% 方差endEY 1,MSE1outf2仁EY1;MSE1; %存放樣本數(shù)字特征%序保存為f21.m運(yùn)行結(jié)果：12%古計積分ex dx0m=f21輸入一元被積函數(shù)如exp(x.A2):exp(x.A2)g1 =exp(x42)輸入積分下界a:0輸入積分上界b:1輸入樣本容量10AV1-10AV2:V1:1V2:4100100010000EY1 =2.07821.65831.50291.4590MSE1 =0.43150.08890.00570.0

18、0082.07821.65831.50291.45900.43150.08890.00570.0008%用matlab指令求積分12ex dxof=i nlin e('exp(x.A2)') f =In li ne fun cti on:f(x) = exp(x42)>> S=quadl(f,O,1)S =1.4627第二類二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)fun cti on I=l2(x,y)if xA2+yA2<=1l=1;elsel=0;end%保存為I2.m%第二類二重積分函數(shù)主體%程序保存為f22.mfun ctio

19、 n outf22=f22()g2=in put('輸入二元被積函數(shù)如 1./(1+x.A4+y.A4).A0.5:','s')%入被積函數(shù)g2=i nlin e(g2,'x','y');fprin tf('輸入樣本容量 10AV1*10AV1-10AV2*10AV2:r')V=zeros(1,2);V(1)=i nput('V1:');% 輸入樣本容量V(2)=i nput('V2:');for m=V(1):V(2)% 樣本容量 10八口1-10八口2n=10八口for j=1:

20、10x=ra ndn (1, n);y=ra ndn (1, n);for i=1: nt2(i)=l2(x(i),y(i);%示性及求和向量endy2=g2(x,y)*(2*pi).*exp(x.A2+y.A2)/2);Y 2(j)=y2*t2'/n; %單次實(shí)驗樣本均值endt=o nes(1,10);EY二Y 2*t'/10; %十次均值d=0;for i=1:10d=d+( Y2(i)-E 丫)八2;endd=d/(10-1);EY 2(m-V(1)+1)=E Y;%樣本容量為10八口時的樣本均值MSE2(m-V(1)+1)=d;% 方差endEY 2,MSE2outf

21、22=EY2;MSE2; %存放樣本數(shù)字特征%第二類二重積分，程序保存為f22.m運(yùn)行結(jié)果:%古計積分1x2,211x4 .y4dXdym=f22輸入二元被積函數(shù)如1./(1+x.A4+y.A4).A0.5:1./(1+x.A4+y.A4).A0.5g2 =1./(1+乂.八4+.八4).八0.5輸入樣本容量10AV1*10AV1-10AV2*10AV2:V1:1V2:410n 二100100010000EY2 =3.07592.96992.85662.8269MSE2 =1.32670.09000.00600.00143.07592.96992.85662.82691.32670.09000

22、.00600.0014實(shí)驗結(jié)果整理：第一類一重積分:估計積分2xsin xdx0積分真值:1積分估計值：1.0018樣本容量：10 100 1000 10000 100000樣本均值：1.26351.00881.00661.01091.0018絕對誤差：0.26350.00880.00660.01090.0018相對誤差：0.26350.00880.00660.01090.0018均方誤差：0.64390.02050.00280.00060.0001估計積分e-Xdx0積分真值：0.8862積分估計值：0.8871樣本容量：10100100010000樣本均值：0.93330.90770.88

23、730.8871絕對誤差：0.04700.02150.00100.0009相對誤差：0.05310.02430.00120.0010均方誤差：0.19270.01120.00160.0000第一類二重積分: 估計積分| | e" "dxdyX2為2藝積分真值：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第二類一重積分

24、：1估計積分 ex dx0積分估計值：1.4590樣本容量：10100100010000樣本均值：2.07821.65831.50291.4590樣本方差：0.43150.08890.00570.0008用matlab指令求得積分結(jié)果1.4627第二類二重積分：估計積分"dxdyx2 4yyj1：x4 匚 y4積分估計值：2.8269樣本容量：10100100010000樣本均值：3.07592.96992.85662.8269樣本方差：1.32670.09000.00600.0014實(shí)驗結(jié)果分析:"2從第一類積分看，以估計積分.xs inxdx為例: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隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，絕對誤 差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大，樣本均 值有接近積分真值的趨勢，說