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 visualise 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) = \pi \Bigl(\prod_{i=1}^n \exp\bigl((Q-\Lambda)x_i\bigr)\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, \(\pi\) 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\).
nll = function(par, timediff, N){
lambda = exp(par[1:N]) # state specific rates
Q = generator(par[N+1:(N*(N-1))])
Pi = stationary_cont(Q)
Qube = tpm_cont(Q - diag(lambda), timediff) # exp((Q-Lambda) * dt)
allprobs = matrix(lambda, nrow = length(timediff + 1), ncol = N, byrow = T)
allprobs[1,] = 1
-forward_g(Pi, Qube, allprobs)
}
