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

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

2、樣本均值以評(píng)價(jià)估計(jì)的無(wú)偏性；通過(guò)計(jì)算均方誤差(針對(duì)第1類(lèi)題)或樣本方差(針又t第2類(lèi)題)以評(píng)價(jià)估計(jì)結(jié)果的精度。目的：(1)能通過(guò)MATLAB或其他數(shù)學(xué)軟件了解隨機(jī)變量的概率密度、分布函數(shù)及其期望、方差、協(xié)方差等；(2)熟練使用MATLAB對(duì)樣本進(jìn)行基本統(tǒng)計(jì)，從而獲取數(shù)據(jù)的基本信息；(3)能用MATLAB熟練進(jìn)行樣本的一元回歸分析。實(shí)驗(yàn)原理：蒙特卡洛方法估計(jì)積分值，總的思想是將積分改寫(xiě)為某個(gè)隨機(jī)變 量的數(shù)學(xué)期望，借助相應(yīng)的隨機(jī)數(shù)，利用樣本均值估計(jì)數(shù)學(xué)期望，從 而估計(jì)相應(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)一隨機(jī)變量 X的概率密度函數(shù)，且f(x)的支持域x|f(x) 0 (a,b), h(x) 幽)；令丫牛兇，則積分 S=E(Y);利用 f(x)matlab軟件，編程產(chǎn)生隨機(jī)變量 X的隨機(jī)數(shù)，在由y h(x)I(x), I(x) 1 ,x (a,b),得到隨機(jī)變量Y的隨機(jī)數(shù)，求出樣本均 0 ,x (a,b)值，以此估計(jì)積分值。積分S g (x, y) dxdy的求法與上述方法類(lèi)似，在此不贅述。A概率密度函數(shù)的選取：0(a,b),為使方法普一重積分，由于要求f(x)的支持域x| f(x)x2遍適用，考慮到標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)f(x) 1 e萬(wàn)支持域?yàn)镽 ,故

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

5、unction 函數(shù)內(nèi)容。編程完整解決用蒙特卡洛方法估計(jì)一重、二重積分值問(wèn)題。程序代碼及運(yùn)行結(jié)果：第一類(lèi)一重積分程序代碼：%構(gòu)造示性函數(shù)%function I=I1（x,a,b）if x>=a&&x<=bI=1;elseI=0;end%保存為I1.m%第一類(lèi)一重積分, 程序主體：%保存為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; % 單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y1*t'/10; % 十次均值D=abs(EY-Real); % 絕對(duì)誤差RD=D/Real; % 絕對(duì)誤差d=0;for i=1:10d=d+(Y1(i)-Real2;endd=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10 Am時(shí)的樣本均值D1(m-V(1)+1)=D;% 絕對(duì)誤差RD1(m-V(1)+1)=RD;% 絕對(duì)誤差MSE1(m-V(1)+1)=d; % 方差endReal,EY1,D1,RD1,MSE1outf11=EY1;

8、D1;RD1;MSE1; % 存放樣本數(shù)字特征%保存為f11.m歡迎下載運(yùn)行結(jié)果:%估計(jì)積分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%估計(jì)積分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歡迎下載第一類(lèi)二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)%funct

11、ion I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存為I2.m%第一類(lèi)二重積分程序主體%保存為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; % 單次實(shí)驗(yàn)樣本均值endt=ones(1,10);EY=Y2*t'/10; % 十次均值D=abs(EY-Real); % 絕對(duì)誤差RD

13、=D/Real; % 絕對(duì)誤差d=0;for i=1:10d=d+(Y2(i)-Real2;endd=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為10 Am時(shí)的樣本均值D2(m-V(1)+1)=D;% 絕對(duì)誤差RD2(m-V(1)+1)=RD;% 絕對(duì)誤差MSE2(m-V(1)+1)=d; % 方差歡迎下載endReal,EY2,D2,RD2,MSE2outf12=EY2;D2;RD2;MSE2; % 存放樣本數(shù)字特征%保存為f12.m運(yùn)行結(jié)果：%估計(jì)積分 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歡迎下載第二類(lèi)一重積分程序代碼：%構(gòu)造示性函數(shù)%function I=I1(x,a,b)if x>=a&&x<=bI=1;elseI=0;end%保存為I1.m%第二類(lèi)一重積分程序主體%程序保存為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; % 單次實(shí)驗(yà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時(shí)的樣本均值MSE1(m-V(1)+1)=d;% 方差endEY1,MSE1outf21=EY1;MSE1; % 存放樣本數(shù)字特征%程序保存為f21.m運(yùn)行結(jié)果：12%估計(jì)積分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第二類(lèi)二重積分程序代碼：%構(gòu)造示性函數(shù)

19、，求不同區(qū)域上積分只需更改示性函數(shù)%function I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存為I2.m%第二類(lèi)二重積分函數(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; % 單次實(shí)驗(yà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時(shí)的樣本均值MSE2(m-V(1)+1)=d; % 方差 endEY2,MSE2outf22=EY2;MSE2; % 存放樣本數(shù)字特征%第二類(lèi)二重積分，程序保存為f22.m運(yùn)行結(jié)果:%估計(jì)積分 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實(shí)驗(yàn)結(jié)果整理: 第一類(lèi)一重積分:估計(jì)積分2xsin xdx0積分真值:積分估計(jì)值:1.0018樣本容量:10100100010000樣本均值:1.26351.00881.00661.01091000001.0018絕對(duì)誤差：0.26350.00880.00660.01090.0018相對(duì)誤差：0.2635均方誤差：0.64392估計(jì)積分e-x dx0積分真值：0.8862樣本容

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

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

25、992.85662.8269樣本方差：1.32670.09000.00600.0014實(shí)驗(yàn)結(jié)果分析：22從第一類(lèi)積分看，以估計(jì)積分xsin xdx為例:0積分真值：1積分估計(jì)值：1.0018樣本容量：10樣本均值：1.2635絕對(duì)誤差：0.2635相對(duì)誤差：0.2635均方誤差：0.643910010001.00881.00660.00880.00660.00880.00660.02050.0028100001000001.01091.00180.01090.00180.01090.00180.00060.0001隨著樣本容量的增大，樣本均值有接近積分真值的趨勢(shì)，絕對(duì)誤 差、相對(duì)誤差、均方誤差呈減小趨勢(shì)；隨著樣本容量的增大