概率實(shí)驗(yàn)報(bào)告-蒙特卡洛積分

本科實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱(chēng):《概率與統(tǒng)計(jì)》隨機(jī)模擬實(shí)驗(yàn)隨機(jī)模擬實(shí)驗(yàn)實(shí)驗(yàn)設(shè)隨機(jī)變量X的分布律為P{X=i}=2-i,i=1,2,3 試產(chǎn)生該分部的隨機(jī)數(shù)1000個(gè)，并作出頻率直方圖。一、實(shí)驗(yàn)原理采用直接抽樣法：定理：設(shè)U是服從[0,1]上的均勻分布的隨機(jī)變量，則隨機(jī)變量y=F-1(U)與X有相同的分布函數(shù)Y=F-1(U)(為F(x)的逆函數(shù))，即Y=F-1(U)的分部函數(shù)為F(x).二、題目分析易得題中X的分布函數(shù)為F(x)=1--,i<x<i+1,i=0,1,2,3,……2i若用ceil表示對(duì)小數(shù)向正無(wú)窮方向取整，則F(x)的反函數(shù)為F-1(〃)=ceilI叱")I

[產(chǎn)生服從[0,1]上的均勻分布的隨機(jī)變量a，則m=F-1(a)則為題中需要產(chǎn)生的隨機(jī)數(shù)。三、MATLAB實(shí)現(xiàn)f=[];i=1;whilei<=1000a=unifrnd(0,1); %產(chǎn)生隨機(jī)數(shù)a，服從【0,1】上的均勻分布m=log(1-a)/log(1/2);b=ceil(m);%對(duì)m向正無(wú)窮取整f=[f,b];i=i+1;enddisplay(f);[n,xout]=hist(f);bar(xout,n/1000,1)產(chǎn)生的隨機(jī)數(shù)(取1000個(gè)中的20個(gè))如下：

i12345678910f(i)2233124133i11121314151617181920f(i)1221111112頻率分布直方圖0.60.70.60.60.40.60.70.60.60.4□.3.0.20.10 .2. 46 8 1012實(shí)驗(yàn)二設(shè)隨機(jī)變量X的密度函數(shù)為，/、f4xe-2x,x>0,

f(x)=<0,x<0試產(chǎn)生該分布的隨機(jī)數(shù)1000個(gè)，并作出頻率直方圖一、實(shí)驗(yàn)原理取舍抽樣方法，當(dāng)分布函數(shù)的逆函數(shù)難以求出時(shí)，可采用此方法。取舍抽樣算法的流程為：(1)選取一個(gè)參考分布，其選取原則，一是該分布的隨機(jī)樣本容易產(chǎn)生；二是存在常數(shù)C，使得f(x)<Cg(x)。

產(chǎn)生參考分布g(x)的隨機(jī)樣本x；0獨(dú)立產(chǎn)生[0,1]產(chǎn)生參考分布g(x)的隨機(jī)樣本x；0獨(dú)立產(chǎn)生[0,1]上的均勻分布隨機(jī)數(shù)u；0若uCg(x)<f(x)00二、題目分析選取參考分布則保留x0,作為所需的隨機(jī)樣本；否則舍棄。Ie-x,x>0,g(x)=LcI0,x<0則有f(x)<2g(x)，密度函數(shù)為g(x)隨機(jī)變量的產(chǎn)生采用題1中的方法即可三、MATLAB實(shí)現(xiàn)f=[];i=1;whilei<=1000a=rand(1); %產(chǎn)生[0,1]上的隨機(jī)分布的數(shù)ab=-log(1-a); %利用直接抽樣法產(chǎn)生密度函數(shù)為g(x)的隨機(jī)數(shù)bc=rand(1); %產(chǎn)生[0,1]上的隨機(jī)分布的數(shù)aifc*2*(1-a)<=4*b*exp(-2*b)%判斷c*2*g(b)是否小于f(b),若滿足，則將b作為所需的隨機(jī)樣本保留f=[f,b];i=i+1;endenddisplay(f);[n,xout]=hist(f);bar(xout,n/sum(n),1)產(chǎn)生的部分隨機(jī)數(shù)(取1000個(gè)中的20個(gè))如下：i12345678910f(i)0.76420.40641.55131.08330.07680.41271.34141.44290.33720.7056i11121314151617181920f(i)0.56730.17380.97900.93140.21520.19140.24050.26000.88501.2865頻率分布直方圖

實(shí)驗(yàn)三兩汽車(chē)司機(jī)每天早上8:00—8:01駕車(chē)到達(dá)某十字路口。若到達(dá)的時(shí)間間隔小于15秒，他們必須停車(chē)相互避讓。問(wèn)他們過(guò)此十字路口需要停車(chē)避讓的概率？如果是三輛、四輛汽車(chē)呢？此概率與汽車(chē)數(shù)的關(guān)系怎樣？有多少汽車(chē)時(shí)有必要安裝紅綠燈呢？分析求解：首先我們考慮樣本空間。如果把n臺(tái)車(chē)的到達(dá)時(shí)間看為一個(gè)序列｛Xn｝?？紤]兩臺(tái)車(chē)到達(dá)時(shí)間分別為了,%+At在同一個(gè)時(shí)間點(diǎn)內(nèi)到達(dá)的概率limJtAtdx=limtAt=0At-00 At-0則可認(rèn)為序列｛X｝幾乎兩兩不相等。n若考慮2臺(tái)車(chē)在60秒內(nèi)，間隔時(shí)間超過(guò)15秒算不擁堵，設(shè)A(2,60,15)為不會(huì)擁堵,則：J45dxJ60dx9P(A(2,60,15))=升—1產(chǎn)2=-J60dxJ60dx 16x112x1則由此可知2臺(tái)車(chē)在60秒內(nèi)，間隔時(shí)間超過(guò)15秒算不擁堵，不擁堵的概率為2。16若考慮三臺(tái)車(chē)，四臺(tái)車(chē)只需將2臺(tái)車(chē)的積分公式進(jìn)行拓展便可求出答案。現(xiàn)在考慮n臺(tái)車(chē)在t秒內(nèi)，間隔時(shí)間超過(guò)s秒算不擁堵，設(shè)A(n,t,s)為事件“不會(huì)擁堵“,則：上述積分式Jt—(n—1)sdxJt—(n—2)sd上述積分式Jt—(n—1)sdxJt—(n—2)sdx...Jt—s20 s . x2+stdx...Jtdx1 2 n—10x x1 —20,elsedxn—1dx1+s n,t>(n—1)sx—1dxnJtdxJtdx...JtdxJt1 2 n—10x x x1 n—2 n—1dxnJt—(n—1)sdxJt—(n—2)sdx...Jt—s dx Jt dx1x1+sxn—2+sn—1x +sn—1可用數(shù)學(xué)歸納法求得：JtdxJt0 1x1dx...Jt2dxJtn—1dxntnJt—(n—1)s Jt—(n—2)s Jt—sJt (t—(n—1)s)nx C/v,^xJ C/v,^x...J C/v,^x J C/v,^xTOC\o"1-5"\h\z1 2 n—1 n n!0 x+s x+s x+s n1 n—2 n—1則有：t>(n—1)selse(1(nt>(n—1)selseP(A(n,t,s))=<1—t—/°，為了驗(yàn)證上述公式的正確性，我們用計(jì)算機(jī)進(jìn)行了模擬(模擬107組)。模擬的思路是產(chǎn)生服從給定時(shí)間范圍內(nèi)的均勻分布的大量隨機(jī)數(shù)，通過(guò)判斷條件來(lái)計(jì)數(shù)需要停車(chē)避讓的點(diǎn)，求得其出現(xiàn)頻率，由大數(shù)定理可得需要停車(chē)避讓的概率。隨機(jī)模擬C++代碼如下：#include<cmath>#include<dlgoritlmi>using-namespacestd;iLitTLfk;doubles,t;doublea二口口口口口二，；doublegetrand(doiibley)return(irand()-EtAND_MAX-iand(J)/(double)(EtAHD_MAX-RAND_MAX))*yintmain()Brand(tiir.e(0));cout<<,rn=,r;cin>>n;aout<<,rt=,r;cin>>t;cont<<"s=*';cin>>s;k=_ooooaoo；intent(0);for(intL=L;i<=k;--i)foriintj=L;j<=n;——j)a'j'=getrand.(t);與口r七(己一二ra-L-n);boolf口;for(intj=L;j<n.;——j)i￡(a；j-1；-a；j；<5)■-：￡=Q;I>reak;；if(fi-一ent：printf「pf二i=W. pow((L.0-(ti-L)"b"L.0/t>,n));printfi,rpl=^.WXX'n",cnt"L.(?/k);return0;設(shè)p0為上述公式計(jì)算值，p1為計(jì)算機(jī)模擬值。ntsP0P1260150.56250.562617360150.1250.1249852460150.003906250.0039277101000100.389416120.3893441301000100.000034490.00003780501000050.289310160.2893139由上表觀察到P0與p1符合程度較好?？苫菊J(rèn)為推導(dǎo)公式正確。兩汽車(chē)司機(jī)需要停車(chē)避讓的概率為0.5625,三輛汽車(chē)時(shí)概率為0.125，四輛汽車(chē)時(shí)0.0039。關(guān)于是否安裝紅綠燈，可根據(jù)

1一(n^s、Vt70,else若要使不堵車(chē)的概率大于a,則有：1-（n^Vt7故當(dāng)V1故當(dāng)V1-7）n<a時(shí)，需安裝紅綠燈進(jìn)行調(diào)控。實(shí)驗(yàn)四試采用兩種方法，模擬計(jì)算下定積分卜sx0.9e-xdx0并分析其誤差。要計(jì)算的定積分是在無(wú)窮區(qū)間上。考慮到被積函數(shù)收斂速度很快（在x為40時(shí)已經(jīng)近似收斂到0），我們用自適應(yīng)辛普森公式（目前精度最高的一種差值積分算法）計(jì)算J2000x0.9e-xdx,求得其近似實(shí)際值作為隨機(jī)模擬誤差分析的依0據(jù),C++程序如下：#include<stdio.h>#include<math.h>#include<stdlib.h>#include<algorithm>#include<iostream>usingnamespacestd;doublex[2100000];doublef(doublex){returnexp(log(x)*0.9-x);}intmain(){intn=1000000;for(inti=1;i<=n;++i)x[i]=i*2000.0/n;doubleans=0;for(inti=1;i<=n/2;++i)ans+=2000.0/(3*n)*(f(x[2*i-2])+4*f(x[2*i-1])+f(x[2*i]));printf("%.8lf\n",ans);return0;}執(zhí)行后得到積分近似值0.961766.下面采用蒙特卡洛法的思想進(jìn)行隨機(jī)模擬計(jì)算定積分：(一)平均值法：原理：若％~U(a,b)，則有E(g(%))=fbg(x)f(x)dx=一--fbg(x)dx，因此要a b-aa計(jì)算fbg(x)dx，只需求隨機(jī)變量y=g(x)的平均值，再乘以區(qū)間長(zhǎng)度即可。a求解：Matlab程序如下：x=unifrnd(a,b,[1n]);fori=1:nm(i)=fun(x(i));endmean(m);p=mean(m)*(b-a).當(dāng)a=100,b=200,n=106時(shí)p=0，故可認(rèn)為在0到100上取x計(jì)算得到的積分值即為0到+8上的積分值；.當(dāng)a=0,b=100,n=106時(shí)重復(fù)五次實(shí)驗(yàn)得序號(hào)12345P0.96740.96670.96820.96970.9643平均值P11-0.96726,|0.96726一0.961766|百分誤差8 x100%“0.571241%.11 0.9617663.當(dāng)a=0,b=100,n=107時(shí)重復(fù)五次實(shí)驗(yàn)得序號(hào)12345P0.95820.96410.96380.96190.9621平均值P12-0.96202,|0.96202一0.961766|百分誤差8 x100%h0.026410%.12 0.961766分析：可以看出，所投點(diǎn)越多，所求值越接近實(shí)際值，誤差越小(二)隨機(jī)投點(diǎn)法原理：設(shè)r.v.(XY服從矩形他<%<b,c<y<d｝上的均勻分布，則％~U(a,b),y~U(c,d),且X,Y相互獨(dú)立。記事件A=｛Y<f(x)｝,則P(A)=P(Y<f(x))=Jbff(x)dydx-Jb(f(%)-c)dx-Jbf(%)dx一c(b一a)-uac a a右側(cè)積分值u即為A出現(xiàn)的頻率。由伯努利大數(shù)定律，用重復(fù)試驗(yàn)中A出現(xiàn)的頻率作為u的估計(jì)值。將(X,Y)看成矩形內(nèi)的隨機(jī)投點(diǎn)，用隨機(jī)點(diǎn)落在區(qū)域Y<f(%)中的頻率作為定積分的近似值。求解：(1)為了確定y的取值范圍使得投點(diǎn)范圍更小，先求f(%)-%0.9e-%在0到100上的最大值：symsx;y=@(x)(-xA0.9*exp(-x))>>c=fminbnd(y,0,100)c=0.9000>>0.9A0.9*exp(-0.9)ans=0.3698得求0到100上的最大值為0.3698(2)在｛100<%<OQ0 <03嵌｝ 上投107個(gè)點(diǎn)，計(jì)算被f(%)包含的點(diǎn)出現(xiàn)的頻率N=10A7;m=0;x=unifrnd(100,200,[1,N]);

y=unifrnd(0,0.3698,[1,N]);m=y<((x.A0.9).*exp(-x));n=sum(m);p=n/N*100*0.3698得p=0故可認(rèn)為在0到100上取x計(jì)算得到的積分值即為0到上的積分值.(3)在｛0<%<100,0<y<0.3698｝上投