Package {funHMM}

funHMM: Hidden Markov Models for Functional Data

Description

Fits hidden Markov models to sequences of observations that live in a function space, such as sample paths of a stochastic process or smooth curves. The emission "densities" are Onsager-Machlup functionals of a Gaussian measure on a locally convex topological vector space, so no finite-dimensional projection of the curves is needed. The Baum-Welch (EM) and Viterbi algorithms are implemented in C.

Details

The main function is thmm(). Simulated data can be produced with rthmm(), rbm(), rou() and rbridge(). See vignette("funHMM") for a worked introduction that reproduces the simulations of the reference below.

Author(s)

Maintainer: Adam B Kashlak kashlak@ualberta.ca

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.

Adjusted Rand index

Description

Measures the agreement between two labellings of the same objects, e.g. the decoded states of a thmm() fit and the true states of a simulation. A value of one indicates identical partitions (up to relabelling) and values near zero indicate agreement no better than chance.

Usage

ari(x, y)

Arguments

Value

A single number.

References

Hubert, L. and Arabie, P. (1985). Comparing partitions. Journal of Classification, 2, 193–218.

Examples

ari(c(1, 1, 2, 2, 3), c(2, 2, 1, 1, 3))   # 1: identical up to relabelling
ari(c(1, 1, 2, 2), c(1, 2, 1, 2))         # negative: worse than chance

Plot a fitted THMM

Description

Plots the observed curves coloured by their decoded (Viterbi) state. For the non-parametric model the fitted mean curves are overlaid; for Brownian motion with drift the fitted linear drifts are overlaid.

Usage

## S3 method for class 'thmm'
plot(
  x,
  which = seq_len(dim(x$data)[3L]),
  col = NULL,
  lwd.mean = 3,
  legend = TRUE,
  ...
)

Arguments

Value

x, invisibly.

Examples

set.seed(3)
tt <- seq(0, 1, length.out = 40)
sim <- rthmm(60, init.prob = c(1, 0), trans = matrix(c(.8, .2, .2, .8), 2),
             type = "nonpar", par = rbind(tt, 1 - tt), sigma = 0.3)
fit <- thmm(sim$data, 2, type = "nonpar", norm = "L2", start = "kmeans")
plot(fit)

Decode new curves with a fitted THMM

Description

Runs the forward-backward and Viterbi algorithms on a new sequence of curves using the parameters of a fitted model.

Usage

## S3 method for class 'thmm'
predict(object, newdata = NULL, type = c("states", "posterior", "loglik"), ...)

Arguments

Value

An integer vector of states, a matrix of posterior probabilities, or a single number.

Examples

set.seed(2)
A <- matrix(c(.9, .1, .1, .9), 2)
sim <- rthmm(80, init.prob = c(1, 0), trans = A, type = "bmwd",
             par = c(-3, 3), len = 50)
fit <- thmm(sim$data, 2, type = "bmwd")
new <- rthmm(40, init.prob = c(1, 0), trans = A, type = "bmwd",
             par = c(-3, 3), len = 50)
table(predict(fit, new$data), new$states)

Simulate sample paths of (fractional) Brownian motion with drift

Description

Simulates n sample paths of Y_\tau = c\tau + \sigma W^H_\tau on the grid \tau_k = k/\code{len}, k = 1, \dots, \code{len}, where W^H is fractional Brownian motion with Hurst parameter hurst (standard Brownian motion for hurst = 0.5).

Usage

rbm(n, len = 100L, drift = 0, hurst = 0.5, sigma = 1)

Arguments

Details

For hurst != 0.5 the paths are generated from the covariance function (\tau_1^{2H} + \tau_2^{2H} - |\tau_1 - \tau_2|^{2H})/2 using a Cholesky factor of the ⁠len x len⁠ covariance matrix, so large values of len are slow.

Value

An ⁠n x len⁠ matrix with one path per row.

Examples

x <- rbm(5, len = 100, drift = c(-2, -1, 0, 1, 2))
matplot(t(x), type = "l", lty = 1)
y <- rbm(3, len = 100, hurst = 0.8)

Simulate smooth Brownian bridge noise

