用蒙特卡洛方法估計(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)容：1用蒙特卡洛方法估計(jì)積分Jxsinxdx, Je-xdx和JJ e*“dxdy的00x2：；y2<1值，并將估計(jì)值與真值進(jìn)行比較。 12用家特卡洛方法估計(jì)積分(e dx和1 , 4 4 dxdy的值,0x2 72 j . 1 x y并對(duì)誤差進(jìn)行估計(jì)。要求：(1)針對(duì)要估計(jì)的積分選擇適當(dāng)?shù)母怕史植荚O(shè)計(jì)蒙特卡洛方法；(2 )利用計(jì)算機(jī)產(chǎn)生所選分布的隨機(jī)數(shù)以估計(jì)積分值；(3)進(jìn)行重

2、復(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é)果的精度。目的：(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)的積分值。具體操作如下：般地，積分S=Jg(x)dx改寫(xiě)成S=Jbf (x)dx = h(x) f

3、(x)dx 的a形式，(其中為f(x)一隨機(jī)變量X的概率密度函數(shù)，且f(x)的支持域x | f (x) >0 n(a,b) ) , h(x)=g(x)f(x);令Y=h(X),貝用分S二E(Y);利用matlab軟件，編程產(chǎn)生隨機(jī)變量X的隨機(jī)數(shù),在由 y=h(x)I(x), I(x)=1,z(a，b)得到隨機(jī)變量 y的隨機(jī)數(shù), 0 ,x“a,b)''求生樣本均值，以此估計(jì)積分值。積分S = g(x,y)dxdy的求法與上述方法類(lèi)似，在此不贅述。 A概率密度函數(shù)的選?。簒222o一重積分，由于要求f(x)的支持域x|f(x)>03(a,b),為使 方法普遍適用，考慮到

4、標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)x2f(x)=忑”2支持域?yàn)镽,故選用f(x)=72e1類(lèi)似的，二重積分選用2 y2二，支持域?yàn)镽2o估計(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)造，求不

5、同區(qū)域二重積分，只需改變function函數(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)一重積分，程序主體：%保存為fll.mfunction outf11=f11()g1=input('輸入一元被積函數(shù)如 x.*sin(x):','s')%輸入被積函數(shù)g1=inline(g1);a=input('輸入積分下界a:');%輸入積分上

6、下限b=input('輸入積分上界 b:');Real=input('積分真值：');輸入積分真值fprintf('輸入樣本容量10AV1-10AV2:r')V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');for m=V(1):V(2)%樣本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);for i=1:nt1(i)=I1(x(i),a,b);%示性及求和向量end y1=g1(x)*(pi*2)A0.5).*ex

7、p(x.A2/2);Y1(j)=y1*t1'/n; % 單次實(shí)驗(yàn)樣本均值 end t=ones(1,10);EY=Y1*t'/10;%十次均值D=abs(EY-Real); % 絕對(duì)誤差 RD=D/Real;%絕對(duì)誤差d=0;foh=1:10d=d+(Y1(i)-Real)A2;%樣本容量為 10Am%絕對(duì)誤差%絕對(duì)誤差時(shí)的end d=d/(10-1); EY1(m-V(1)+1)=EY;樣本均值D1(m-V(1)+1)=D;RD1(m-V(1)+1)=RD;MSE1(m-V(1)+1)=d;% 方差end Real,EY1,D1,RD1,MSE1 outf11=EY1;D1;

8、RD1;MSE1; % 存放樣本數(shù)字特征 %保存為fll.m運(yùn)行結(jié)果:%估計(jì)積分fxsinxdx,積分真值為1m=f11輸入一元被積函數(shù)如x.*sin(x):x.*sin(x) g1 = x.*sin(x)輸入積分下界a:0輸入積分上界b:pi/2積分真值:1輸入樣本容量10AV1-10AV2:V1:1V2:510100 n =1000n =10000n =100000Real =11.0018EY1 =1.26351.00881.00661.0109D1 =0.26350.00880.00660.01090.00180.00180.00011.00180.00180.00180.00010.

9、64390.02050.00280.0006m=1.26351.00881.00661.01090.26350.00880.00660.01090.26350.00880.00660.01090.64390.02050.00280.0006RD1 =0.26350.00880.00660.0109MSE1 =%估計(jì)積分 ；e-x2dx真值為0.8862 0M=f11輸入一元被積函數(shù)如x.*sin(x):exp(-x.A2)g1 =exp(-xd2)輸入積分下界a:0輸入積分上界b:+inf積分真值：piA0.5/2%0.8862輸入樣本容量 10AV1-10AV2:V1:1V2:4101001

