matlab蒙特卡洛法估計(jì)積分值[6頁(yè)]

1、西安交通大學(xué)實(shí)驗(yàn)報(bào)告 課程：概率論與數(shù)理統(tǒng)計(jì)實(shí)驗(yàn)日期：報(bào)告日期：專業(yè)班級(jí)： 姓 名： 學(xué) 號(hào)：實(shí)驗(yàn)內(nèi)容：用蒙特卡洛方法估計(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ì)算樣本均值以評(píng)價(jià)估計(jì)的無(wú)偏性；通過(guò)計(jì)算均方誤差（針對(duì)第1類題）或樣本方差（針對(duì)第2類題）以評(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)行樣本

2、的一元回歸分析。1用蒙特卡洛方法估計(jì)積分 ，和的值，并將估計(jì)值與真值進(jìn)行比較。1)仍是用均勻分布來(lái)估計(jì)此積分的大小，g(x)=xsinx,=1/().x0.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對(duì)照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;r=0;n=50;h(1:10)=0;for j=1:10 for i=1:na=unifrnd(0,pi/2,n,1); x=sort(a); y=pi/2*mean(x.*sin(x); s=s+y; end b=s./n; fprintf(b=%.4fn,b); h(j)=b; s=0; m=m+b;e

3、ndp=m./10 z=1for j=1:10 r=(h(j)-z).2; f=f+r;endf=f./10;fprintf(f=%.6fn,f)結(jié)果顯示f=0.000221，表明估計(jì)結(jié)果與理論值非常接近。2) I=1/2* g(x)=e 為標(biāo)準(zhǔn)正態(tài)分布的概率密度.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對(duì)照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;n=50;r=0;h(1:10)=0;for j=1:10 for i=1:n a=normrnd(0,1,1,n); x=sort(a); z=(sqrt(2.*pi).*exp(-x(i).2./

4、2); s=s+z; end b=(s./n)./2; fprintf(b=%.4fn,b); h(j)=b; s=0; m=m+b;endp=m./10z=sqrt(pi)./2for j=1:10 r=(h(j)-z).2; f=f+r;endf=f./10;fprintf(f=%.6fn,f)結(jié)果如下：結(jié)果顯示估計(jì)結(jié)果與真實(shí)值的方差為f=0.00322，估計(jì)結(jié)果與真實(shí)值非常接近。3) m=10000;sum=0;n=50;D=0;X=unifrnd(-1,1,n,m);Y=unifrnd(-1,1,n,m);for i=1:n a=0; for j=1:m if(X(i,j)2+Y(i,

5、j)20.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對(duì)照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;r=0;n=50;h(1:10)=0;for j=1:10 for i=1:n a=unifrnd(0,1,n,1); x=sort(a); y=exp(x(i).2); s=s+y; end b=s./n; fprintf(b=%.4fn,b); h(j)=b; s=0; m=m+b;endp=m./10for j=1:10 r=(h(j)-p).2; f=f+r;endf=f./9;fprintf(f=%.6fn,f)結(jié)果如下：結(jié)果顯示，誤差為0.000322，以平均值作為真實(shí)值，均方誤差也比較小。2) n=1000;m=100;sum=0;S=0;I=0;x=unifrnd(-2,2,m,n);y=unifrnd(-2,2,m,n);for j=1:m s=0; for i=1:n if x(j,i)2+y(j,i)2=4 s=s+16/sqrt(1+x(j,i)4+y(j,i)2);