{"word_count": 580, "character_count": 4798}
Consider a Poisson process , with non-homogeneous intensity . Here, we consider a deterministic function, not a stochastic intensity. Define the cumulated intensityin 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)$valueThe 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)&#62;=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)]&#60;=Tmax){     u=runif(1)     s=s-log(u)     t=min(v[which(Vectorize(Lambda)(v)&#62;=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)&#60;=u){binf=(a+b)/2;bsup=b}       if(Ft((a+b)/2)&#62;=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)]&#60;=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)&#60;=u){binf=(a+b)/2;bsup=b}        if(Ft((a+b)/2)&#62;=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 'sThis 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)&#60;=u){binf=(a+b)/2;bsup=b}       if(Ft((a+b)/2)&#62;=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)*tThe code to generate a Poisson process is   t=0   X=numeric(0)   while(X[length(X)]&#60;=Tmax){     u=runif(1)     t=t-log(u)/lambdau     if(runif(1)&#60;=lambda(t)/lambdau) X=c(X,t)  }The histogram is hereFinally, the last one is also based on a rejection technique, mixed with the second one. I.e. defineThe 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)]&#60;=Tmax){     Ftinvu=function(u) -log(1-x)/lambdau     x=Ftinvu(runif(1))     t=t+x     if(runif(1)&#60;=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...&#160;