## Simulate a NON-Homogeneous Poisson Processes (NHPP)
## Here we follow the computationally slower but more intuitive approach
## from Lews and Shedler, 1979. See also "Algorithm 2" in the Thinning notes
## online at https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf
## Those notes also have deeper coverage of these topics, if you're interested.
## The idea here is to first simulate a HPP with rate rmax
## which is the maximal value of the intensity function r(t) over [0,L]
## Then for each event generated under that HPP, we keep an event at time t
## with probability p=r(t)/rmax. Thinking about this in terms of our
## "small discrete bins each with probability p=r*dt" approximation of an HPP
## all we're doing for the NHPP analog is replacing the
## single event probability p(t)=r(t)*dt with a two stage process of first
## picking whether or not an event happens with prob1=rmax*dt, then
## using another Bernoilli to decide whether or not to keep it with prob2=r(t)/rmax.
## This gives us events with total prob p(t)=prob1*prob2=r(t)*dt, as desired.
## Set random number generator seed...
set.seed(1880)
## First, we simulate our HPP...
# Suppose we had intensity function rvar(t) that was reven or rodd
# depending on whether floor(t) was even or odd...
rvar <- function(s,reven=10,rodd=1) { ifelse(floor(s) %% 2 == 0, reven, rodd) }
## Consider the interval [0,L]
L=6 # interval length
## Define rmax (to reuse old code, call it r...)
# r = 5 is too easy ;) In general, let's let optimize() find our maximum! :)
ropt = optimize(rvar, interval = c(0,L), maximum = TRUE)
ropt
r <- ropt$objective
r
## HPP with rate r=rmax over interval [0,L]
Ts=c() # initialize vector of inter-event times
X=c() # initialize vector event times
i=1 # index for while loop
xlast = 0 # initialize most recent event time (absolute time)
while(xlast < L) { # while our last event is not yet past L...
## Step 1, generate a vector of inter-event times
## by drawing iid values from an exponential dist. with rate=r.
## Ts = that vector of times.
Ts[i] <- rexp(1,rate=r)
## Step 2, get "absolute" even times from inter-event times
## by addition...
## X[1] = Ts[1]
## X[2] = Ts[1] + Ts[2]
## X[3] = Ts[1] + Ts[2] + Ts[3]
## ...
## This can be achieved with a cummulative sum:
## X = cumsum(Ts)
##
## In this while loop, just keep track of the most recent.
xlast <- xlast + Ts[i]
i <- i + 1; # update our counter for the while loop...
} # end while loop.
## Fill out absolute even times
X0 = cumsum(Ts)
X0
## remove the last one, which is > L
X0 <- X0[-length(X0)] # if length(x) is 10, X0[-10] is the same as X0[1:9]
################
## STEP 2 -- Thin our HPP events t1, t2, ... with probability rvar(ti)/rmax
prob <- function(s) {rvar(s)/r} # our probability function rvar(s)/rmax
# Draw Bernoulli r.v.s one at a time via rbinom(1,size=1,p=prob(X[i])) or
# Draw Bernoulli r.v.s all at once via rbinom(length(X),size=1,p=prob(X[i])) or
# Draw Bernoilli r.v.s all at once by sampling from a uniform then seeing
# if those values fall above or below "prob" for each event... we'll do #2
keep <- as.logical(rbinom(length(X0),1,prob(X0))) # make 0/1 into False/True...
keep
X <- X0[keep]
## Plot the results
curve(rvar(x),0,L,ylim=c(0,1.1*r), n=1000)
abline(h=0)
points(X0,X0*0, pch='|', col="gray", lwd=2) # Plot our HPP events in gray on [0,L]
points(X, X*0, pch='|', col="black", lwd=2) # Plot our NHPP events in black on [0,L]
#################################################################################
## An alternative way to simulate NHPPs is to integrate the intensity function
## as part of calculating the CDF for the first even time, then sample that
## using a trick: Step 1 sample u=runif(1) from a uniform[0,1], then
## Step 2 calculate the inverse CDF value of u, which
## is equivalent to sampling from the r.v. with that CDF.
## This can be computationally faster, but is a bit more complicated to code.