Gaussian HMM

Hidden Markov model with Gaussian emissions, fitted by expectation-maximisation.

Description

Each hidden state $k$ emits a Gaussian $\mathcal{N}(\mu_k, \operatorname{diag}(\sigma_k^2))$, and the states follow a first-order Markov chain with initial probabilities $\pi$ and transition matrix $A$. All parameters are estimated by maximum likelihood with the Baum-Welch algorithm, an expectation-maximisation procedure whose E-step is the forward-backward algorithm (Rabiner, 1989). Each time step then receives the state of the most probable path (Viterbi).

Initialisation: k-means on the (z-normalised) observations gives the means and variances; $\pi$ and $A$ start uniform. n_init initialisations with different seeds are run and the one with the highest log-likelihood is kept.
Regularisation: reg_covar (positive) is added to every variance at each M-step, which limits the collapse of a state onto a few points. On integer or coarsely quantised data a state can still collapse onto a single value; raise reg_covar then.
Number of states: given by n_states, or, with n_states=None, chosen by the Bayesian information criterion among 1, ..., max_states. The fitted model has $(K-1) + K(K-1) + 2Kd$ free parameters.

The emissions are independent given the state, so the model has no notion of duration or of temporal shape within a state. It suits regimes that differ in level or spread, not regimes that differ only in frequency or waveform: on oscillating signals the states follow the oscillation and the labels change within a regime. For the same reason BIC often selects max_states on long autocorrelated series (a warning says so), and the number of states is then better given.

Type: state detection
Supervision: unsupervised (BIC) or semi-supervised (n_states)
Scope: univariate and multivariate
Complexity: \(O(nK^2 + nKd)\) per EM iteration; the BIC search
multiplies it by the number of candidates
Requires: nothing beyond the core dependencies; numba
(tsseg[accelerators]) speeds up the recursions considerably

Parameters

Name	Type	Default	Description
`n_states`	int / None	`None`	Number of hidden states; `None` selects it by BIC.
`max_states`	int	`8`	Largest number of states tried by BIC.
`n_iter`	int	`100`	Maximum number of EM iterations per initialisation.
`tol`	float	`1e-4`	EM stops when the log-likelihood per time step improves by less.
`n_init`	int	`3`	Number of k-means initialisations.
`reg_covar`	float	`1e-3`	Added to the variances at each M-step (positive; units of the normalised data).
`normalize`	bool	`True`	Z-normalise each channel before fitting.
`random_state`	int / None	`None`	Seed of the k-means initialisations.
`axis`	int	`0`	Time axis.

Usage

from tsseg.algorithms import GaussianHMMDetector

# number of states chosen by BIC
labels = GaussianHMMDetector(random_state=0).fit_predict(X)

# known number of states
labels = GaussianHMMDetector(n_states=4, random_state=0).fit_predict(X)

Implementation: written for tsseg in NumPy, with numba for the forward-backward and Viterbi recursions. With the same starting point, the parameters and log-likelihood after each EM iteration match those of hmmlearn (GaussianHMM with covariance_type="diag" and covars_prior=0) to $10^{-10}$ or better. It replaces HMM, which only decoded with user-supplied parameters.

