## Simple usage
Gamma = array(c(0.9, 0.2, 0.1, 0.8), dim = c(2,2,10))
delta = c(0.5, 0.5)
allprobs = matrix(0.5, 10, 2)
forward_g(delta, Gamma, allprobs)
# \donttest{
## Full model fitting example
## negative log likelihood function
nll = function(par, step, Z) {
# parameter transformations for unconstrained optimisation
beta = matrix(par[1:6], nrow = 2)
Gamma = tpm_g(Z, beta) # multinomial logit link for each time point
delta = stationary(Gamma[,,1]) # stationary HMM
mu = exp(par[7:8])
sigma = exp(par[9:10])
# calculate all state-dependent probabilities
allprobs = matrix(1, length(step), 2)
ind = which(!is.na(step))
for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
# simple forward algorithm to calculate log-likelihood
-forward_g(delta, Gamma, allprobs)
}
## fitting an HMM to the trex data
par = c(-1.5,-1.5, # initial tpm intercepts (logit-scale)
rep(0, 4), # initial tpm slopes
log(c(0.3, 2.5)), # initial means for step length (log-transformed)
log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:500], Z = trigBasisExp(trex$tod[1:500]))
# }
Run the code above in your browser using DataLab