复合泊松过程的实现

电子信息与通信工程学院实验报告实验名称非其次泊松过程课程名称随机信号分析姓名顾康学号U201413323 日期 6.13 地点南一楼成绩教师董燕1.题目Consider the nonhomogeneous Poisson process with its intensity function spectified in Example2.3.6. (a) Write a MATLAB program to generate (stimulate) the first eighty arrival times. (b) Given t=8(hours),write a Matlab program to generate N(8) and then the arrival times in the interval(0,8],draw the respective histograms showing hour5y arrival counts.(a)由定理设λ（ｔ）≤λ，其中λ为一常数，而ｓ１，ｓ２，…，ｓｎ，…为参数λ的齐次泊松过程的事件发生的时刻，对每个ｓｉ，以概率λ（ｓｉ）／λ进行保留，以概率１－λ（ｓｉ）／λ舍弃，由此得到的序列ｓ（１），ｓ（２），…，ｓ（ｎ），…是强度为λ（ｔ）的非齐次泊松过程事件发生的时刻。

证明显然，ｓ（１），ｓ（２），…，ｓ（ｎ），…是ｓ１，ｓ２，…，ｓｎ，…的稀疏。

设Ａ＝｛非齐次泊松过程Ｎ（ｔ）在（ｔ，ｔ＋ｈ］中有一个事件发生｝，Ｂ＝｛齐次泊松过程Ｎ（ｔ）在（ｔ，ｔ＋ｈ］中有一个事件发生｝，则有Ｐ（ＡＢ）＝Ｐ（Ｂ）Ｐ（Ａ｜Ｂ）＝（λｈ＋ｏ（ｈ））λ（ｔ）/λ= λ（ｔ）ｈ＋ｏ（ｈ），由此可知从ｓ１，ｓ２，…，ｓｎ，…中选出的序列ｓ（１），ｓ（２），…，ｓ（ｎ），…．满足非其次泊松过程的性质。

根据定理，先产生齐次泊松过程事件发生的时刻，再按概率稀疏就得到非齐次泊松过程事件发生时刻，步骤如下．（１）产生参数λ的齐次泊松过程的Ｔ前事件发生的时刻ｓ１，ｓ２，…，ｓｎ．( 2 )产生（０，１）上的随机数ｘｉ，若ｘｉ≤λ（ｓｉ）／λ，保留ｓｉ，否则舍弃ｓｉ．（３）将保留的ｓｉ，分别记为ｓ（１），ｓ（２），…，ｓ（ｋ）并输出即可．(a). CODEsyms tnamdanamda=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;%产生{s}的多少times=80; %到达次数y=zeros(1,size);z=zeros(1,times);T=zeros(1,times);mu=34;for i=1:1:sizex=rand(1);y(i)=-log(x)./mu;%产生{s}endfor i=1:1:timesfor j=1:1:sizex=rand(1);temp=subs(namda,'t',8+y(j));if x< temp/mu%筛选过程z(i)=y(j);Break;endendendT(1)=0;for k=1:1:timesfor i=2:kT(i)=T(i-1)+z(i);endendplot(T)X=1:1:80;(b)关于产生N（8），只需应用公式：P{N(t)=n}=exp（-λt)* (λt)^n/n!而关于在（0,8]内的到达次数，原理与（a）相同，只需修改代码的边界条件。

(b).code partⅠtimes=8;z=zeros(1,100);for j=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);(b).code part Ⅱsyms tnamdanamda=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;%产生{s}的多少times=300;%更改到达次数y=zeros(1,size);z=zeros(1,times);T=zeros(1,times);mu=20;for i=1:1:sizex=rand(1);y(i)=-log(x)./mu;%产生{s}endfor i=1:1:timesfor j=1:1:sizex=rand(1);temp=subs(namda,'t',8+y(j)); if x< temp/mu%筛选过程z(i)=y(j);breakendendendT(1)=0;for k=1:1:timesfor i=2:kT(i)=T(i-1)+z(i);endendplot(T)axis([0 200 0 8]);%限制时间[0,8]2.题目Consider the problem described in Example 2.3.9.Suppose mow that we have two identical HP computers to handle the incoming traffic.Assume that the service time of each computer is exponential with a rate of 3.5 per hour(so the aggregate total service rate is still 7 per hour ).Again we assume that there are 3 waiting spaces.A waiting customer will be served ny the first computer that becomes free on a first-come-first-served pute the loss probabilities as a function of time t over the interval (0,8].Plot your result s and compsre them agianst those shown inFigure 2.7.2.1基于之前的例题，可以确定解题思路是构造关于Pn(t)的隐式方程组。

注意到该题条件的特殊性在于有两台处理器同时工作。

例2.3.6的推导过程可以借鉴：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)[μ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:当1<n≤S,P n(t+h)=P{X(t+h)=n}=∑P{X(t)=k,X(t+h)=n} (3)展开3式P n(t+h)=P n-1(t)[λ(t)h+o(h)]+P n(t) [1-λ(t)h-nμh+o(h)]+P n+1(t) [(n+1)μh+o(h)]+o(h).再应用2式相同的方法P’n(t+h)=λ(t)P n-1(t)+(n+1)μhP n+1(t)- [λ(t)+nμ]P n(t) 1<n≤S,P’n(t+h)=λ(t)P n-1(t)+sμhP n+1(t)- [λ(t)+sμ]P n(t) n>S最终得到了我需要的用以构建隐式方程组的递推公式基于以上结论，我以矩阵形式构造了方程组，并利用matlab ode45 解出了{P}.2.2CODEfunction y=random3()%主函数y0=[1,0,0,0,0,0];[t,y]=ode45(@odefun,[0,8],y0);%四阶-五阶Runge-Kutta算法plot(t,y(:,6));%从矩阵中取得我关心的P5xlabel('t');ylabel('loss probability');title('P5');Endfunction dx=odefun(t,x)%构造隐式方程组，子函数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;%exponetial rateB=[ 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+2*mu),2*mu,0;0, 0 , 0, namda,-(namda+2*mu),2*mu;0, 0 , 0, 0, namda,-2*mu ];dx=C*B';%复杂的方程组简化为矩阵运算end8am 4pm8pm 4pm对比λ(t) 与P5，峰值同样在相近的时间达到参考文献：1.《非齐次泊松过程的仿真方法》-宁如云。