Skip to contents

Forward algorithm to calculate the HMM log-likelihood

forward.Rd

Calculates the log-likelihood of a sequence of observations under a hidden Markov model using the forward algorithm (Zucchini, MacDonald & Langrock, 2016).

Usage

forward(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  logspace = FALSE,
  bw = NULL,
  report = TRUE,
  ad = NULL,
  tod = NULL
)

Arguments

delta

initial distribution; either

  • a vector of length nStates, or

  • a matrix of dimension c(nTracks, nStates), if trackID is provided, or

  • "stationary", for the (periodically) stationary distribution of Gamma at the start of each track (see ‘Details’).

Gamma

transition probability matrix; either

  • a matrix of dimension c(nStates, nStates),

  • an array of dimension c(nStates, nStates, nTracks) if trackID is provided,

  • an array of dimension c(nStates, nStates, nObs) for time-varying transition probabilities, or

  • an array of dimension c(nStates, nStates, L) for periodic transition probabilities, if tod is provided.

allprobs

matrix of state-dependent probabilities or density values of dimension c(nObs, nStates)

trackID

optional vector of length nObs containing nTracks unique IDs that separate tracks (see ‘Details’). The observations of each track need to be consecutive.

logspace

logical; if TRUE, allprobs is assumed to be on the log-scale, improving numerical stability for small probabilities.

bw

optional positive integer specifying the bandwidth for a banded approximation of the forward algorithm, inducing a banded Hessian w.r.t. the observations. Defaults to NULL (exact algorithm). Approximation error decays geometrically in bw. Only used with automatic differentiation; otherwise the exact likelihood is computed.

report

logical; if TRUE (default), delta, Gamma, allprobs, and trackID are reported from the fitted model. Requires ad = TRUE.

ad

logical; whether to use automatic differentiation. Determined automatically — for debugging only.

tod

optional integer vector of length nObs with values in 1, ..., L, indexing the slices of Gamma for periodic models (e.g. time of day with L = 24 for hourly data). Gamma[,,tod[t]] is the transition into time t, hence the periodically stationary initial distribution is stationary_p(Gamma, tod[1] %% L + 1).

Value

HMM log-likelihood for given data and parameters

Details

The likelihood of observations \(x_1, \dots, x_T\) is $$\mathcal{L}(\boldsymbol{\theta}) = \boldsymbol{\delta}^{(1)} \mathbf{P}(x_1) \boldsymbol{\Gamma}^{(2)} \mathbf{P}(x_2) \boldsymbol{\Gamma}^{(3)} \cdots \boldsymbol{\Gamma}^{(T)} \mathbf{P}(x_T) \mathbf{1}^\top,$$ where \(\mathbf{P}(x_t)\) is a diagonal matrix with the \(t\)-th row of allprobs on its diagonal, \(\boldsymbol{\Gamma}^{(t)}\) is the transition probability matrix into time \(t\), and \(\mathbf{1}\) is a row vector of ones. The forward algorithm evaluates this product recursively and on the log-scale.

If trackID is provided, the total log-likelihood is the sum of each track's likelihood contribution. In this case, Gamma can be

  • a matrix (same transition probabilities for each track),

  • an array of dimension c(nStates, nStates, nTracks) (track-specific transition probability matrix), or

  • an array of dimension c(nStates, nStates, nObs) for time-varying transition probabilities, or

  • an array of dimension c(nStates, nStates, L) for periodic transition probabilities, if tod is provided.

Additionally, delta can be a vector (same initial distribution for each track) or a matrix of dimension c(nTracks, nStates) (track-specific initial distributions).

With delta = "stationary", the initial distribution of each track is computed from Gamma:

  • the stationary distribution of Gamma (or of its slice for each track, if track-specific),

  • the periodically stationary distribution at the first time point of each track, if tod is provided (see stationary_p),

  • the stationary distribution of the first slice of each track, for time-varying transition probabilities, i.e. at the covariate values of the first observation.

This requires Gamma to contain transition probability matrices (not, e.g., sub-stochastic matrices), which is checked without automatic differentiation. The computed delta is reported as a matrix of dimension c(nTracks, nStates).

Note: forward should only be called once in the likelihood function, with a trackID argument for multiple tracks. Downstream functions like lama, viterbi, stateprobs or pseudo_res rely on the quantities it reports, and only those of the last call are kept.

References

Zucchini, W., MacDonald, I.L., & Langrock, R. (2016). Hidden Markov Models for Time Series: An Introduction Using R (2nd ed.). Chapman & Hall/CRC.

Mews, Koslik, & Langrock (2025). How to build your latent Markov model: The role of time and space. Statistical Modelling, 25(6), 481–507.

See also

Other forward algorithms: forward_g(), forward_hsmm(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

## negative log likelihood function
nll = function(par, step) {
 # parameter transformations for unconstrained optimisation
 Gamma = tpm(par[1:2]) # multinomial logit link
 delta = stationary(Gamma) # stationary HMM
 mu = exp(par[3:4])
 sigma = exp(par[5:6])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward(delta, Gamma, allprobs)
}

## fitting an HMM to the trex data
par = c(-2,-2,            # initial tpm params (logit-scale)
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:1000])