復(fù)合泊松過程的實(shí)現(xiàn)

-.z.電子信息與通信工程學(xué)院實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱非其次泊松過程課程名稱隨機(jī)信號(hào)分析顧康**U201413323日期6.13地點(diǎn)南一樓成績(jī)教師董燕題目ConsiderthenonhomogeneousPoissonprocesswithitsintensityfunctionspectifiedinE*ample2.3.6.(a)WriteaMATLABprogramtogenerate(stimulate)thefirsteightyarrivaltimes.(b)Givent=8(hours),writeaMatlabprogramtogenerateN(8)andthenthearrivaltimesintheinterval(0,8],drawtherespectivehistogramsshowinghour5yarrivalcounts.(a)由定理設(shè)λ〔ｔ〕≤λ，其中λ為一常數(shù)，而ｓ１，ｓ２，…，ｓｎ，…為參數(shù)λ的齊次泊松過程的事件發(fā)生的時(shí)刻，對(duì)每個(gè)ｓｉ，以概率λ〔ｓｉ〕／λ進(jìn)展保存，以概率１－λ〔ｓｉ〕／λ舍棄，由此得到的序列ｓ〔１〕，ｓ〔２〕，…，ｓ〔ｎ〕，…是強(qiáng)度為λ〔ｔ〕的非齊次泊松過程事件發(fā)生的時(shí)刻。證明顯然，ｓ〔１〕，ｓ〔２〕，…，ｓ〔ｎ〕，…是ｓ１，ｓ２，…，ｓｎ，…的稀疏。設(shè)Ａ＝｛非齊次泊松過程Ｎ〔ｔ〕在〔ｔ，ｔ＋ｈ］中有一個(gè)事件發(fā)生｝，Ｂ＝｛齊次泊松過程Ｎ〔ｔ〕在〔ｔ，ｔ＋ｈ］中有一個(gè)事件發(fā)生｝，則有Ｐ〔ＡＢ〕＝Ｐ〔Ｂ〕Ｐ〔Ａ｜Ｂ〕＝〔λｈ＋ｏ〔ｈ〕〕λ〔ｔ〕/λ=λ〔ｔ〕ｈ＋ｏ〔ｈ〕，由此可知從ｓ１，ｓ２，…，ｓｎ，…中選出的序列ｓ〔１〕，ｓ〔２〕，…，ｓ〔ｎ〕，…．滿足非其次泊松過程的性質(zhì)。根據(jù)定理，先產(chǎn)生齊次泊松過程事件發(fā)生的時(shí)刻，再按概率稀疏就得到非齊次泊松過程事件發(fā)生時(shí)刻，步驟如下．產(chǎn)生參數(shù)λ的齊次泊松過程的Ｔ前事件發(fā)生的時(shí)刻ｓ１，ｓ２，…，ｓｎ．(2)產(chǎn)生〔０，１〕上的隨機(jī)數(shù)ｘｉ，假設(shè)ｘｉ≤λ〔ｓｉ〕／λ，保存ｓｉ，否則舍棄ｓｉ．〔３〕將保存的ｓｉ，分別記為ｓ〔１〕，ｓ〔２〕，…，ｓ〔ｋ〕并輸出即可．.CODEsymstnamdanamda=8.924-1.584*cos(pi*t/1.51)+7.897*sin(pi*t/3.02)-10.434*cos(pi*t/4.53)+4.293*cos(pi*t/6.04);size=1000;%產(chǎn)生{s}的多少times=80;%到達(dá)次數(shù)y=zeros(1,size);z=zeros(1,times);T=zeros(1,times);mu=34;fori=1:1:size*=rand(1);y(i)=-log(*)./mu;%產(chǎn)生{s}endfori=1:1:timesforj=1:1:size*=rand(1);temp=subs(namda,'t',8+y(j));if*<temp/mu%篩選過程z(i)=y(j);Break;endendendT(1)=0;fork=1:1:timesfori=2:kT(i)=T(i-1)+z(i);endendplot(T)*=1:1:80;(b)關(guān)于產(chǎn)生N〔8〕，只需應(yīng)用公式：P{N(t)=n}=e*p〔-λt)*(λt)^n/n!而關(guān)于在〔0,8]的到達(dá)次數(shù)，原理與〔a〕一樣，只需修改代碼的邊界條件。(b).codepartⅠtimes=8;z=zeros(1,100);forj=1:1:80;mu=int(namda,0,j/10);z(j)=e*p(-mu)*(mu)^times/factorial(times);endplot(z)set(handles,'*tick',0:0.1:10);(b).codepartⅡsymstnamdanamda=8.924-1.584*cos(pi*t/1.51)+7.897*sin(pi*t/3.02)-10.434*cos(pi*t/4.53)+4.293*cos(pi*t/6.04);size=1000;%產(chǎn)生{s}的多少times=300;%更改到達(dá)次數(shù)y=zeros(1,size);z=zeros(1,times);T=zeros(1,times);mu=20;fori=1:1:size*=rand(1);y(i)=-log(*)./mu;%產(chǎn)生{s}endfori=1:1:timesforj=1:1:size*=rand(1);temp=subs(namda,'t',8+y(j));if*<temp/mu%篩選過程z(i)=y(j);breakendendendT(1)=0;fork=1:1:timesfori=2:kT(i)=T(i-1)+z(i);endendplot(T)a*is([020008]);%限制時(shí)間[0,8]題目ConsidertheproblemdescribedinE*ample2.3.9.SupposemowthatwehavetwoidenticalHPputerstohandletheiningtraffic.Assumethattheservicetimeofeachputerise*ponentialwitharateof3.5perhour(sotheaggregatetotalservicerateisstill7perhour).Againweassumethatthereare3waitingspaces.Awaitingcustomerwillbeservednythefirstputerthatbeesfreeonafirst-e-first-servedbasis.putethelossprobabilitiesasafunctionoftimetovertheinterval(0,8].PlotyourresultsandpsrethemagianstthoseshowninFigure2.7.2.1基于之前的例題，可以確定解題思路是構(gòu)造關(guān)于Pn(t)的隱式方程組。注意到該題條件的特殊性在于有兩臺(tái)處理器同時(shí)工作。例的推導(dǎo)過程可以借鑒：Po(t+h)=P{*(t+h)=0}=∑P{*(t)=k,*(t+h)=0}.......1隨后將右式展開Po(t+h)=Po(t)[1-λ(t)+o(h)]+P1(t)[μh+o(h)]+o(h)=Po(t)[1-λ(t)]+P1(t)μh+o(h)隨后等式兩端同減Po(t)，并除以h得Po’(t)=-λ(t)Po(t)+μP1(t)....2注意到這里S=2,將1式推廣至n:當(dāng)1<n≤S,Pn(t+h)=P{*(t+h)=n}=∑P{*(t)=k,*(t+h)=n}.....3展開3式Pn(t+h)=Pn-1(t)[λ(t)h+o(h)]+Pn(t)[1-λ(t)h-nμh+o(h)]+Pn+1(t)[(n+1)μh+o(h)]+o(h).再應(yīng)用2式一樣的方法P’n(t+h)=λ(t)Pn-1(t)+(n+1)μhPn+1(t)-[λ(t)+nμ]Pn(t)1<n≤S,P’n(t+h)=λ(t)Pn-1(t)+sμhPn+1(t)-[λ(t)+sμ]Pn(t)n>S最終得到了我需要的用以構(gòu)建隱式方程組的遞推公式基于以上結(jié)論，我以矩陣形式構(gòu)造了方程組，并利用matlabode45解出了{(lán)P}.2.2CODEfunctiony=random3()%主函數(shù)y0=[1,0,0,0,0,0];[t,y]=ode45(odefun,[0,8],y0);%四階-五階Runge-Kutta算法plot(t,y(:,6));%從矩陣中取得我關(guān)心的P5*label('t');ylabel('lossprobability');title('P5');Endfunctiond*=odefun(t,*)%構(gòu)造隱式方程組，子函數(shù)namda=8.924-1.584*cos(pi*t/1.51)+7.897*sin(pi*t/3.02)-10.434*cos(pi*t/4.53)+4.293*cos(pi*t/6.04);mu=3.5;%e*ponetialrateB=[*(1),*(2),*(3),*(4),*(5),*(6)];C=[-namda,mu,0,0,0,0;namda,-(namda+mu),2*mu,0,0,0;0,namda,-(namda+2*mu),2*mu,0,0;0,0,namda,-(namda+2*mu),2*mu,0;0,0,0,namda,-(namda+2*mu),2*mu;0,0,0,0