# 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://web.ics.purdue.edu/~pasupath/…. In the paper, non-homogeneous spatial Poisson processes are also mentioned…

# Actuariat IARD

Cet hiver (même si la nouvelle ne sera officielle qu’à la rentrée), je devrais donner le cours ACT2040, actuariat IARD. Le plan de cours sera bientôt en ligne, mais je peux déjà dire que le cours sera basé sur le Tome 2 du livre écrit avec Michel Denuit il y a quelques années, mathématiques de l’assurance non-vie. Le cours est une suite du cours ACT6420 méthodes de prévisions, donné cet automne (qui est un prérequis indiqué sur le site du registrariat http://websysinfo.uqam.ca/…): je partirais donc du fait que le modèle linéaire de régression est connu (et compris) et que tout le monde sait utiliser R, et lire des sorties de régression. Mais les premières démonstrations reviendront sur l’utilisation de R, et sur l’analyse de la variance, que l’on n’a pas vraiment eu le temps d’aborder dans le cours de régression. Pour des références sur R, je conseille

• “R pour les débutants” d’Emmanuel Paradis, (PDF)
• “Introduction à la programmation en S” par Vincent Goulet, (PDF)

pour les documents en français, ou pour des documents plus complets, mais en anglais

• “R for Beginners” d’Emmanuel Paradis (PDF),
• “An Introduction to R” par Longhow Lam (PDF)
• “The R language — a short companion” par Marc Vandemeulebroecke (PDF),
• “The R Guide” par Jason Owen (PDF),
• “Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
• “Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
• “Statistics with R and S-Plus” d’Hugo Quené (PDF)
• “Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
• “Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Sinon, les transparents du premier cours sont en ligne ici, et le plan de cours est

et je mettrais bientôt en ligne des liens vers des bases de données que l’on utilisera tout au long du cours, ou en démonstration.

Pour les références, je citerais deux livres sur lesquels je m’appuierai beaucoup car je les connais presque par cœur. Ils sont disponibles à la Coopuqam