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

1. start with
2. generate 
3. set 
4. set  denote 
5. deliver
6. 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 ,

1. start with
2. generate 
3. set 
4. deliver
5. 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.

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

1. start with
2. generate 
3. set 
4. generate  (independent of )
5. if  then deliver
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

The good thing is that this function can easily be inverted

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