10、00010000Real =0.8862EY1 =0.9333D1 =0.0470RD1 =0.0531MSE1 =0.1927M =0.93330.04700.05310.19270.90770.88730.88710.02150.00100.00090.02430.00120.00100.01120.00160.00000.90770.88730.88710.02150.00100.00090.02430.00120.00100.01120.00160.0000第一類(lèi)二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)function I=I2(x,y)if xA2+yA2&

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

12、%輸入樣本容量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(xA2+yA2)/2);Y2(j)=y2*t2'/n; %單次實(shí)驗(yàn)樣本均值 endt=ones(1,10);EY=Y2*t'/10;D=abs(EY-Real); % RD=D/Real;d=0;foh=1:10d=d+(Y2(i)-Real)A

13、2; end d=d/(10-1);EY2(m-V(1)+1)=EY;樣本均值D2(m-V(1)+1)=D;RD2(m-V(1)+1)=RD;MSE2(m-V(1)+1)=d;%十次均值絕對(duì)誤差%絕對(duì)誤差%樣本容量為 10Am 時(shí)的%絕對(duì)誤差%絕對(duì)誤差%方差end存放樣本數(shù)字特征Real,EY2,D2,RD2,MSE2 outf12=EY2;D2;RD2;MSE2; % %保存為f12.m運(yùn)行結(jié)果:%估計(jì)積分 ex*dxdy,真值為 pi*(exp(1)-1)%5.3981 22x y也m=f12輸入二元被積函數(shù)如 exp(x.A2+y.A2):exp(x.A2+y.A2)g2 =exp(x.

14、A2+yd2)積分真值：pi*(exp(1)-1)%5.3981輸入樣本容量 10AV1*10AV1-10AV2*10AV2：V1:1V2:4 n = 101001000n =10000Real =5.3981EY2 =4.77025.12505.43175.4041D2 =0.62790.27320.03350.0060RD2 =0.11630.05060.00620.0011MSE2 =3.89650.55640.02470.00174.77025.12505.43175.40410.62790.27320.03350.00600.11630.05060.00620.00113.89650

