dat <- list(step = trex$step[1:500], angle = trex$angle[1:500])
nll <- function(par) {
getAll(par, dat)
Gamma <- tpm(eta)
allprobs <- statedens(step, dgamma2, exp(log_mu), exp(log_sigma)) *
statedens(angle, dvm, 0, exp(log_kappa))
-forward(stationary(Gamma), Gamma, allprobs)
}
par <- list(eta = c(-2, -2), log_mu = log(c(0.3, 2.5)), log_sigma = log(c(0.2, 1.5)),
log_kappa = log(c(0.2, 1.5)))
obj <- MakeADFun(nll, par, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
mod <- lama(obj)
plot(mod) # pseudo-residuals
plot(mod, type = "statedens") # state-dependent densities of both streams
plot(mod, type = "statedens", stream = "step", weights = "viterbi",
breaks = 40, xlab = "step length")
Run the code above in your browser using DataLab