## negative log likelihood function
nll = function(par, step) {
# parameter transformations for unconstrained optimisation
Gamma = tpm(par[1:2]) # multinomial logit link
delta = stationary(Gamma) # stationary HMM
mu = exp(par[3:4])
sigma = exp(par[5:6])
# 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(delta, Gamma, allprobs)
}
## fitting an HMM to the trex data
par = c(-2,-2, # initial tpm params (logit-scale)
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:1000])
Run the code above in your browser using DataLab