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

1、用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)專業(yè)班級：材料43學(xué)生姓名：王宏輝學(xué) 號：2140201060指導(dǎo)教師：李耀武完成時間：2016年6月8日用蒙特卡洛方法估計積分方法及matlab編程實現(xiàn)實驗內(nèi)容：22 一221用家特卡洛方法估計積分xsinxdx, e-xdx和ex y dxdy的值，00x2 y2 1并將估計值與真值進行比較。2用蒙特卡洛方法估計積分 1ex2dx和, 11dxdy的值,并對誤0x2 y2 1 1 x4 y4差進行估計。要求：(1)針對要估計的積分選擇適當?shù)母怕史植荚O(shè)計蒙特卡洛方法；(2)利用計算機產(chǎn)生所選分布的隨機數(shù)以估計積分值；(3)進行重復(fù)試驗，通過計算

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

3、aa f(x)a式，(其中為f(x)一隨機變量 X的概率密度函數(shù)，且f(x)的支持域x|f(x) 0 (a,b), h(x) 幽)；令丫牛兇，則積分 S=E(Y);利用 f(x)matlab軟件，編程產(chǎn)生隨機變量 X的隨機數(shù)，在由y h(x)I(x), I(x) 1 ,x (a,b),得到隨機變量Y的隨機數(shù)，求出樣本均 0 ,x (a,b)值，以此估計積分值。積分S g (x, y) dxdy的求法與上述方法類似，在此不贅述。A概率密度函數(shù)的選?。?(a,b),為使方法普一重積分，由于要求f(x)的支持域x| f(x)x2遍適用，考慮到標準正態(tài)分布概率密度函數(shù)f(x) 1 e萬支持域為R ,故

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

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

6、ut(' 輸入積分上界b:');Real=input(' 積分真值:');% 輸入積分真值fprintf('輸入樣本容量 1071-10"2:')V=zeros(1,2);V(1)=input('V1:');% 輸入樣本容量V(2)=input('V2:');for m!=丫(1)»出 樣本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);% 示性及求和向量endy1=g1(x)*(pi*2)八0.5).

7、*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;for i=1:10d=d+(Y1(i)-Real2;endd=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10 Am時的樣本均值D1(m-V(1)+1)=D;% 絕對誤差RD1(m-V(1)+1)=RD;% 絕對誤差MSE1(m-V(1)+1)=d; % 方差endReal,EY1,D1,RD1,MSE1outf11=EY1;

8、D1;RD1;MSE1; % 存放樣本數(shù)字特征%保存為f11.m歡迎下載運行結(jié)果:%估計積分xsin xdx ,積分真值為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 =歡迎下載1EY1 =歡迎下載1.26351.00881.00661.01091.0018D1 =0.26350.00880.00660.01090.0018RD1 =0.26350.00880.00660.0

9、1090.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.00012%估計積分e-x dx真值為0.88620M=f11輸入一元被積函數(shù)如x.*sin(x):exp(-x.A2)g1 =exp(-x.A2)輸入積分下界a:0輸入積分上界b:+inf積分真值:piA0.5/2%0.8862輸入樣本容量10AV1-10AV2:V1:1

10、歡迎下載V2:4歡迎下載10100100010000Real =0.8862EY1 =0.93330.90770.88730.8871D1 =0.04700.02150.00100.0009RD1 =0.05310.02430.00120.0010MSE1 =歡迎下載0.19270.01120.00160.00000.93330.90770.04700.02150.05310.02430.19270.01120.88730.88710.00100.00090.00120.00100.00160.0000歡迎下載第一類二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)%funct

11、ion I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存為I2.m%第一類二重積分程序主體%保存為f12.mfunction outf12=f12()g2=input('輸入二元被積函數(shù)如exp(x.A2+y.A2):','s')%輸入被積函數(shù)g2=inline(g2,'x','y');Real=input(' 積分真值:');% 輸入積分真值fprintf(' 輸入樣本容量10AV1*10AV1-10AV2*10AV2:r')V=zeros(1,2);V(1)

12、=input('V1:');% 輸入樣本容量V(2)=input('V2:');for m=V(1):V(2)% 樣本容量10Am1-10Am2n=10Amfor 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.A2+y.A2)/2);Y2(j)=y2*t2'/n; % 單次實驗樣本均值endt=ones(1,10);EY=Y2*t'/10; % 十次均值D=abs(EY-Real); % 絕對誤差RD