Reference: Rabiner (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE 77(2); Schwarz (1978), Estimating the dimension of a model, The Annals of Statistics 6(2).

API reference

Gaussian hidden Markov model fitted by Baum-Welch (EM) and decoded by Viterbi.

class tsseg.algorithms.gaussian_hmm.detector.GaussianHMMDetector(*, n_states=None, max_states=8, n_iter=100, tol=0.0001, n_init=3, reg_covar=0.001, normalize=True, random_state=None, axis=0)[source]

Bases: BaseSegmenter

State detection with a Gaussian hidden Markov model.

Each hidden state emits a Gaussian with diagonal covariance. The parameters (initial and transition probabilities, means, variances) are estimated by maximum likelihood with the Baum-Welch algorithm, an expectation-maximisation procedure whose E-step is the forward-backward algorithm [1]. Each time step then gets the state of the most probable path (Viterbi).

With n_states=None the number of states is chosen by the Bayesian information criterion [2] among 1, ..., max_states.

Parameters:

n_states (int | None) – Number of hidden states. None selects it by BIC.
max_states (int) – Largest number of states tried when n_states is None.
n_iter (int) – Maximum number of EM iterations per initialisation.
tol (float) – EM stops when the log-likelihood per time step improves by less than tol.
n_init (int) – Number of initialisations (k-means with different seeds); the fit with the highest log-likelihood is kept.
reg_covar (float) – Added to every variance at each M-step (must be positive), which limits the collapse of a state onto a few points. Expressed in units of the (normalised) data.
normalize (bool) – Z-normalise each channel before fitting.
random_state (int | None) – Seed of the k-means initialisations. Pass an int for reproducible results.
axis (int) – Time axis.

n_states_

Number of states $K$ of the fitted model.

Type:: int

startprob_

Initial state probabilities.

Type:: ndarray of shape (K,)

transmat_

transmat_[i, j] is the probability of moving from state i to j.

Type:: ndarray of shape (K, K)

means_

Emission means.

Type:: ndarray of shape (K, n_channels)

variances_

Emission variances.

Type:: ndarray of shape (K, n_channels)

loglik_

Log-likelihood of the fitted model.

Type:: float

bic_

BIC of each candidate number of states (only when n_states=None).

Type:: dict

n_iter_

Number of EM iterations of the kept initialisation.

Type:: int

``means_``, ``variances_``, ``loglik_`` and ``bic_`` refer to the normalised

data when ``normalize=True``. The labels returned by ``predict`` are numbered

by order of first appearance, so label ``k`` is not necessarily row ``k`` of

these arrays.

Notes

One EM iteration costs $O(nK^2 + nKd)$ for $n$ time steps, $K$ states and $d$ channels. Selecting $K$ by BIC multiplies this by the number of candidates. The forward-backward and Viterbi recursions use numba when it is installed (tsseg[accelerators]) and run, much slower, in pure Python otherwise.

The emissions are i.i.d. given the state, so this model has no notion of duration or of temporal shape within a state: it suits regimes that differ in level or spread, not regimes that differ only in frequency or waveform. On oscillating signals, the states follow the oscillation and the labels change within a regime. For the same reason, BIC often selects max_states on long autocorrelated series; a warning says so.

On integer or coarsely quantised data, a state can collapse onto a single value, with a variance of reg_covar; raise reg_covar in that case.

The initial probabilities are estimated from one series and are often close to one-hot: on another series that starts in a different regime, the first labels can be wrong. A state that receives no posterior mass during EM stays empty, so fewer than n_states labels may appear. reg_covar is added to the variances rather than used as a floor, so EM is not strictly monotone; a decrease of the log-likelihood stops it.

References

Examples

>>> import numpy as np
>>> from tsseg.algorithms import GaussianHMMDetector
>>> rng = np.random.default_rng(0)
>>> X = np.concatenate([rng.normal(0, 1, 200), rng.normal(5, 1, 200)])
>>> labels = GaussianHMMDetector(n_states=2, random_state=0).fit_predict(X)
>>> int(labels[:200].mean() != labels[200:].mean())
1

set_fit_request(*, axis: bool | None | str = '$UNCHANGED$') → GaussianHMMDetector

Configure whether metadata should be requested to be passed to the fit method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:: axis (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for axis parameter in fit.
Returns:: self – The updated object.
Return type:: object

set_predict_request(*, axis: bool | None | str = '$UNCHANGED$') → GaussianHMMDetector

Configure whether metadata should be requested to be passed to the predict method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to predict if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to predict.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:: axis (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for axis parameter in predict.
Returns:: self – The updated object.
Return type:: object