Hidden Markov models for functional data with 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:

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\).

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.

fit.low <- thmm(low$data, nstates = 5, type = "bmwd", tol = 1e-8)
fit.low
#> Topological hidden Markov model
#> Call: thmm(data = low$data, nstates = 5, type = "bmwd", tol = 1e-08)
#> Model: Brownian motion with drift; 5 states; 200 curves on 100 grid points
#> sigma = 0.9987 (estimated); 219 iterations (converged); log-likelihood -258.3
#> 
#> Initial state probabilities:
#> [1] 0 1 0 0 0
#> 
#> Transition matrix:
#>        [,1]   [,2]   [,3]   [,4]   [,5]
#> [1,] 0.7623 0.1700 0.0000 0.0000 0.0677
#> [2,] 0.0000 0.5501 0.2468 0.2031 0.0000
#> [3,] 0.0000 0.0000 0.5788 0.2066 0.2146
#> [4,] 0.0396 0.0321 0.0934 0.6035 0.2314
#> [5,] 0.1433 0.0440 0.1106 0.0000 0.7021
#> 
#> State parameters:
#>  drift1  drift2  drift3  drift4  drift5 
#>  2.7840 -4.1842  0.4828  4.1921 -2.0432 
#> 
#> State counts (Viterbi):
#> 
#>  1  2  3  4  5 
#> 37 25 39 35 64

The decoded states are compared with the truth through a confusion table and the adjusted Rand index (ari(), 1 is perfect agreement):

table(estimated = fit.low$states, truth = low$states)
#>          truth
#> estimated  1  2  3  4  5
#>         1  0  0  2 21 14
#>         2 22  3  0  0  0
#>         3  0  1 25 13  0
#>         4  0  0  0  3 32
#>         5 12 48  2  2  0
ari(fit.low$states, low$states)
#> [1] 0.5213842

With larger separation between the drifts the states are almost perfectly recovered:

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)
#> [1] -7.8618570 -3.8590856  0.3130001  4.2131815  8.0365126
ari(fit.med$states, med$states)
#> [1] 0.9572509
plot(fit.med, legend = FALSE)

plot() colours the curves by decoded state and overlays the fitted drift lines.

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)
#>    single best.of.5 
#> -314.1303 -314.1303

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)\).

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).

fit.ou <- thmm(ou$data, nstates = 5, type = "ou", nstart = 3, tol = 1e-8)
round(fit.ou$par, 2)
#>       mean  rate
#> [1,] -2.05  3.55
#> [2,]  4.07  6.89
#> [3,]  2.05  1.96
#> [4,]  1.01 16.46
#> [5,]  0.03  4.39
table(estimated = fit.ou$states, truth = ou$states)
#>          truth
#> estimated  1  2  3  4  5
#>         1 34  0  0  0  0
#>         2  0  0 29  0  0
#>         3  0  0  0 38  0
#>         4  0  0  0  0 46
#>         5  0 52  0  1  0
ari(fit.ou$states, ou$states)
#> [1] 0.986116

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:

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)
#>      hurst  ari sigma                             
#> [1,]  0.25 0.70  0.05 -11.56 -6.91 -2.21 0.72 3.05
#> [2,]  0.50 0.71  0.17  -9.81 -5.86 -1.85 0.68 2.65
#> [3,]  0.80 0.71  1.00  -9.32 -5.51 -1.50 1.06 2.48

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:

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:

fit.l2 <- thmm(sinu$data, 5, type = "nonpar", norm = "L2",
               nstart = 3)
ari(fit.l2$states, sinu$states)
#> [1] 1
fit.w21 <- thmm(sinu$data, 5, type = "nonpar", norm = "W21",
                nstart = 3)
ari(fit.w21$states, sinu$states)
#> [1] 1
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:

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

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:

c(bmwd = fit.med$sigma, ou = fit.ou$sigma)
#>      bmwd        ou 
#> 0.9987456 1.0403373

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:

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)
})
#>             [,1]      [,2]      [,3]
#> sigma       0.01     0.050    0.2000
#> ari         1.00     1.000    1.0000
#> loglik -25036.91 -1216.987 -286.3981

Decoding new sequences

predict() runs the forward-backward and Viterbi algorithms on a new sequence of curves with the fitted parameters held fixed:

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)
#> [1] 0.9066035
head(round(predict(fit.med, new$data, type = "posterior"), 3))
#>       [,1]  [,2] [,3]  [,4]  [,5]
#> [1,] 0.000 0.000    0 0.000 1.000
#> [2,] 0.000 0.001    0 0.999 0.000
#> [3,] 0.000 0.005    0 0.995 0.000
#> [4,] 0.000 0.000    0 1.000 0.000
#> [5,] 0.011 0.000    0 0.000 0.989
#> [6,] 0.989 0.000    0 0.011 0.000

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.

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)
#> [1] 150 100   2
fit2 <- thmm(bm2$data, 3, type = "bmwd")
round(fit2$par, 2)
#>       [,1]  [,2]
#> [1,] -3.00  2.94
#> [2,] -0.05 -2.97
#> [3,]  2.83  3.01
ari(fit2$states, bm2$states)
#> [1] 1

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. https://jmlr.org/papers/v24/22-0685.html

Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2), 257-286.