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

璋平科技大字電子信息與通信工程學(xué)院實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱 非其次泊松過程課程名稱—隨機(jī)信號分析姓名 顧康 學(xué)號 U201413323日期 6.13 地點(diǎn) 南一樓成績 教師 董燕1.題目ConsiderthenonhomogeneousPoissonprocesswithitsintensityfunctionspectifiedinExample2.3.6.(a)WriteaMATLABprogramtogenerate(stimulate)thefirsteightyarrivaltimes.(b)Givent=8(hours),writeaMatlabprogramtogenerateN(8)andthenthearrivaltimesintheinterval(0,8］,drawtherespectivehistogramsshowinghour5yarrivalcounts.(a)由定理 設(shè)入(t)<A，其中入為一常數(shù)，而sl,s2,...,sn，…為參數(shù)入的齊次泊松過程的事件發(fā)生的時刻，對每個si,以概率入(si)/入進(jìn)行保留，以概率1-A(si)/A舍棄，由此得到的序列s(l),s(2),...,s(n),…是強(qiáng)度為入(t)的非齊次泊松過程事件發(fā)生的時刻。證明顯然，s(1),s(2),...,s(n),..是s1,s2,...,sn,…的稀疏。設(shè)入=｛非齊次泊松過程N(yùn)(t)在(t,t+h］中有一個事件發(fā)生｝，B=｛齊次泊松過程N(yùn)(t)在(t,t+h］中有一個事件發(fā)生｝，則有P(AB)=P(B)P(A|B)=Ah+o(h))A(t)/A(t)h+o(h),由此可知從sl,s2,...,sn,…中選出的序列s(l),s(2),...,s(n)，….滿足非其次泊松過程的性質(zhì)。根據(jù)定理，先產(chǎn)生齊次泊松過程事件發(fā)生的時刻，再按概率稀疏就得到非齊次泊松過程事件發(fā)生時刻，步驟如下.(1)產(chǎn)生參數(shù)入的齊次泊松過程的T前事件發(fā)生的時刻si,s2,..,sn.(2)產(chǎn)生(0,1)上的隨機(jī)數(shù)xi,若xi<A(si)/入，保留si,否則舍棄si.(3)將保留的si，分別記為s(1),s(2),...,s(k)并輸出即可.(a).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:sizex=rand(1);y(i)=-log(x)./mu%產(chǎn)生｛s｝endfori=1:1:timesforj=1:1:sizex=rand(1);temp=subs(namda't',8+y(j));ifx<temp/mu%篩選過程z(i)=y(j);Break;endendendT(1)=0;fork=1:1:timesfori=2:kT(i)=T(i-1)+z(i);endendplot(T)X=1:1:80;(b)關(guān)于產(chǎn)生N(8)，只需應(yīng)用公式:

P{N(t)=n}=exp(-入t)*(入t)人n/n!而關(guān)于在（0,8］內(nèi)的到達(dá)次數(shù)，原理與（a）相同，只需修改代碼的邊界條件。(b).codepartItimes=8;z=zeros(1,100);forj=1:1:80;mu=int(namda,0,j/10);z(j)=exp(-mu)*(mu)八times/factorial(times);endplot(z)set(handles,'xtick',0:0.1:10);0.14fcozW輯EQ0.14fcozW輯EQ」d(b).codepartnsymstnamdanamda=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:sizex=rand(1);y(i)=-log(x)./mu%產(chǎn)生{s}endfori=1:1:timesforj=1:1:sizex=rand(1);temp=subs(namda',t',8+y(j));ifx<temp/mu%篩選過程z(i)=y(j);breakendendendT(1)=0;fork=1:1:timesfori=2:kT(i)=T(i-1)+z(i);endendplot(T)axis([020008]);%限制時間[0,8]

arrivaltimesarrivaltimes2.題目ConsidertheproblemdescribedinExample2.3.9.SupposemowthatwehavetwoidenticalHPcomputerstohandletheincomingtraffic.Assumethattheservicetimeofeachcomputerisexponentialwitharateof3.5perhour(sotheaggregatetotalservicerateisstill7perhour).Againweassumethatthereare3waitingspaces.Awaitingcustomerwillbeservednythefirstcomputerthatbecomesfreeonafirst-come-first-servedbasis.Computethelossprobabilitiesasafunctionoftimetovertheinterval(0,8].PlotyourresultsandcompsrethemagianstthoseshowninFigure基于之前的例題，可以確定解題思路是構(gòu)造關(guān)于Pn(t)的隱式方程組。注意到該題條件的特殊性在于有兩臺處理器同時工作。例2.3.6的推導(dǎo)過程可以借鑒：Po(t+h)=P{X(t+h)=0}='P{X(t)=k,X(t+h)=0}.......1隨后將右式展開Po(t+h)=Po(t)[1-入(t)+o(h)]+P1(t)[ph+o(h)]+o(h)=Po(t)[1-入(t)]+P1(t)ph+o(h)隨后等式兩端同減Po(t)，并除以h得Po’(t)=-入(t)Po(t)+MP1(t)....2注意到這里S=2,將1式推廣至n:當(dāng)1<n<S,Pn(t+h)=P{X(t+h)=n}=JP{X(t)=k,X(t+h)=n}.....3展開3式Pn(t+h)=Pn-1(t)[入(t)h+o(h)]+Pn(t)[1-A(t)h-nph+o(h)]+Pn+1(t)[(n+1)Mh+o(h)]+o(h).再應(yīng)用2式相同的方法P’n(t+h)=A(t)Pn-1(t)+(n+1)MhPn+1(t)-[A(t)+np]Pn(t) 1<n<S,P’n(t+h)=A(t)Pn-1(t)+SMhPn+1(t)-[A(t)+sp]Pn(t) n>S最終得到了我需要的用以構(gòu)建隱式方程組的遞推公式基于以上結(jié)論，我以矩陣形式構(gòu)造了方程組，并利用matlabode45解出T{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)心的P5xlabel('t');ylabel('lossprobability');title('P5');Endfunctiondx=odefun(t,x)%構(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;%exponetialrateB=[x(1),x(2),x(3), x(4),x(5),x(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+