.. currentmodule:: hmmlearn

hmmlearn implements the Hidden Markov Models (HMMs). The HMM is a generative probabilistic model, in which a sequence of observable \mathbf{X} variables is generated by a sequence of internal hidden states \mathbf{Z}. The hidden states are not observed directly. The transitions between hidden states are assumed to have the form of a (first-order) Markov chain. They can be specified by the start probability vector \boldsymbol{\pi} and a transition probability matrix \mathbf{A}. The emission probability of an observable can be any distribution with parameters \boldsymbol{\theta} conditioned on the current hidden state. The HMM is completely determined by \boldsymbol{\pi}, \mathbf{A} and \boldsymbol{\theta}.

There are three fundamental problems for HMMs:

• Given the model parameters and observed data, estimate the optimal sequence of hidden states.
• Given the model parameters and observed data, calculate the likelihood of the data.
• Given just the observed data, estimate the model parameters.

The first and the second problem can be solved by the dynamic programming algorithms known as the Viterbi algorithm and the Forward-Backward algorithm, respectively. The last one can be solved by an iterative Expectation-Maximization (EM) algorithm, known as the Baum-Welch algorithm.

.. autosummary::
   :nosignatures:

   hmm.GaussianHMM
   hmm.GMMHMM
   hmm.MultinomialHMM

:ref:`Read on <customizing>` for details on how to implement a HMM with a custom emission probability.

You can build a HMM instance by passing the parameters described above to the constructor. Then, you can generate samples from the HMM by calling :meth:`~base._BaseHMM.sample`.

>>> import numpy as np
>>> from hmmlearn import hmm
>>> np.random.seed(42)

>>> model = hmm.GaussianHMM(n_components=3, covariance_type="full")
>>> model.startprob_ = np.array([0.6, 0.3, 0.1])
>>> model.transmat_ = np.array([[0.7, 0.2, 0.1],
...                             [0.3, 0.5, 0.2],
...                             [0.3, 0.3, 0.4]])
>>> model.means_ = np.array([[0.0, 0.0], [3.0, -3.0], [5.0, 10.0]])
>>> model.covars_ = np.tile(np.identity(2), (3, 1, 1))
>>> X, Z = model.sample(100)

The transition probability matrix need not to be ergodic. For instance, a left-right HMM can be defined as follows:

>>> lr = hmm.GaussianHMM(n_components=3, covariance_type="diag",
...                      init_params="cm", params="cmt")
>>> lr.startprob_ = np.array([1.0, 0.0, 0.0])
>>> lr.transmat_ = np.array([[0.5, 0.5, 0.0],
...                          [0.0, 0.5, 0.5],
...                          [0.0, 0.0, 1.0]])

If any of the required parameters are missing, :meth:`~base._BaseHMM.sample` will raise an exception:

>>> hmm.GaussianHMM(n_components=3)
>>> X, Z = model.sample(100)
Traceback (most recent call last):
    ...
sklearn.utils.validation.NotFittedError: This GaussianHMM instance is not fitted yet. Call 'fit' with appropriate arguments before using this method.

You can train an HMM by calling the :meth:`~base._BaseHMM.fit` method. The input is a matrix of concatenated sequences of observations (aka samples) along with the lengths of the sequences (see :ref:`Working with multiple sequences <multiple_sequences>`).

Note, since the EM algorithm is a gradient-based optimization method, it will generally get stuck in local optima. You should in general try to run fit with various initializations and select the highest scored model.

The score of the model can be calculated by the :meth:`~base._BaseHMM.score` method.

The inferred optimal hidden states can be obtained by calling :meth:`~base._BaseHMM.predict` method. The predict method can be specified with a decoder algorithm. Currently the Viterbi algorithm ("viterbi"), and maximum a posteriori estimation ("map") are supported.

This time, the input is a single sequence of observed values. Note, the states in remodel will have a different order than those in the generating model.

>>> remodel = hmm.GaussianHMM(n_components=3, covariance_type="full", n_iter=100)
>>> remodel.fit(X)  # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
GaussianHMM(algorithm='viterbi',...
>>> Z2 = remodel.predict(X)

>>> remodel.monitor_  # doctest: +ELLIPSIS
ConvergenceMonitor(history=[...],
          iter=12, n_iter=100, tol=0.01, verbose=False)
>>> remodel.monitor_.converged
True

>>> X1 = [[0.5], [1.0], [-1.0], [0.42], [0.24]]
>>> X2 = [[2.4], [4.2], [0.5], [-0.24]]

>>> X = np.concatenate([X1, X2])
>>> lengths = [len(X1), len(X2)]

>>> hmm.GaussianHMM(n_components=3).fit(X, lengths)  # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
GaussianHMM(algorithm='viterbi', ...

After training, a HMM can be easily persisted for future use with the standard :mod:`pickle` module or its more efficient replacement in the :mod:`joblib` package:

>>> from sklearn.externals import joblib
>>> joblib.dump(remodel, "filename.pkl")
["filename.pkl"]
>>> joblib.load("filename.pkl")  # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
GaussianHMM(algorithm='viterbi',...

If you want to implement a custom emission probability (e.g. Poisson), you have to subclass :class:`~base._BaseHMM` and override the following methods

.. autosummary::

   base._BaseHMM._init
   base._BaseHMM._check
   base._BaseHMM._generate_sample_from_state
   base._BaseHMM._compute_log_likelihood
   base._BaseHMM._initialize_sufficient_statistics
   base._BaseHMM._accumulate_sufficient_statistics
   base._BaseHMM._do_mstep