--- title: "Hidden Markov models for functional data with funHMM" author: "Adam B Kashlak" output: rmarkdown::html_vignette: toc: true fig_width: 7 fig_height: 4.5 vignette: > %\VignetteIndexEntry{Hidden Markov models for functional data with funHMM} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(funHMM) ``` ## Introduction A hidden Markov model (HMM) describes a sequence of observations $O_1, \dots, O_n$ that are conditionally independent given a hidden state sequence $s_1, \dots, s_n$ that evolves as a Markov chain on $\{1, \dots, p\}$ with initial probabilities $\eta_j$ and transition matrix $A = (a_{ij})$. Each state $j$ has an *emission function* $b_j$ that scores how plausible an observation is under state $j$; in the classic HMM for multivariate data, $b_j$ is a Gaussian density. The **topological hidden Markov model** (THMM) of Kashlak, Loliencar and Heo (2023) lets the observations be *curves*: sample paths of a stochastic process, or smooth functional data. There is no Lebesgue measure in infinite dimensions and hence no density to use for $b_j$. Instead, each state is a Gaussian measure on the function space, shifted to a state-specific centre $h_j$ in its Cameron-Martin space $H$, and $b_j$ is replaced by the **Onsager-Machlup functional**, the limit of the ratio of small-ball probabilities $$ \log b_j(O_t) = \lim_{\varepsilon\to 0} \log \frac{P(\|Y - O_t\| < \varepsilon)}{P(\|W\| < \varepsilon)} = -\tfrac{1}{2\sigma^2}\, |O_t - h_j|_H^2 \quad (\text{up to terms not depending on } h_j). $$ Everything else in the Baum-Welch (EM) and Viterbi algorithms carries over unchanged, and the re-estimated centres are posterior-weighted averages of the observations. The `funHMM` package implements three emission models with the forward, backward, re-estimation and Viterbi steps written in C: | `type` | model | state parameters | |------------|-------------------------------------------------------|------------------| | `"bmwd"` | Brownian motion with drift $dY = c_j\,d\tau + \sigma\,dW$ (fractional Brownian motion via `hurst`) | drift $c_j$ | | `"ou"` | Ornstein-Uhlenbeck $dY = \theta_j(\mu_j - Y)\,d\tau + \sigma\,dW$ | mean $\mu_j$, rate $\theta_j$ | | `"nonpar"` | non-parametric mean curve $h_j$ under the `L2`, `W21` or `W22` norm | curve $h_j$ | Throughout, curves are stored **one per row** of a matrix, observed on `len` equally spaced points of $[0, 1]$, and the rows are in time order. The examples below reproduce the simulation studies in Section 5 of the paper. ## Brownian motion with drift We simulate 200 Brownian sample paths with five states whose drifts are $-4, -2, 0, 2, 4$ ("low separation"). The transition matrix $A_1$ of the paper keeps the chain in its current state with probability $0.64$. ```{r 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") ``` `thmm()` fits the model. With `type = "bmwd"` the only per-curve statistic that matters is the increment $O_t(1) - O_t(0)$, so the fit is essentially instantaneous. ```{r bmwd-fit} fit.low <- thmm(low$data, nstates = 5, type = "bmwd", tol = 1e-8) fit.low ``` The decoded states are compared with the truth through a confusion table and the adjusted Rand index (`ari()`, 1 is perfect agreement): ```{r bmwd-eval} table(estimated = fit.low$states, truth = low$states) ari(fit.low$states, low$states) ``` With larger separation between the drifts the states are almost perfectly recovered: ```{r 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) ``` `plot()` colours the curves by decoded state and overlays the fitted drift lines. ### Starting values and multiple starts As with any EM algorithm the result depends on the starting values. By default (`start = "kmeans"`) the per-curve statistics are clustered with k-means and one re-estimation step from that hard clustering provides the initial parameters; `start = "random"` instead picks `nstates` curves at random as the initial centres. `nstart` runs the whole algorithm several times and keeps the fit with the largest final log-likelihood: ```{r 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) ``` ## Ornstein-Uhlenbeck process The OU process reverts towards a state-specific mean $\mu_j$ at rate $\theta_j$. The paper uses five states with means $(-2, 0, 4, 2, 1)$ and rates $(4, 4, 8, 2, 20)$. ```{r 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") ``` The Onsager-Machlup functional of the OU process is quadratic in $(\theta\mu, \theta)$, so `thmm` re-estimates both parameters in closed form (the paper used a numerical optimiser for this step). ```{r 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) ``` Both the means and the rates are recovered well. Note that the rate is constrained to be non-negative; a fitted rate of exactly zero means the state behaves like Brownian motion with drift (and its `mean` is reported as `NA`). ## Fractional Brownian motion For `type = "bmwd"` the argument `hurst` selects the Hurst parameter of the driving fractional Brownian motion. Values above $1/2$ give smoother paths with positively correlated increments, values below $1/2$ rougher paths. Following the paper we simulate with Hurst parameter $0.8$ and drifts $(-10, -6, -2, 0, 2)$ and then fit the model with three different values of `hurst`: ```{r 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) ``` The Hurst parameter enters both the drift statistic (a weighted integral of the increments) and the automatic estimate of `sigma`, so the correctly specified fit also reports a diffusion coefficient close to one. ## Non-parametric mean curves When no parametric form is assumed, `type = "nonpar"` estimates a mean curve for each state. The norm used to compare a curve with a mean is chosen with `norm`: * `"L2"`: $\int (O - h)^2\,d\tau$; * `"W21"`: $\int (\dot O - \dot h)^2\,d\tau$, the Cameron-Martin norm of the standard Wiener measure; * `"W22"`: $\int (\ddot O - \ddot h)^2\,d\tau$. The paper's example uses phase-shifted sinusoids corrupted by smooth Brownian-bridge noise, which `rbridge()` generates: ```{r 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) ``` The default k-means start matters most for the non-parametric model, where random starts frequently converge to local optima in which two states share a mean curve: ```{r 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) ``` The fitted mean curves (black) sit on top of the true sinusoids. The non-parametric model can also be applied to the Brownian data of the first section; it loses only a little accuracy for not knowing the parametric form: ```{r nonpar-bm} fit.np <- thmm(med$data, 5, type = "nonpar", norm = "L2", start = "kmeans", nstart = 3) ari(fit.np$states, med$states) ``` ## The scale parameter `sigma` Every log-emission function is divided by $\sigma^2$, the variance (diffusion coefficient) of the driving Gaussian measure, which must be common to all states. It therefore acts as a *temperature*: a small `sigma` makes the emissions dominate the transition probabilities and the fit behaves like a time-aware k-means clustering, while a large `sigma` lets the Markov chain smooth over the emissions. By default `sigma` is estimated from the realised quadratic variation of the curves (second differences, so drifts and smooth means do not affect it). For data generated from the standard models this recovers the true coefficient: ```{r sigma} c(bmwd = fit.med$sigma, ou = fit.ou$sigma) ``` For smooth curves under the `L2` norm the estimate is tiny and the fit is almost a hard clustering. Supplying `sigma` explicitly is the way to trade off emissions against transitions: ```{r 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) }) ``` ## Decoding new sequences `predict()` runs the forward-backward and Viterbi algorithms on a new sequence of curves with the fitted parameters held fixed: ```{r 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)) ``` `simulate()` draws a new sequence from a fitted model, and `logLik()` returns the final log-likelihood with the number of free parameters as its `df` attribute. ## Multivariate curves All models except `"ou"` accept `d`-dimensional curves supplied as an `n x len x d` array. For `"bmwd"` each coordinate is an independent Brownian motion with its own drift, and for `"nonpar"` the squared norms are summed over coordinates. ```{r 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) ``` ## Practical notes * **Preprocessing.** The `"bmwd"`, `"ou"` and `"W21"` models differentiate the curves, so noisy discrete data should be smoothed beforehand (the paper uses kernel smoothing for EEG spectra and snowfall curves). The `"L2"` norm is the most forgiving choice for rough data. * **Grid.** Curves are assumed to be observed on the grid $k/\mathrm{len}$, $k = 1, \dots, \mathrm{len}$; only the spacing matters for the fits. * **Convergence.** The log-likelihood built from Onsager-Machlup functionals is not a true likelihood and can be positive; the EM iterations still increase it monotonically (Theorem 2 of the paper). Iterations stop when its relative change falls below `tol`. * **Local optima.** Keep the default `start = "kmeans"` and use `nstart > 1` for difficult problems. The paper fits each model 20 times and keeps the best. With `start = "random"` the Ornstein-Uhlenbeck model in particular often ends in a poor local optimum. ## References Kashlak, A. B., Loliencar, P. and Heo, G. (2023). Topological Hidden Markov Models. *Journal of Machine Learning Research*, 24(340), 1-49. Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. *Proceedings of the IEEE*, 77(2), 257-286.