Generating a non-homogeneous Poisson process

Consider a Poisson process $(N_t)_{t\geq 0}$, with non-homogeneous intensity $\lambda(t)$. Here, we consider a deterministic function, not a stochastic intensity. Define the cumulated intensity

$\Lambda(t)=\int_0^t\lambda(s)ds$

in the sense that the number of events that occurred between time $0$ and $t$ is a random variable that is Poisson distributed with parameter  $\Lambda(t)$.

For example, consider here a cyclical Poisson process, with intensity

   lambda=function(x) 100*(sin(x*pi)+1)

To compute the cumulated intensity, consider a very general function

   Lambda=function(t) integrate(f=lambda,lower=0,upper=t)\$value

The idea is to generate a Poisson process on a finite interval $[0,T]$.

The first code is based on a proposition from Çinlar (1975),

1. start with $s=0$
2. generate $u\sim\mathcal{U}([0,1])$
3. set $s\leftarrow s-\log(u)$
4. set $t$ denote $\inf\{v;\Lambda(v)>s\}$
5. deliver $t$
6. go to step 2.

In order to get the infinimum of $\Lambda$, consider a code as

   v=seq(0,Tmax,length=1000)
t=min(v[which(Vectorize(Lambda)(v)>=s)])

(it might not be very efficient…. but it should work). Here, the code to generate that Poisson process is

   s=0; v=seq(0,Tmax,length=1000)
X=numeric(0)
while(X[length(X)]<=Tmax){
u=runif(1)
s=s-log(u)
t=min(v[which(Vectorize(Lambda)(v)>=s)])
X=c(X,t)
}

Here, we get the following histogram,

   hist(X,breaks=seq(0,max(X)+1,by=.1),col="yellow")
u=seq(0,max(X),by=.02)
lines(u,lambda(u)/10,lwd=2,col="red")

Consider now another strategy. The idea is to use the conditional distribution before the next event, given that one occurred at time $t$,

$F_t(x)=1-\exp\left(\Lambda(x)-\Lambda(x+t)\right)=1-\exp\left(-\int_t^{t+x}\lambda(s)ds\right)$

1. start with $t=0$
2. generate $x\sim F_t$
3. set $t\leftarrow t+x$
4. deliver $t$
5. go to step 2.

Here the algorithm is simple. For the computational side, at each step, we have to compute $F_t$ and then $formdata=F_t^{-1}$. To do so, since $F_t$ is increasing with values in $[0,1]$, we can use a dichotomic algorithm,

   Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t))
Ftinv=function(u){
a=0
b=Tmax
for(j in 1:20){
if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
a=binf
b=bsup
}
return((a+b)/2)
}

Here the code is the following

   t=0; X=t
while(X[length(X)]<=Tmax){
Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t))
Ftinv=function(u){
a=0
b=Tmax
for(j in 1:20){
if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
a=binf
b=bsup
}
return((a+b)/2)
}
x=Ftinv(runif(1))
t=t+x
X=c(X,t)
}

The third code is based on a classical algorithm to generate an homogeneous Poisson process on a finite interval: first, we generate the number of events, then, we draw uniform variates, and we sort them. Here, the strategy is closed, except that is won’t be uniform any longer.

1. generate the number of events on the time interval $n\sim\mathcal{P}(\Lambda(T))$
2. generate independently $z_1,\cdots,z_n\sim F$ where $F(t)=\Lambda(t)/\Lambda(T)$
3. set $t_i = z_{i:n}$ i.e. the ordered values  $t_1\leq t_2\leq \cdots\leq t_n$
4. deliver $formdata=t_i$‘s

This algorithm is extremely simple, and also very fast. This is one function to inverse, and it is not in the loop,

   n=rpois(1,Lambda(Tmax))
Ft=function(x) Lambda(x)/Lambda(Tmax)
Ftinv=function(u){
a=0
b=Tmax
for(j in 1:20){
if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
a=binf
b=bsup
}
return((a+b)/2)
}
X0=rep(NA,n)
for(i in 1:n){
X0[i]=Ftinv(runif(1))
}
X=sort(X0)

Here is the associated histogram,

An alternative is based on a rejection technique. Actually, it was the algorithm mentioned a few years ago on this blog (well, the previous one). Here, we need an upper bound for the intensity, so that computations might be much faster. Here, consider

1. start with $t=0$
2. generate $u\sim\mathcal{U}([0,1])$
3. set $t\leftarrow t-\log(x)/\lambda_u$
4. generate $v\sim\mathcal{U}([0,1])$ (independent of $formdata=u$)
5. if $v\leq\lambda(t)/\lambda_u$ then deliver $formdata=t$
6. go to step 2.

Here, consider a constant upper bound,

   lambdau=function(t) 200
Lambdau=function(t) lambdau(t)*t

The code to generate a Poisson process is

   t=0
X=numeric(0)
while(X[length(X)]<=Tmax){
u=runif(1)
t=t-log(u)/lambdau
if(runif(1)<=lambda(t)/lambdau) X=c(X,t)
}

The histogram is here

Finally, the last one is also based on a rejection technique, mixed with the second one. I.e. define

$F_{t,u}(x)=1-\exp\left(\Lambda_u(x)-\Lambda_u(x+t)\right)=1-\exp\left(-x\lambda_u\right)$

The good thing is that this function can easily be inverted

$F_{t,u}(x)^{-1}=-\log(1-x)/\lambda_u$

1. start (as usual) with $t=0$
2. generate $x\sim F_{t,u}$
3. set $t\leftarrow t+x$
4. generate $u\sim\mathcal{U}([0,1])$
5. if $u\leq \lambda(t+x)/\lambda_u$ then deliver $formdata=t$
6. goto step 2.

Here, the algorithm is simply

   t=0
while(X[length(X)]<=Tmax){
Ftinvu=function(u) -log(1-x)/lambdau
x=Ftinvu(runif(1))
t=t+x
if(runif(1)<=lambda(t+x)/lambdau(t+x)) X=c(X,t)
}

Obviously those five codes work, the first one being much slower than the other three. But it might be because my strategy to seek the infimum is not great. And the latter worked well since there were not much rejection, I guess it can be worst…

All those algorithms were mentioned in a nice survey written by Raghu Pasupathy and can be downloaded from http://filebox.vt.edu/… . In the paper, non-homogeneous spatial Poisson processes are also mentioned…

Print This Post