matlab蒙特卡洛法估計(jì)積分值

課程：概率論與數(shù)理統(tǒng)計(jì)實(shí)驗(yàn)日期：報(bào)告日期：專業(yè)班級：姓名：學(xué)號：實(shí)驗(yàn)內(nèi)容：用蒙特卡洛方法估計(jì)積分值要求：（1）針對要估計(jì)的積分選擇適當(dāng)?shù)母怕史植荚O(shè)計(jì)蒙特卡洛方法；（2）利用計(jì)算機(jī)產(chǎn)生所選分布的隨機(jī)數(shù)以估計(jì)積分值；（3）進(jìn)行重復(fù)試驗(yàn)，通過計(jì)算樣本均值以評價(jià)估計(jì)的無偏性；通過計(jì)算均方誤差（針對第1類題）或樣本方差（針對第2類題）以評價(jià)估計(jì)結(jié)果的精度。目的：（1）能通過MATLAB或其他數(shù)學(xué)軟件了解隨機(jī)變量的概率密度、分布函數(shù)及其期望、方差、協(xié)方差等；（2）熟練使用MATLAB對樣本進(jìn)行基本統(tǒng)計(jì)，從而獲取數(shù)據(jù)的基本信息；（3）能用MATLAB熟練進(jìn)行樣本的一元回歸分析。1用蒙特卡洛方法估計(jì)積分，和的值，并將估計(jì)值與真值進(jìn)行比較。1)仍是用均勻分布來估計(jì)此積分的大小，g(x)=xsinx,=1/().x>0.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;r=0;n=50;h(1:10)=0;forj=1:10fori=1:na=unifrnd(0,pi/2,n,1);x=sort(a);y=pi/2*mean(x.*sin(x));s=s+y;endb=s./n;fprintf('b=%.4f\n',b);h(j)=b;s=0;m=m+b;endp=m./10z=1forj=1:10r=(h(j)-z).^2;f=f+r;endf=f./10;fprintf('f=%.6f\n',f)221，表明估計(jì)結(jié)果與理論值非常接近。2)I==1/2*g(x)=e為標(biāo)準(zhǔn)正態(tài)分布的概率密度.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;n=50;r=0;h(1:10)=0;forj=1:10fori=1:na=normrnd(0,1,1,n);x=sort(a);z=(sqrt(2.*pi)).*exp(-x(i).^2./2);s=s+z;endb=(s./n)./2;fprintf('b=%.4f\n',b);h(j)=b;s=0;m=m+b;endp=m./10z=sqrt(pi)./2forj=1:10r=(h(j)-z).^2;f=f+r;endf=f./10;fprintf('f=%.6f\n',f)結(jié)果如下：322，估計(jì)結(jié)果與真實(shí)值非常接近。3)m=10000;sum=0;n=50;D=0;X=unifrnd(-1,1,n,m);Y=unifrnd(-1,1,n,m);fori=1:na=0;forj=1:mif(X(i,j)^2+Y(i,j)^2<=1)Z(i,j)=exp(X(i,j)^2+Y(i,j)^2);a=a+Z(i,j);endendS(i)=a/m;sum=sum+S(i);endI=sum/n*4fori=1:nD=D+(S(i)*4-pi*(exp(1)-1))^2;endd=D/n2用蒙特卡洛方法估計(jì)積分和的值，并對誤差進(jìn)行估計(jì)。1)此積分采用的是均勻分布。g(x)=,=1.x>0.分別取10個(gè)估計(jì)值h(j)，求得估計(jì)值的均值p，對照積分的真實(shí)值求得估計(jì)均方誤差f。Matlab程序代碼如下：s=0;m=0;f=0;r=0;n=50;h(1:10)=0;forj=1:10fori=1:na=unifrnd(0,1,n,1);x=sort(a);y=exp(x(i).^2);s=s+y;endb=s./n;fprintf('b=%.4f\n',b);h(j)=b;s=0;m=m+b;endp=m./10forj=1:10r=(h(j)-p).^2;f=f+r;endf=f./9;fprintf('f=%.6f\n',f)結(jié)果如下：22，以平均值作為真實(shí)值，均方誤差也比較小。2)n=1000;m=100;sum=0;S=0;I=0;x=unifrnd(-2,2,m,n);y=unifrnd(-2,2,m,n);forj=1:ms=0;fori=1:nifx(j,i)^2+y(j,i)^2<=4s=s+16/sqrt(1+x(j,i)^4+y(j,i)^2);