## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(funHMM)

## ----bmwd-sim-----------------------------------------------------------------
set.seed(137)
A1 <- matrix(0.09, 5, 5) + 0.55 * diag(5)
low <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A1,
             type = "bmwd", par = c(-4, -2, 0, 2, 4), len = 100)
matplot(t(low$data), type = "l", lty = 1, col = low$states + 1,
        xlab = "grid point", ylab = "", main = "Brownian motion with drift")

## ----bmwd-fit-----------------------------------------------------------------
fit.low <- thmm(low$data, nstates = 5, type = "bmwd", tol = 1e-8)
fit.low

## ----bmwd-eval----------------------------------------------------------------
table(estimated = fit.low$states, truth = low$states)
ari(fit.low$states, low$states)

## ----bmwd-med-----------------------------------------------------------------
set.seed(137)
med <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A1,
             type = "bmwd", par = c(-8, -4, 0, 4, 8), len = 100)
fit.med <- thmm(med$data, nstates = 5, type = "bmwd", tol = 1e-8)
sort(fit.med$par)
ari(fit.med$states, med$states)
plot(fit.med, legend = FALSE)

## ----bmwd-nstart--------------------------------------------------------------
fit.best <- thmm(med$data, nstates = 5, type = "bmwd", nstart = 5, tol = 1e-8)
c(single = fit.med$loglik, best.of.5 = fit.best$loglik)

## ----ou-sim-------------------------------------------------------------------
set.seed(137)
ou.par <- cbind(mean = c(-2, 0, 4, 2, 1), rate = c(4, 4, 8, 2, 20))
ou <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A1,
            type = "ou", par = ou.par, len = 100)
matplot(t(ou$data), type = "l", lty = 1, col = ou$states + 1,
        xlab = "grid point", ylab = "", main = "Ornstein-Uhlenbeck")

## ----ou-fit-------------------------------------------------------------------
fit.ou <- thmm(ou$data, nstates = 5, type = "ou", nstart = 3, tol = 1e-8)
round(fit.ou$par, 2)
table(estimated = fit.ou$states, truth = ou$states)
ari(fit.ou$states, ou$states)

## ----fbm----------------------------------------------------------------------
set.seed(137)
fbm <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A1, type = "bmwd",
             par = c(-10, -6, -2, 0, 2), len = 100, hurst = 0.8)
res <- sapply(c(0.25, 0.5, 0.8), function(h) {
  set.seed(1)
  f <- thmm(fbm$data, 5, type = "bmwd", hurst = h, nstart = 3, tol = 1e-8)
  c(hurst = h, ari = ari(f$states, fbm$states), sigma = f$sigma,
    sort(f$par))
})
round(t(res), 2)

## ----nonpar-sim---------------------------------------------------------------
set.seed(137)
tt <- seq_len(100) / 100
mu <- rbind(sin(2 * pi * tt), sin(2 * pi * (tt + 0.2)), sin(2 * pi * (tt + 0.4)),
            sin(2 * pi * (tt + 0.6)), sin(2 * pi * (tt + 0.8)))
sinu <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A1,
              type = "nonpar", par = mu, sigma = 0.4)

## ----nonpar-fit---------------------------------------------------------------
fit.l2 <- thmm(sinu$data, 5, type = "nonpar", norm = "L2",
               nstart = 3)
ari(fit.l2$states, sinu$states)
fit.w21 <- thmm(sinu$data, 5, type = "nonpar", norm = "W21",
                nstart = 3)
ari(fit.w21$states, sinu$states)
plot(fit.w21, legend = FALSE)

## ----nonpar-bm----------------------------------------------------------------
fit.np <- thmm(med$data, 5, type = "nonpar", norm = "L2", start = "kmeans", nstart = 3)
ari(fit.np$states, med$states)

## ----sigma--------------------------------------------------------------------
c(bmwd = fit.med$sigma, ou = fit.ou$sigma)

## ----sigma-l2-----------------------------------------------------------------
sapply(c(0.01, 0.05, 0.2), function(s) {
  set.seed(1)
  f <- thmm(sinu$data, 5, type = "nonpar", norm = "L2", sigma = s)
  c(sigma = s, ari = ari(f$states, sinu$states), loglik = f$loglik)
})

## ----predict------------------------------------------------------------------
set.seed(2)
new <- rthmm(50, init.prob = c(1, 0, 0, 0, 0), trans = A1,
             type = "bmwd", par = c(-8, -4, 0, 4, 8), len = 100)
s <- predict(fit.med, new$data)
ari(s, new$states)
head(round(predict(fit.med, new$data, type = "posterior"), 3))

## ----multi--------------------------------------------------------------------
set.seed(4)
drift2 <- rbind(c(-3, 3), c(3, 3), c(0, -3))   # three states in two dimensions
A3 <- matrix(0.1, 3, 3) + 0.6 * diag(3)
bm2 <- rthmm(150, init.prob = c(1, 0, 0), trans = A3, type = "bmwd",
             par = drift2, len = 100)
dim(bm2$data)
fit2 <- thmm(bm2$data, 3, type = "bmwd")
round(fit2$par, 2)
ari(fit2$states, bm2$states)