Generating a non-homogeneous Poisson process

Consider a Poisson process gif.latex (54×20), 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 gif.latex (8×13) and gif.latex (6×12) 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 https://latex.codecogs.com/gif.latex?s=0
  2. generate gif.latex (96×19)
  3. set gif.latex (112×19)
  4. set gif.latex (6×12) denote gif.latex (124×19)
  5. deliver
  6. go to step 2.

In order to get the infinimum of gif.latex (12×13), 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 gif.latex (51×16)
  3. set gif.latex (74×14)
  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 http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=F_t%5E%7B-1%7D. 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 gif.latex (101×19)
  2. generate independently gif.latex (114×17) where 
  3. set gif.latex (60×15) i.e. the ordered values  gif.latex (136×16)
  4. deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?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
  2. generate gif.latex (96×19)
  3. set gif.latex (137×19)
  4. generate gif.latex (95×19) (independent of http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=u)
  5. if gif.latex (90×19) then deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?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

gif.latex (433×20)

The good thing is that this function can easily be inverted

gif.latex (215×21)

  1. start (as usual) with
  2. generate gif.latex (63×19)
  3. set gif.latex (74×14)
  4. generate gif.latex (96×19)
  5. if gif.latex (124×19) then deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?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://web.ics.purdue.edu/~pasupath/…. In the paper, non-homogeneous spatial Poisson processes are also mentioned…

 


3 thoughts on “Generating a non-homogeneous Poisson process”

  1. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.