Simulation example
We begin by considering homogeneous HSMMs. In such models, each state has an associated state dwell-time distribution. The transition probability matrix of a regular HMM is replaced by these distributions and the conditional transition probabilities given the state is left.
Setting parameters
Here we choose the simplest case of dwell times that are shifted Poisson distributed. We have to specify the Poisson mean for each state, the conditional transition probability matrix called \(\Omega\) and the parameters of the state-dependent process.
lambda = c(7, 4, 4)
omega = matrix(c(0, 0.7, 0.3,
0.5, 0, 0.5,
0.7, 0.3, 0), nrow = 3, byrow = TRUE)
mu = c(10, 40, 100)
sigma = c(5, 20, 50)
color = c("orange", "deepskyblue", "seagreen2")
curve(dnorm(x, mu[1], sigma[1]), lwd = 2, col = color[1], bty = "n",
xlab = "x", ylab = "density", xlim = c(0, 150), n = 300)
curve(dnorm(x, mu[2], sigma[2]), lwd = 2, col = color[2], add = T)
curve(dnorm(x, mu[3], sigma[3]), lwd = 2, col = color[3], add = T)
Simulating data
We simulate data by drawing dwell times from the dwell-time distribution of the current state and then draw the next state using the conditional transition probabilities. The state-dependent process is drawn conditional on the current state.
set.seed(123)
k = 50 # number of stays
s = rep(NA, k)
s[1] = sample(1:3, 1) # uniform initial distribution
staylength = rpois(1, lambda[s[1]]) + 1 # drawing dwell time from shifted Poisson
C = rep(s[1], staylength)
x = rnorm(staylength, mu[s[1]], sigma[s[1]])
for(t in 2:k){
# conditionally drawing state
s[t] = sample(c(1:3)[-s[t-1]], 1, prob = omega[s[t-1], -s[t-1]])
staylength = rpois(1, lambda[s[t]]) + 1 # drawing dwell time from shifted Poisson
C = c(C, rep(s[t], staylength))
x = c(x, rnorm(staylength, mu[s[t]], sigma[s[t]]))
}
plot(x, pch = 20, col = color[C], bty = "n")
legend("topright", col = color, pch = 20,
legend = paste("state", 1:3), box.lwd = 0)
Writing the negative log-likelihood function
We now write the negative log-likelihood function for an
approximating HMM. As a semi-Markov chain is specified in terms of
state-specific dwell-time distributions and conditional transition
probabilities given that the current state is left, we have to compute
both (here called dm and omega). For the
latter, we can use the function tpm_emb() that constructs a
transition probability matrix via the inverse multinomial logit link
(softmax), where the diagonal entries are forced to equal zero.
The transition probability matrix of the approxmiating HMM can then
be computed by the function tpm_hsmm() where the exact
procedure is detailed by Langrock and Zucchini
(2011). We need the extra
argument agsizes to specify the aggregate sizes that should
be used to approximate the dwell-time distributions. These should be
chosen such that most of the support of the state-specific dwell-time
distributions is covered.
nll = function(par, x, N, agsizes){
mu = par[1:N]
sigma = exp(par[N+1:N])
lambda = exp(par[2*N+1:N])
omega = if(N==2) tpm_emb() else tpm_emb(par[3*N+1:(N*(N-2))])
dm = list() # list of dwell-time distributions
for(j in 1:N) dm[[j]] = dpois(1:agsizes[j]-1, lambda[j]) # shifted Poisson
Gamma = tpm_hsmm(omega, dm, sparse = FALSE)
delta = stationary(Gamma)
allprobs = matrix(1, length(x), N)
ind = which(!is.na(x))
for(j in 1:N){
allprobs[ind,j] = dnorm(x[ind], mu[j], sigma[j])
}
-forward_s(delta, Gamma, allprobs, agsizes)
}Fitting an HSMM (as an approxiating HMM) to the data
# intial values
par = c(10, 40, 100, log(c(5, 20, 50)), # state-dependent
log(c(7,4,4)), # dwell time means
rep(0, 3)) # omega
agsizes = qpois(0.95, lambda)+1
system.time(
mod <- nlm(nll, par, x = x, N = 3, agsizes = agsizes, stepmax = 2)
)
#> user system elapsed
#> 0.455 0.006 0.462Fitting HSMMs is rather slow (even using C++) as we translate the additional model complexity into a higher computational overhead (31 states here).
Results
N = 3
(mu = mod$estimate[1:N])
#> [1] 10.16569 39.06161 107.66034
(sigma = exp(mod$estimate[N+1:N]))
#> [1] 4.78882 19.35639 48.56115
(lambda = exp(mod$estimate[2*N+1:N]))
#> [1] 6.942983 4.595469 3.354765
(omega = tpm_emb(mod$estimate[3*N+1:(N*(N-2))]))
#> [,1] [,2] [,3]
#> [1,] 0.0000000 0.5541030 0.4458970
#> [2,] 0.5040937 0.0000000 0.4959063
#> [3,] 0.6654703 0.3345297 0.0000000