Description

Simulates smooth random curves from a truncated Karhunen-Loeve expansion of the Brownian bridge, \epsilon(\tau) = \sqrt{2}\sum_{k=1}^{K} Z_k \sin(k\pi\tau)/(k\pi) with Z_k independent N(0, \sigma^2). This is the error process used in the non-parametric simulations of Kashlak, Loliencar and Heo (2023).

Usage

rbridge(n, len = 100L, sigma = 1, nterms = 16L)

Arguments

Value

An ⁠n x len⁠ matrix with one curve per row.

Examples

e <- rbridge(5, len = 100, sigma = 0.4)
matplot(t(e), type = "l", lty = 1)

Simulate a Markov chain

Description

Simulate a Markov chain

Usage

rmarkov(n, init.prob, trans)

Arguments

Value

An integer vector of states in 1:length(init.prob).

Examples

rmarkov(20, c(1, 0), matrix(c(.9, .1, .1, .9), 2))

Simulate sample paths of the Ornstein-Uhlenbeck process

Description

Simulates dY = \theta(\mu - Y)\,d\tau + \sigma\,dW, Y_0 = 0, by an Euler scheme on the grid \tau_k = k/\code{len}.

Usage

rou(n, len = 100L, mean = 0, rate = 1, sigma = 1)

Arguments

Value

An ⁠n x len⁠ matrix with one path per row.

Examples

x <- rou(4, len = 100, mean = c(-2, 0, 2, 4), rate = c(2, 4, 8, 20))
matplot(t(x), type = "l", lty = 1)

Simulate data from a topological hidden Markov model

Description

Generates a hidden state sequence from a Markov chain and, conditional on the states, a sequence of curves from one of the emission models of thmm().

Usage

rthmm(
  n,
  init.prob,
  trans,
  type = c("bmwd", "ou", "nonpar"),
  par,
  len = 100L,
  hurst = 0.5,
  sigma = 1,
  nterms = 16L
)

Arguments

Value

A list with components data (an ⁠n x len⁠ matrix, or an ⁠n x len x d⁠ array) and states (integer vector of true states).

Examples

A <- matrix(0.09, 5, 5) + 0.55 * diag(5)   # matrix A1 of the paper
sim <- rthmm(200, init.prob = c(1, 0, 0, 0, 0), trans = A,
             type = "bmwd", par = c(-4, -2, 0, 2, 4), len = 100)
matplot(t(sim$data), type = "l", lty = 1, col = sim$states + 1)

Simulate from a fitted THMM

Description

Simulate from a fitted THMM

Usage

## S3 method for class 'thmm'
simulate(object, nsim = 1L, seed = NULL, ...)

Arguments

Value

A list with components data and states, see rthmm().

Examples

set.seed(4)
sim <- rthmm(60, init.prob = c(1, 0), trans = matrix(c(.8, .2, .2, .8), 2),
             type = "bmwd", par = c(-3, 3), len = 30)
fit <- thmm(sim$data, 2, type = "bmwd")
new <- simulate(fit, nsim = 10)
dim(new$data)

Fit a topological hidden Markov model

Description

Fits a hidden Markov model to a time-ordered sequence of curves (functional observations) using the Baum-Welch algorithm, with emission functions given by Onsager-Machlup functionals as described in Kashlak, Loliencar and Heo (2023). The most likely state sequence is then decoded with the Viterbi algorithm. The heavy lifting is done in C.

Usage

thmm(
  data,
  nstates,
  type = c("bmwd", "ou", "nonpar"),
  norm = c("W21", "L2", "W22"),
  hurst = 0.5,
  sigma = NULL,
  par = NULL,
  trans = NULL,
  init.prob = NULL,
  start = c("kmeans", "random"),
  nstart = 1L,
  max.iter = 500L,
  min.iter = 10L,
  tol = 1e-06,
  verbose = FALSE
)

Arguments

Details

Let O_1, \dots, O_n be the observed curves and let b_j(O_t) denote the emission function of state j. Writing \dot O_t for the derivative of the curve, the models are:

"bmwd"

