Consider a Poisson process , with non-homogeneous intensity . Here, we consider a deterministic function, not a stochastic intensity. Define the cumulated intensity

in the sense that the number of events that occurred between time and is a random variable that is Poisson distributed with parameter .

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 .

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

- start with
- generate
- set
- set denote
- deliver
- go to step 2.

In order to get the infinimum of , 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 ,

- start with
- generate
- set
- deliver
- go to step 2.

Here the algorithm is simple. For the computational side, at each step, we have to compute and then . To do so, since is increasing with values in , 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.

- generate the number of events on the time interval
- generate independently where
- set i.e. the ordered values
- deliver ‘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

- start with
- generate
- set
- generate (independent of )
- if then deliver
- 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

The good thing is that this function can easily be inverted

- start (as usual) with
- generate
- set
- generate
- if then deliver
- 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://web.ics.purdue.edu/~pasupath/…. In the paper, non-homogeneous spatial Poisson processes are also mentioned…

Hi, thank you for the post. I read it couple of months ago, and was using one of the algorithms in one of my simulation studies. However, I was able to find a minor mistake in one of the codes you provided (which led to problematic results in my simulation study, and that’s why I found it). Below is the code I am talking about:

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){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)

The minor problem with this code is that if n happens to be 0, then the “for(i in 1:n)” loop does not give what we need.

Just though I should post it in case others are in need. 🙂

Thanks for the post again.