## generating data from homogeneous 2-state HSMM
mu = c(0, 6)
beta = matrix(c(log(4),log(6),-0.2,0.2,-0.1,0.4), nrow=2)
# time varying mean dwell time
Lambda = exp(cbind(1, trigBasisExp(1:24, 24, 1))%*%t(beta))
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# simulation
# for a 2-state HSMM the embedded chain always alternates between 1 and 2
s = rep(1:2, 100)
C = x = numeric(0)
tod = rep(1:24, 50) # time of day variable
time = 1
for(t in 1:100){
dt = rpois(1, Lambda[tod[time], s[t]])+1 # dwell time depending on time of day
time = time + dt
C = c(C, rep(s[t], dt))
x = c(x, rnorm(dt, mu[s[t]], 1.5)) # fixed sd 2 for both states
}
tod = tod[1:length(x)]
## negative log likelihood function
mllk = function(theta.star, x, sizes, tod){
# parameter transformations for unconstraint optimization
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE) # omega fixed (2-states)
mu = theta.star[1:2]
sigma = exp(theta.star[3:4])
beta = matrix(theta.star[5:10], nrow=2)
# time varying mean dwell time
Lambda = exp(cbind(1, trigBasisExp(1:24, 24, 1))%*%t(beta))
dm = list()
for(j in 1:2){
dm[[j]] = sapply(1:sizes[j]-1, dpois, lambda = Lambda[,j])
}
Gamma = tpm_phsmm(omega, dm)
delta = stationary_p(Gamma, tod[1])
# calculate all state-dependent probabilities
allprobs = matrix(1, length(x), 2)
for(j in 1:2){ allprobs[,j] = dnorm(x, mu[j], sigma[j]) }
# return negative for minimization
-forward_sp(delta, Gamma, allprobs, sizes, tod)
}
## fitting an HSMM to the data
theta.star = c(1, 4, log(2), log(2), # state-dependent parameters
log(4), log(6), rep(0,4)) # state process parameters dm
mod = nlm(mllk, theta.star, x = x, sizes = c(10, 15), tod = tod, stepmax = 5)
Run the code above in your browser using DataLab