Non-homogeneous Poisson processes

Introduction

For many applications, it is obvious that one should consider a rate \(\lambda = \lambda(t)\) that depends on time; the expected number arrivals during to a coffee shop will most likely depend on the time. Here we are interested in simulating such Poisson processes by thinning as constructed by Lewis and Shedler, 1979. Recall that if we are given a Poisson process of intensity \(4\), and we keep each arrival independently with probability \(\tfrac{2}{3}\), we are left with a Poisson process of intensity \(\tfrac{8}{3}\). This same idea can be used to generate a non-homogeneous Poisson point process of intensity \(\lambda(t)\), where \(0\leq \lambda(t) \leq M\) is uniformly bounded by \(M\).

• First generate a homogeneous Poisson process of intensity \(M\)
• At each point/arrival \(t\) that is generated, we keep it independently with probability \[ \frac{\lambda(t)}{M}\]

Exercise

Suppose that we are given \(\lambda(t) = \sin(t) + 4\). Generate a nonhomogeneous Poisson point process on \([0, 200]\) with intensity function \(\lambda\).

Solution

First, we generate the homogeneous Poisson point process of intensity \(5\).

arr <- function(end){
T = rexp(1,5)
while(T[length(T)] < end){
  T <- c(T,  (T[length(T)] + rexp(1,5))  )
}
T <- T[-length(T)]   # remove the last arrival, as it will be larger than 10
T
}

z=arr(10)  # outputs the arrival times of a Poisson process up to time 10

Next, we write a function that applies the thinning

thinning <- function(x){
  del=NULL
  for( i in 1:length(x) ){
    t= x[i]
    prop = (sin(t)+4)/5
    coin = rbinom(1,1,prop)
    if(coin==0){
      del = c(del, i)
    }
  }
  del    # outputs the coordinates to be deleted.
}
z

##  [1] 0.06111844 0.20237217 0.28589158 0.41150918 0.47386211 0.51888097
##  [7] 0.79226588 1.22207402 1.32100734 1.33561335 1.44895489 1.46404218
## [13] 1.69482899 1.71427886 1.83646057 1.94468428 2.14769464 2.15689143
## [19] 2.31178626 2.31838113 2.39242229 3.08138784 3.12571870 3.21549634
## [25] 3.29419737 3.82298970 4.19496224 4.40147581 4.55191990 4.86309012
## [31] 5.02860231 5.21780134 5.23140555 5.38115596 5.84352102 5.97550751
## [37] 6.08430644 6.62265717 6.74499760 6.98281261 7.01010646 7.05563083
## [43] 7.32028113 7.53006552 7.70238402 7.72230050 7.77401133 7.85485737
## [49] 8.00314313 8.09975031 8.72892179 8.81497241 9.68991499

thin=thinning(z)
nonhomg = z[-thin]    # the corresponding entries deleted.
thin

##  [1]  3  4  5 19 22 26 29 36 37 40 53

nonhomg

##  [1] 0.06111844 0.20237217 0.51888097 0.79226588 1.22207402 1.32100734
##  [7] 1.33561335 1.44895489 1.46404218 1.69482899 1.71427886 1.83646057
## [13] 1.94468428 2.14769464 2.15689143 2.31838113 2.39242229 3.12571870
## [19] 3.21549634 3.29419737 4.19496224 4.40147581 4.86309012 5.02860231
## [25] 5.21780134 5.23140555 5.38115596 5.84352102 6.62265717 6.74499760
## [31] 7.01010646 7.05563083 7.32028113 7.53006552 7.70238402 7.72230050
## [37] 7.77401133 7.85485737 8.00314313 8.09975031 8.72892179 8.81497241

length(z)

## [1] 53

length(thin)

## [1] 11

length(nonhomg)

## [1] 42