13、=D/Real; % 絕對誤差d=0;for i=1:10d=d+(Y2(i)-Real2;endd=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為10 Am時的樣本均值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運行結(jié)果：%估計積分 ex2 y2dxdy ,真值為 pi*(exp(1)-1)%5.3981 x2 y2 1m=f12輸入二元被積函數(shù)如

14、exp(x.A2+y.A2):exp(x.A2+y.A2)g2 =exp(x.A2+y.A2)積分真值:pi*(exp(1)-1)%5.3981輸入樣本容量10AV1*10AV1-10AV2*10AV2:V1:1V2:4n =10歡迎下載100100010000Real =5.3981歡迎下載EY2 =4.7702 5.1250 5.4317D2 =0.6279 0.2732 0.0335RD2 =0.1163 0.0506 0.0062MSE2 =3.8965 0.5564 0.0247m =4.7702 5.1250 5.43170.6279 0.2732 0.03350.1163 0.0

15、506 0.00623.8965 0.5564 0.02475.40410.00600.00110.00175.40410.00600.00110.0017歡迎下載第二類一重積分程序代碼：%構(gòu)造示性函數(shù)%function I=I1(x,a,b)if x>=a&&x<=bI=1;elseI=0;end%保存為I1.m%第二類一重積分程序主體%程序保存為f21.mfunction outf21=f21()輸入被積函g1=input('輸入一元被積函數(shù)如exp(x.A2):','s')%數(shù)g1=inline(g1);a=input('

16、;輸入積分下界a:');% 輸入積分上下限b=input('輸入積分上界b:');fprintf('輸入樣本容量10AV1-10AV2:r')V=zeros(1,2);V(1)=input('V1:');% 輸入樣本容量V(2)=input('V2:');for m=V(1):V(2)% 樣本容量 10Am110Am2 n=10Amfor j=1:10 x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);% 示性及求和向量endy1=g1(x)*(pi*2)A0.5).*exp(x.A2/2)

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

18、1輸入樣本容量10AV1-10AV2:V1:1V2:4n =n =10歡迎下載100100010000EY1 =2.07821.65831.50291.4590MSE1 =0.43150.08890.00570.0008歡迎下載2.07821.65831.50291.45900.43150.08890.00570.000812%用 matlab 指令求積分ex dx0f=inline('exp(x.A2)')f =Inline function:f(x) = exp(x.A2)>> S=quadl(f,0,1)S =1.4627第二類二重積分程序代碼：%構(gòu)造示性函數(shù)

19、，求不同區(qū)域上積分只需更改示性函數(shù)%function I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存為I2.m%第二類二重積分函數(shù)主體%，程序保存為f22.mfunction outf22=f22()g2=input('輸入二元被積函數(shù)如1./(l+x.M+y.M)八0.5:','s')%輸入被積函數(shù)g2=inline(g2,'x','y');fprintf('輸入樣本容量 1071*1071-1072*10"2:丁)V=zeros(1,2);V(1)=input('

20、;V1:');% 輸入樣本容量V(2)=input('V2:');for m=V(1):V(2)% 樣本容量 10Am1-10Am2n=10Amfor 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.A2+y.A2)/2);Y2(j)=y2*t2'/n; % 單次實驗樣本均值 endt=ones(1,10);EY=Y2*t'/10; % 十次均值d=0;for i=1:10d=d+(Y2(i)-EY2;end

21、d=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為10 Am時的樣本均值MSE2(m-V(1)+1)=d; % 方差 endEY2,MSE2outf22=EY2;MSE2; % 存放樣本數(shù)字特征%第二類二重積分，程序保存為f22.m運行結(jié)果:%估計積分 1dxdyx2 y21 .1 x4 y4m=f22輸入二元被積函數(shù)如 1./(1+x.A4+y.M)八0.5:1./(1+x.M+y.A4)八0.5 g2 =1./(1+x.A4+y.A4).A0.5輸入樣本容量 10AV1*10AV1-10AV2*10AV2:V1:1V2:4 n = 10 n = 100100010000E

22、Y2 =3.07592.96992.85662.8269MSE2 =1.32670.09000.00600.0014歡迎下載3.07592.96992.85662.82691.32670.09000.00600.0014實驗結(jié)果整理: 第一類一重積分:估計積分2xsin xdx0積分真值:積分估計值:1.0018樣本容量:10100100010000樣本均值:1.26351.00881.00661.01091000001.0018絕對誤差：0.26350.00880.00660.01090.0018相對誤差：0.2635均方誤差：0.64392估計積分e-x dx0積分真值：0.8862樣本容

23、量：10樣本均值：0.9333絕對誤差：0.0470相對誤差：0.0531均方誤差：0.19270.00880.00660.02050.00280.01090.00180.00060.0001積分估計值：10010000.90770.88730.02150.00100.02430.00120.01120.00160.8871100000.88710.00090.00100.0000第一類二重積分：估計積分x222e y dxdy y2 1積分真值：5.3981積分估計值：5.4041樣本容量：10100100010000樣本均值：4.77025.1250 5.4317 5.4041絕對誤差：0

24、.6279 0.2732 0.0335 0.0060相對誤差：0.11630.0506 0.0062 0.0011均方誤差：3.8965 0.5564 0.0247 0.0017第二類一重積分:12估計積分ex dx0積分估計值：1.4590樣本容量：10100100010000樣本均值：2.07821.65831.50291.4590樣本方差：0.43150.08890.00570.0008用matlab指令求得積分結(jié)果1.4627第二類二重積分:估計積分f 1 dxdyx2 y2 1 J x4 y4積分估計值：2.8269樣本容量：10100100010000樣本均值：3.07592.96

25、992.85662.8269樣本方差：1.32670.09000.00600.0014實驗結(jié)果分析：22從第一類積分看，以估計積分xsin xdx為例:0積分真值：1積分估計值：1.0018樣本容量：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隨著樣本容量的增大，樣本均值有接近積分真值的趨勢，絕對誤 差、相對誤差、均方誤差呈減小趨勢；隨著樣本容量的增大