Brownian motion with drift, dY = c_j\,d\tau + \sigma\,dW. Up to a term that does not depend on c_j, \log b_j(O_t) = -V (D_t - c_j)^2 / (2\sigma^2) where D_t = \Gamma(2\upsilon+2)/\Gamma(\upsilon+1) \int_0^1 \tau^\upsilon \dot O_t(\tau)\,d\tau, \upsilon = 1/2 - \code{hurst} and V = \Gamma(\upsilon+1)^2/\{\Gamma(2\upsilon+1)\Gamma(2\upsilon+2)\}. For hurst = 0.5 this reduces to D_t = O_t(1) - O_t(0) and V = 1. The drift is re-estimated as the posterior-weighted mean of D_t. For d-dimensional curves the coordinates are treated as independent Brownian motions with a d-vector of drifts per state.

"ou"

Ornstein-Uhlenbeck process dY = \theta_j(\mu_j - Y)\,d\tau + \sigma\,dW, parametrised internally by b_0 = \theta_j\mu_j and b_1 = \theta_j. Then \log b_j(O_t) = \{b_0 A_t - b_1 E_t - b_0^2/2 + b_0 b_1 B_t - b_1^2 C_t/2\}/\sigma^2 + b_1/2 with A_t = O_t(1)-O_t(0), B_t = \int O_t, C_t = \int O_t^2 and E_t = \int O_t \circ dO_t (Stratonovich). Because this is quadratic in (b_0, b_1) the re-estimation step has a closed form, with b_1 constrained to be non-negative.

"nonpar"

Non-parametric mean curves h_j with \log b_j(O_t) = -|O_t - h_j|_H^2 / (2\sigma^2) where |\cdot|_H is the chosen Cameron-Martin norm, approximated by a Riemann sum on the observation grid: \int (O-h)^2 for "L2", \int (\dot O - \dot h)^2 for "W21" and \int (\ddot O - \ddot h)^2 for "W22". The mean curves are re-estimated as posterior-weighted averages of the observed curves.

All forward and backward probabilities are computed on the log scale. Iteration stops when the relative change in log-likelihood drops below tol (after at least min.iter iterations) or after max.iter iterations. Note that the "log-likelihood" is built from Onsager-Machlup functionals rather than densities and may be positive.

Automatic choice of sigma. When sigma = NULL the diffusion coefficient is estimated from second differences of the curves (the realised quadratic variation, which does not depend on the drift or mean curves). For "bmwd" and "ou" data generated from the standard model this returns a value close to one; for the non-parametric model it returns the scale of the noise in the chosen norm. For smooth curves under the "L2" norm the estimate can be very small, in which case the emissions dominate the transition probabilities and the fit behaves like a k-means clustering that respects time ordering. Supplying a larger sigma smooths the state assignments.

Value

An object of class "thmm": a list with components

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

See Also

rthmm() to simulate data, predict.thmm() to decode new sequences, plot.thmm().

Examples

set.seed(1)
A <- matrix(0.1, 3, 3); diag(A) <- 0.8
sim <- rthmm(100, init.prob = c(1, 0, 0), trans = A,
             type = "bmwd", par = c(-4, 0, 4), len = 50)
fit <- thmm(sim$data, nstates = 3, type = "bmwd")
fit
table(fit$states, sim$states)
ari(fit$states, sim$states)

## Ornstein-Uhlenbeck curves with two states
sim <- rthmm(100, init.prob = c(1, 0), trans = matrix(c(.9, .1, .1, .9), 2),
             type = "ou", par = cbind(mean = c(-2, 2), rate = c(4, 4)),
             len = 50)
fit <- thmm(sim$data, nstates = 2, type = "ou")
fit$par

## non-parametric mean curves
tt <- seq(0, 1, length.out = 50)
mu <- rbind(sin(2 * pi * tt), cos(2 * pi * tt))
sim <- rthmm(100, init.prob = c(1, 0), trans = matrix(c(.9, .1, .1, .9), 2),
             type = "nonpar", par = mu, sigma = 0.4)
fit <- thmm(sim$data, nstates = 2, type = "nonpar", norm = "L2")
ari(fit$states, sim$states)
plot(fit)