15、.55640.02470.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('輸入積分下界a:');%輸入積分上下限 b=input('輸入積分上界 b:');fprintf('

16、;輸入樣本容量10AV1-10AV2:r')V=zeros(1,2);V(1)=input('V1:');% 輸入樣本容量V(2)=input('V2:');for m=V(1):V(2)%樣本容量 10Am1-10Am2n=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);Y1(j)=y1*t1'/n; % 單次實(shí)驗(yàn)樣本均值 end t=ones(1,10);EY=Y1*t/10；%十次均值d

17、=0; fo門(mén)=1:10d=d+(Y1(i)-EY)A2;end d=d/(10-1);EY1(m-V(1)+1)=EY;%樣本容量為10Am 時(shí)的樣本均值MSE1(m-V(1)+1)=d;% 方差end EY1,MSE1outf21=EY1;MSE1; % 存放樣本數(shù)字特征 %程序保存為f21.m運(yùn)行結(jié)果：1 %估計(jì)積分exdx0m=f21輸入一元被積函數(shù)如 exp(x.A2):exp(x.A2)g1 =exp(xd2)輸入積分下界a:0 輸入積分上界b:1 輸入樣本容量 10AV1-10AV2:V1:1V2:4 n =10100n = 1000n = 10000EY1 =1.45900.0

18、0082.07821.65831.5029MSE1 =0.43150.08890.0057m =2.07821.65831.50291.45900.43150.08890.00570.000812%用matlab指令求積分fe dx0f=inline('exp(x.A2)')f =Inline function:f(x) = exp(x-2)>> S=quadl(f,0,1)S =1.4627第二類(lèi)二重積分程序代碼：%構(gòu)造示性函數(shù)，求不同區(qū)域上積分只需更改示性函數(shù)function I=I2(x,y)if xA2+yA2<=1I=1;elseI=0;end%保存

19、為I2.m%第二類(lèi)二重積分函數(shù)主體%,程序保存為f22.mfunction outf22=f22()g2=input(' 輸 入 二 元 被 積 函 數(shù) 如1./(1+xA4+yA4)A0.5：','s')%輸入被積函數(shù)g2=inline(g2,'x','y');fprintf('輸入樣本容量10AVl*l0AVi-i0AV2*i0AV2：r')V=zeros(1,2);V(1)=input('V1:');%輸入樣本容量V(2)=input('V2:');for m=V(1):V(2

20、)% 樣本容量 10Am1-10Am2n=10Amfor j=1:10x=randn(1,n);y=randn(1,n);for i=1:nt2(i)=I2(x(i),y(i);% 示性及求和向量 end y2=g2(x,y)*(2*pi).*exp(xA2+yA2)/2);Y2(j)=y2*t2'/n; % 單次實(shí)驗(yàn)樣本均值 end t=ones(1,10);EY=Y2*t'/10;%十次均值d=0; fo門(mén)=1:10d=d+(Y2(i)-EY)A2;end d=d/(10-1);EY2(m-V(1)+1)=EY;%樣本容量為 10Am 時(shí)的樣本均值 MSE2(m-V(1)+

21、1)=d;% 方差end EY2,MSE2 outf22=EY2;MSE2; % 存放樣本數(shù)字特征 %第二類(lèi)二重積分，程序保存為f22.m運(yùn)行結(jié)果：%估計(jì)積分一1-dxdyx>yZ:1 x4 -y4m=f22輸 入 二 元 被 積 函 數(shù) 如1./(1+x.A4+yA4)A0.5：1./(1+xA4+yA4)A0.5g2 =1./(1+xA4+yA4)A0.5輸入樣本容量 10AV1*10AV1-10AV2*10AV2:V1:1V2:4 n =10n =100 n =100010000EY2 =3.07592.96992.85662.8269MSE2 =1.32670.09000.006

22、00.0014m =3.07592.96992.85661.32670.09000.00602.82690.0014實(shí)驗(yàn)結(jié)果整理:第一類(lèi)一重積分:估計(jì)積分積分真值： 樣本容量： 樣本均值 1.0018 絕對(duì)誤差 0.0018 相對(duì)誤差 0.0018均方誤差 0.0001n2 xsin xdx01積分估計(jì)值:10100:1.26350.26350.26350.64391000100001000001.00881.00661.01090.00880.00660.01090.00880.00660.01090.02050.00280.00061.0018-be-x2 .e dx樣本容量：101001

23、00010000樣本均值：0.93330.90770.88730.8871絕對(duì)誤差：0.04700.02150.00100.0009相對(duì)誤差：0.05310.02430.00120.0010均方誤差：0.19270.01120.00160.0000估計(jì)積分積分真值積分估計(jì)值:0.88620.88710第一類(lèi)二重積分:22估計(jì)積分ex ydxdy22 dX招色積分真值：5.3981樣本容量：10100積分估計(jì)值： 5.4041100010000樣本均值:絕對(duì)誤差:相對(duì)誤差:均方誤差:4.77020.62790.11633.89655.12500.27320.05060.55645.43170.0

24、3350.00620.02475.40410.00600.00110.0017第二類(lèi)一重積分：12估計(jì)積分exdx0積分估計(jì)值：1.4590樣本容量：10100100010000樣本均值：2.07821.65831.50291.4590樣本方差：0.43150.08890.00570.0008用matlab指令求得積分結(jié)果1.4627第二類(lèi)二重積分：估計(jì)積分14 dxdyx2y2j-.1 X4 y4積分估計(jì)值：2.8269樣本容量：101001000100002.85662.82690.00600.00141000100001.00881.00660.00880.00660.00880.00660.02050.00281000001.01090.01090.01090.0006樣本均值：3.07592.9699樣本方差：1.32670.0900實(shí)驗(yàn)結(jié)果分析:工2從第一類(lèi)積分看，以估計(jì)積分 fxsinxdx為例:0積分真值：1積分估計(jì)值：1.0018樣本容量：10100樣本均值：1.26351.0018絕對(duì)誤差：0.26350.0018相對(duì)誤差：0.26350.0018均方誤差：0.64390.0001隨著樣本容量的增大，樣本均值有接近積分真值的趨勢(shì), 絕對(duì)誤差、相對(duì)誤差、均方誤差呈減小趨勢(shì)；隨著樣本容量 的增大，樣本均值有接近積分真值的趨勢(shì)，說(shuō)明估計(jì)具有