Example 1: Markov-modulated Poisson processes
Setting parameters
We choose to have a considerably higher rate and shorter stays of the underlying Markov chain in state 2, i.e. state 2 is bursty.
Simulating an MMPP
set.seed(123)
k = 200 # number of state switches
trans_times = s = rep(NA, k) # time points where the chain transitions
s[1] = sample(1:2, 1) # initial distribuion c(0.5, 0.5)
# exponentially distributed waiting times
trans_times[1] = rexp(1, -Q[s[1],s[1]])
# in a fixed interval, the number of arrivals is Pois(lambda * interval_length)
n_arrivals = rpois(1, lambda[s[1]]*trans_times[1])
# arrival times within fixed interval are uniformly distributed
arrival_times = runif(n_arrivals, 0, trans_times[1])
for(t in 2:k){
s[t] = c(1,2)[-s[t-1]] # for 2-states, always a state swith when transitioning
# exponentially distributed waiting times
trans_times[t] = trans_times[t-1] + rexp(1, -Q[s[t], s[t]])
# in a fixed interval, the number of arrivals is Pois(lambda * interval_length)
n_arrivals = rpois(1, lambda[s[t]]*(trans_times[t]-trans_times[t-1]))
# arrival times within fixed interval are uniformly distributed
arrival_times = c(arrival_times,
runif(n_arrivals, trans_times[t-1], trans_times[t]))
}
arrival_times = sort(arrival_times)Let’s visualize the simulated MMPP
n = length(arrival_times)
color = c("orange", "deepskyblue")
plot(arrival_times[1:100], rep(0.5,100), type = "h", bty = "n", ylim = c(0,1),
yaxt = "n", xlab = "arrival times", ylab = "")
segments(x0 = c(0,trans_times[1:98]), x1 = trans_times[1:99],
y0 = rep(0,100), y1 = rep(0,100), col = color[s[1:99]], lwd = 4)
legend("top", lwd = 2, col = color, legend = c("state 1", "state 2"), box.lwd = 0)
What makes the MMPP special compared to a regular Poisson point process
is its burstiness when the Markov chain is in the
second state.
Writing the negative log-likelihood function
The likelihood of a stationary MMPP for waiting times \(x_1, \dots, x_n\) is (Meier-Hellstern (1987), Langrock, Borchers, and Skaug (2013)) \[
L(\theta) = \delta \bigl(\prod_{i=1}^n
\exp((Q-\Lambda)x_i)\Lambda\bigr)1,
\] where \(Q\) is the generator
matrix of the continuous-time Markov chain, \(\Lambda\) is a diagonal matrix of
state-dependent Poisson intensities, \(\delta\) is the stationary distribution of
the continuous-time Markov chain, and \(1\) is a column vector of ones. For more
details on continuous-time Markov chains, see the vignette
continuous-time HMMs or also Dobrow (2016). We can easily
calculate the log of the above expression using the standard
implementation of the general forward algorithm forward_g()
when choosing the first matrix of state-dependent densities to be the
identity (i.e.) the first row of the allprobs matrix to be
one and all other matrices of state-dependent density matrices to be
\(\Lambda\).
mllk = function(theta.star, timediff, N=2){
lambda = exp(theta.star[1:N]) # state specific rates
Q = diag(N) # generator matrix
Q[!Q] = exp(theta.star[N+1:(N*(N-1))])
diag(Q) = 0
diag(Q) = -rowSums(Q)
Qube = LaMa::tpm_cont(Q-diag(lambda), timediff) # exp((Q-Lambda)*deltat)
allprobs = matrix(lambda, nrow = length(timediff+1), ncol = N, byrow = T)
allprobs[1,] = 1
delta = solve(t(Q+1), rep(1,N), tol = 1e-20)
-LaMa::forward_g(delta, Qube, allprobs)
}Fitting an MMPP to the data
theta.star = log(c(2, 15, # lambda
2, 0.5)) # off-diagonals of Q
timediff = diff(arrival_times)
t1 = Sys.time()
mod = nlm(mllk, theta.star, timediff=timediff, stepmax = 10)
# we often need the stepmax, as the matrix exponential can be numerically unstable
Sys.time()-t1
#> Time difference of 0.07354617 secs
