Hidden Markov Models#

This is a complete pure-Cython optimized implementation of Hidden Markov Models. It fully supports Discrete, Gaussian, and Mixed Gaussian emissions.

The best references for the basic HMM algorithms implemented here are:

LICENSE: Some of the code in this file is based on reading Kanungo’s GPLv2+ implementation of discrete HMM’s, hence the present code must be licensed with a GPLv2+ compatible license.

AUTHOR:

  • William Stein, 2010-03

class sage.stats.hmm.hmm.DiscreteHiddenMarkovModel[source]#

Bases: HiddenMarkovModel

A discrete Hidden Markov model implemented using double precision floating point arithmetic.

INPUT:

  • A – a list of lists or a square \(N \times N\) matrix, whose \((i,j)\) entry gives the probability of transitioning from state \(i\) to state \(j\).

  • B – a list of \(N\) lists or a matrix with \(N\) rows, such that \(B[i,k]\) gives the probability of emitting symbol \(k\) while in state \(i\).

  • pi – the probabilities of starting in each initial state, i.e., pi[i] is the probability of starting in state \(i\).

  • emission_symbols – None or list (default: None); if None, the emission_symbols are the ints [0..N-1], where \(N\) is the number of states. Otherwise, they are the entries of the list emissions_symbols, which must all be hashable.

  • normalize – bool (default: True); if given, input is normalized to define valid probability distributions, e.g., the entries of \(A\) are made nonnegative and the rows sum to 1, and the probabilities in pi are normalized.

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.5,.5]); m
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[0.4 0.6]
[0.1 0.9]
Emission matrix:
[0.1 0.9]
[0.5 0.5]
Initial probabilities: [0.5000, 0.5000]
sage: m.log_likelihood([0,1,0,1,0,1])
-4.66693474691329...
sage: m.viterbi([0,1,0,1,0,1])
([1, 1, 1, 1, 1, 1], -5.378832842208748)
sage: m.baum_welch([0,1,0,1,0,1])
(0.0, 22)
sage: m  # rel tol 1e-10
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[1.0134345614745788e-70                    1.0]
[                   1.0 3.9974352713558623e-19]
Emission matrix:
[ 7.380221566254936e-54                    1.0]
[                   1.0 3.9974352626002193e-19]
Initial probabilities: [0.0000, 1.0000]
sage: m.sample(10)
[0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
sage: m.graph().plot()                                                          # needs sage.plot
Graphics object consisting of 6 graphics primitives

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.5'),RealNumber('.5')]); m
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[0.4 0.6]
[0.1 0.9]
Emission matrix:
[0.1 0.9]
[0.5 0.5]
Initial probabilities: [0.5000, 0.5000]
>>> m.log_likelihood([Integer(0),Integer(1),Integer(0),Integer(1),Integer(0),Integer(1)])
-4.66693474691329...
>>> m.viterbi([Integer(0),Integer(1),Integer(0),Integer(1),Integer(0),Integer(1)])
([1, 1, 1, 1, 1, 1], -5.378832842208748)
>>> m.baum_welch([Integer(0),Integer(1),Integer(0),Integer(1),Integer(0),Integer(1)])
(0.0, 22)
>>> m  # rel tol 1e-10
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[1.0134345614745788e-70                    1.0]
[                   1.0 3.9974352713558623e-19]
Emission matrix:
[ 7.380221566254936e-54                    1.0]
[                   1.0 3.9974352626002193e-19]
Initial probabilities: [0.0000, 1.0000]
>>> m.sample(Integer(10))
[0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
>>> m.graph().plot()                                                          # needs sage.plot
Graphics object consisting of 6 graphics primitives

A 3-state model that happens to always outputs ‘b’:

sage: m = hmm.DiscreteHiddenMarkovModel([[1/3]*3]*3, [[0,1,0]]*3, [1/3]*3, ['a','b','c'])
sage: m.sample(10)
['b', 'b', 'b', 'b', 'b', 'b', 'b', 'b', 'b', 'b']

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[Integer(1)/Integer(3)]*Integer(3)]*Integer(3), [[Integer(0),Integer(1),Integer(0)]]*Integer(3), [Integer(1)/Integer(3)]*Integer(3), ['a','b','c'])
>>> m.sample(Integer(10))
['b', 'b', 'b', 'b', 'b', 'b', 'b', 'b', 'b', 'b']
baum_welch(obs, max_iter=100, log_likelihood_cutoff=0.0001, fix_emissions=False)[source]#

Given an observation sequence obs, improve this HMM using the Baum-Welch algorithm to increase the probability of observing obs.

INPUT:

  • obs – list of emissions

  • max_iter – integer (default: 100) maximum number of Baum-Welch steps to take

  • log_likelihood_cutoff – positive float (default: 1e-4); the minimal improvement in likelihood with respect to the last iteration required to continue. Relative value to log likelihood.

  • fix_emissions – bool (default: False); if True, do not change emissions when updating

OUTPUT:

changes the model in place, and returns the log likelihood and number of iterations.

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]],
....:                                   [[.5,.5],[0,1]],
....:                                   [.2,.8])
sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0)
(0.0, 4)
sage: m  # rel tol 1e-14
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[1.3515269707707603e-51                    1.0]
[                   1.0                    0.0]
Emission matrix:
[                  1.0 6.462537138850569e-52]
[                  0.0                   1.0]
Initial probabilities: [0.0000, 1.0000]

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.9'),RealNumber('0.1')]],
...                                   [[RealNumber('.5'),RealNumber('.5')],[Integer(0),Integer(1)]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> m.baum_welch([Integer(1),Integer(0)]*Integer(20), log_likelihood_cutoff=Integer(0))
(0.0, 4)
>>> m  # rel tol 1e-14
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[1.3515269707707603e-51                    1.0]
[                   1.0                    0.0]
Emission matrix:
[                  1.0 6.462537138850569e-52]
[                  0.0                   1.0]
Initial probabilities: [0.0000, 1.0000]

The following illustrates how Baum-Welch is only a local optimizer, i.e., the above model is far more likely to produce the sequence [1,0]*20 than the one we get below:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]],
....:                                   [[.5,.5],[.5,.5]],
....:                                   [.5,.5])
sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0)
(-27.725887222397784, 1)
sage: m
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[0.5 0.5]
[0.5 0.5]
Emission matrix:
[0.5 0.5]
[0.5 0.5]
Initial probabilities: [0.5000, 0.5000]

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.5'),RealNumber('0.5')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [[RealNumber('.5'),RealNumber('.5')],[RealNumber('.5'),RealNumber('.5')]],
...                                   [RealNumber('.5'),RealNumber('.5')])
>>> m.baum_welch([Integer(1),Integer(0)]*Integer(20), log_likelihood_cutoff=Integer(0))
(-27.725887222397784, 1)
>>> m
Discrete Hidden Markov Model with 2 States and 2 Emissions
Transition matrix:
[0.5 0.5]
[0.5 0.5]
Emission matrix:
[0.5 0.5]
[0.5 0.5]
Initial probabilities: [0.5000, 0.5000]

We illustrate fixing emissions:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]],
....:                                   [[.5,.5],[.2,.8]],
....:                                   [.2,.8])
sage: set_random_seed(0); v = m.sample(100)
sage: m.baum_welch(v,fix_emissions=True)
(-66.98630856918774, 100)
sage: m.emission_matrix()
[0.5 0.5]
[0.2 0.8]
sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]],
....:                                   [[.5,.5],[.2,.8]],
....:                                   [.2,.8])
sage: m.baum_welch(v)
(-66.782360659293..., 100)
sage: m.emission_matrix()  # rel tol 1e-14
[ 0.5303085748626447 0.46969142513735535]
[ 0.2909775550173978  0.7090224449826023]

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.9'),RealNumber('0.1')]],
...                                   [[RealNumber('.5'),RealNumber('.5')],[RealNumber('.2'),RealNumber('.8')]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> set_random_seed(Integer(0)); v = m.sample(Integer(100))
>>> m.baum_welch(v,fix_emissions=True)
(-66.98630856918774, 100)
>>> m.emission_matrix()
[0.5 0.5]
[0.2 0.8]
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.9'),RealNumber('0.1')]],
...                                   [[RealNumber('.5'),RealNumber('.5')],[RealNumber('.2'),RealNumber('.8')]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> m.baum_welch(v)
(-66.782360659293..., 100)
>>> m.emission_matrix()  # rel tol 1e-14
[ 0.5303085748626447 0.46969142513735535]
[ 0.2909775550173978  0.7090224449826023]
emission_matrix()[source]#

Return the matrix whose \(i\)-th row specifies the emission probability distribution for the \(i\)-th state.

More precisely, the \(i,j\) entry of the matrix is the probability of the Markov model outputting the \(j\)-th symbol when it is in the \(i\)-th state.

OUTPUT: a Sage matrix with real double precision (RDF) entries.

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.5,.5])
sage: E = m.emission_matrix(); E
[0.1 0.9]
[0.5 0.5]

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.5'),RealNumber('.5')])
>>> E = m.emission_matrix(); E
[0.1 0.9]
[0.5 0.5]

The returned matrix is mutable, but changing it does not change the transition matrix for the model:

sage: E[0,0] = 0; E[0,1] = 1
sage: m.emission_matrix()
[0.1 0.9]
[0.5 0.5]

>>> from sage.all import *
>>> E[Integer(0),Integer(0)] = Integer(0); E[Integer(0),Integer(1)] = Integer(1)
>>> m.emission_matrix()
[0.1 0.9]
[0.5 0.5]
generate_sequence(length, starting_state=None)[source]#

Return a sample of the given length from this HMM.

INPUT:

  • length – positive integer

  • starting_state – int (or None); if specified, generate a sequence using this model starting with the given state instead of the initial probabilities to determine the starting state.

OUTPUT:

  • an IntList or list of emission symbols

  • IntList of the actual states the model was in when emitting the corresponding symbols

EXAMPLES:

In this example, the emission symbols are not set:

sage: set_random_seed(0)
sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]],
....:                                   [[1,0],[0,1]],
....:                                   [0,1])
sage: a.generate_sequence(5)
([1, 0, 1, 1, 1], [1, 0, 1, 1, 1])
sage: list(a.generate_sequence(1000)[0]).count(0)
90

>>> from sage.all import *
>>> set_random_seed(Integer(0))
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[Integer(1),Integer(0)],[Integer(0),Integer(1)]],
...                                   [Integer(0),Integer(1)])
>>> a.generate_sequence(Integer(5))
([1, 0, 1, 1, 1], [1, 0, 1, 1, 1])
>>> list(a.generate_sequence(Integer(1000))[Integer(0)]).count(Integer(0))
90

Here the emission symbols are set:

sage: set_random_seed(0)
sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]],
....:                                   [[1,0],[0,1]],
....:                                   [0,1], ['up', 'down'])
sage: a.generate_sequence(5)
(['down', 'up', 'down', 'down', 'down'], [1, 0, 1, 1, 1])

>>> from sage.all import *
>>> set_random_seed(Integer(0))
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.5'),RealNumber('0.5')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[Integer(1),Integer(0)],[Integer(0),Integer(1)]],
...                                   [Integer(0),Integer(1)], ['up', 'down'])
>>> a.generate_sequence(Integer(5))
(['down', 'up', 'down', 'down', 'down'], [1, 0, 1, 1, 1])

Specify the starting state:

sage: set_random_seed(0); a.generate_sequence(5, starting_state=0)
(['up', 'up', 'down', 'down', 'down'], [0, 0, 1, 1, 1])

>>> from sage.all import *
>>> set_random_seed(Integer(0)); a.generate_sequence(Integer(5), starting_state=Integer(0))
(['up', 'up', 'down', 'down', 'down'], [0, 0, 1, 1, 1])
log_likelihood(obs, scale=True)[source]#

Return the logarithm of the probability that this model produced the given observation sequence. Thus the output is a non-positive number.

INPUT:

  • obs – sequence of observations

  • scale – boolean (default: True); if True, use rescaling to overoid loss of precision due to the very limited dynamic range of floats. You should leave this as True unless the obs sequence is very small.

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.2,.8])
sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0])
-7.3301308009370825
sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0], scale=False)
-7.330130800937082
sage: m.log_likelihood([])
0.0

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.2,.8], ['happy','sad'])
sage: m.log_likelihood(['happy','happy'])
-1.6565295199679506
sage: m.log_likelihood(['happy','sad'])
-1.4731602941415523

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> m.log_likelihood([Integer(0), Integer(1), Integer(0), Integer(1), Integer(1), Integer(0), Integer(1), Integer(0), Integer(0), Integer(0)])
-7.3301308009370825
>>> m.log_likelihood([Integer(0), Integer(1), Integer(0), Integer(1), Integer(1), Integer(0), Integer(1), Integer(0), Integer(0), Integer(0)], scale=False)
-7.330130800937082
>>> m.log_likelihood([])
0.0

>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.2'),RealNumber('.8')], ['happy','sad'])
>>> m.log_likelihood(['happy','happy'])
-1.6565295199679506
>>> m.log_likelihood(['happy','sad'])
-1.4731602941415523

Overflow from not using the scale option:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.2,.8])
sage: m.log_likelihood([0,1]*1000, scale=True)
-1433.820666652728
sage: m.log_likelihood([0,1]*1000, scale=False)
-inf

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> m.log_likelihood([Integer(0),Integer(1)]*Integer(1000), scale=True)
-1433.820666652728
>>> m.log_likelihood([Integer(0),Integer(1)]*Integer(1000), scale=False)
-inf
viterbi(obs, log_scale=True)[source]#

Determine “the” hidden sequence of states that is most likely to produce the given sequence seq of observations, along with the probability that this hidden sequence actually produced the observation.

INPUT:

  • seq – sequence of emitted ints or symbols

  • log_scale – bool (default: True) whether to scale the sequence in order to avoid numerical overflow.

OUTPUT:

  • list – “the” most probable sequence of hidden states, i.e., the Viterbi path.

  • float – log of probability that the observed sequence was produced by the Viterbi sequence of states.

EXAMPLES:

sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]],
....:                                   [[0.9,0.1],[0.1,0.9]],
....:                                   [0.5,0.5])
sage: a.viterbi([1,0,0,1,0,0,1,1])
([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...)

>>> from sage.all import *
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.9'),RealNumber('0.1')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [RealNumber('0.5'),RealNumber('0.5')])
>>> a.viterbi([Integer(1),Integer(0),Integer(0),Integer(1),Integer(0),Integer(0),Integer(1),Integer(1)])
([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...)

We predict the state sequence when the emissions are 3/4 and ‘abc’.:

sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]],
....:                                   [[0.9,0.1],[0.1,0.9]],
....:                                   [0.5,0.5], [3/4, 'abc'])

>>> from sage.all import *
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.9'),RealNumber('0.1')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [RealNumber('0.5'),RealNumber('0.5')], [Integer(3)/Integer(4), 'abc'])

Note that state 0 is common below, despite the model trying hard to switch to state 1:

sage: a.viterbi([3/4, 'abc', 'abc'] + [3/4]*10)
([0, 1, 1, 0, 0 ... 0, 0, 0, 0, 0], -25.299405845367794)

>>> from sage.all import *
>>> a.viterbi([Integer(3)/Integer(4), 'abc', 'abc'] + [Integer(3)/Integer(4)]*Integer(10))
([0, 1, 1, 0, 0 ... 0, 0, 0, 0, 0], -25.299405845367794)
class sage.stats.hmm.hmm.HiddenMarkovModel[source]#

Bases: object

Abstract base class for all Hidden Markov Models.

graph(eps=0.001)[source]#

Create a weighted directed graph from the transition matrix, not including any edge with a probability less than eps.

INPUT:

  • eps – nonnegative real number

OUTPUT: a DiGraph

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[.3,0,.7],[0,0,1],[.5,.5,0]],
....:                                   [[.5,.5,.2]]*3,
....:                                   [1/3]*3)
sage: G = m.graph(); G                                                      # needs sage.graphs
Looped digraph on 3 vertices
sage: G.edges(sort=True)                                                    # needs sage.graphs
[(0, 0, 0.3), (0, 2, 0.7), (1, 2, 1.0), (2, 0, 0.5), (2, 1, 0.5)]
sage: G.plot()                                                              # needs sage.graphs sage.plot
Graphics object consisting of 11 graphics primitives

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('.3'),Integer(0),RealNumber('.7')],[Integer(0),Integer(0),Integer(1)],[RealNumber('.5'),RealNumber('.5'),Integer(0)]],
...                                   [[RealNumber('.5'),RealNumber('.5'),RealNumber('.2')]]*Integer(3),
...                                   [Integer(1)/Integer(3)]*Integer(3))
>>> G = m.graph(); G                                                      # needs sage.graphs
Looped digraph on 3 vertices
>>> G.edges(sort=True)                                                    # needs sage.graphs
[(0, 0, 0.3), (0, 2, 0.7), (1, 2, 1.0), (2, 0, 0.5), (2, 1, 0.5)]
>>> G.plot()                                                              # needs sage.graphs sage.plot
Graphics object consisting of 11 graphics primitives
initial_probabilities()[source]#

Return the initial probabilities as a TimeSeries of length \(N\), where \(N\) is the number of states of the Markov model.

EXAMPLES:

sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],
....:                                   [[0.1,0.9],[0.5,0.5]],
....:                                   [.2,.8])
sage: pi = m.initial_probabilities(); pi
[0.2000, 0.8000]
sage: type(pi)
<... 'sage.stats.time_series.TimeSeries'>

>>> from sage.all import *
>>> m = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.4'),RealNumber('0.6')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.5'),RealNumber('0.5')]],
...                                   [RealNumber('.2'),RealNumber('.8')])
>>> pi = m.initial_probabilities(); pi
[0.2000, 0.8000]
>>> type(pi)
<... 'sage.stats.time_series.TimeSeries'>

The returned time series is a copy, so changing it does not change the model:

sage: pi[0] = .1; pi[1] = .9
sage: m.initial_probabilities()
[0.2000, 0.8000]

>>> from sage.all import *
>>> pi[Integer(0)] = RealNumber('.1'); pi[Integer(1)] = RealNumber('.9')
>>> m.initial_probabilities()
[0.2000, 0.8000]

Some other models:

sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]],
....:                                   [(1,1), (-1,1)],
....:                                   [.1,.9])
sage: m.initial_probabilities()
[0.1000, 0.9000]
sage: m = hmm.GaussianMixtureHiddenMarkovModel(
....:         [[.9,.1],[.4,.6]],
....:         [[(.4,(0,1)), (.6,(1,0.1))], [(1,(0,1))]],
....:         [.7,.3])
sage: m.initial_probabilities()
[0.7000, 0.3000]

>>> from sage.all import *
>>> m = hmm.GaussianHiddenMarkovModel([[RealNumber('.1'),RealNumber('.9')],[RealNumber('.5'),RealNumber('.5')]],
...                                   [(Integer(1),Integer(1)), (-Integer(1),Integer(1))],
...                                   [RealNumber('.1'),RealNumber('.9')])
>>> m.initial_probabilities()
[0.1000, 0.9000]
>>> m = hmm.GaussianMixtureHiddenMarkovModel(
...         [[RealNumber('.9'),RealNumber('.1')],[RealNumber('.4'),RealNumber('.6')]],
...         [[(RealNumber('.4'),(Integer(0),Integer(1))), (RealNumber('.6'),(Integer(1),RealNumber('0.1')))], [(Integer(1),(Integer(0),Integer(1)))]],
...         [RealNumber('.7'),RealNumber('.3')])
>>> m.initial_probabilities()
[0.7000, 0.3000]
sample(length, number=None, starting_state=None)[source]#

Return number samples from this HMM of given length.

INPUT:

  • length – positive integer

  • number – (default: None) if given, compute list of this many sample sequences

  • starting_state – int (or None); if specified, generate a sequence using this model starting with the given state instead of the initial probabilities to determine the starting state.

OUTPUT:

  • if number is not given, return a single TimeSeries.

  • if number is given, return a list of TimeSeries.

EXAMPLES:

sage: set_random_seed(0)
sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]],
....:                                   [[1,0],[0,1]],
....:                                   [0,1])
sage: print(a.sample(10, 3))
[[1, 0, 1, 1, 1, 1, 0, 1, 1, 1],
 [1, 1, 0, 0, 1, 1, 1, 1, 1, 1],
 [1, 1, 1, 1, 0, 1, 0, 1, 1, 1]]
sage: a.sample(15)
[1, 1, 1, 1, 0 ... 1, 1, 1, 1, 1]
sage: a.sample(3, 1)
[[1, 1, 1]]
sage: list(a.sample(1000)).count(0)
88

>>> from sage.all import *
>>> set_random_seed(Integer(0))
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.1'),RealNumber('0.9')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[Integer(1),Integer(0)],[Integer(0),Integer(1)]],
...                                   [Integer(0),Integer(1)])
>>> print(a.sample(Integer(10), Integer(3)))
[[1, 0, 1, 1, 1, 1, 0, 1, 1, 1],
 [1, 1, 0, 0, 1, 1, 1, 1, 1, 1],
 [1, 1, 1, 1, 0, 1, 0, 1, 1, 1]]
>>> a.sample(Integer(15))
[1, 1, 1, 1, 0 ... 1, 1, 1, 1, 1]
>>> a.sample(Integer(3), Integer(1))
[[1, 1, 1]]
>>> list(a.sample(Integer(1000))).count(Integer(0))
88

If the emission symbols are set:

sage: set_random_seed(0)
sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]],
....:                                   [[1,0],[0,1]], [0,1],
....:                                   ['up', 'down'])
sage: a.sample(10)
['down', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up']

>>> from sage.all import *
>>> set_random_seed(Integer(0))
>>> a = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.5'),RealNumber('0.5')],[RealNumber('0.1'),RealNumber('0.9')]],
...                                   [[Integer(1),Integer(0)],[Integer(0),Integer(1)]], [Integer(0),Integer(1)],
...                                   ['up', 'down'])
>>> a.sample(Integer(10))
['down', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up']

Force a starting state:

sage: set_random_seed(0); a.sample(10, starting_state=0)
['up', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up']

>>> from sage.all import *
>>> set_random_seed(Integer(0)); a.sample(Integer(10), starting_state=Integer(0))
['up', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up']
transition_matrix()[source]#

Return the state transition matrix.

OUTPUT: a Sage matrix with real double precision (RDF) entries.

EXAMPLES:

sage: M = hmm.DiscreteHiddenMarkovModel([[0.7,0.3],[0.9,0.1]],
....:                                   [[0.5,.5],[.1,.9]],
....:                                   [0.3,0.7])
sage: T = M.transition_matrix(); T
[0.7 0.3]
[0.9 0.1]

>>> from sage.all import *
>>> M = hmm.DiscreteHiddenMarkovModel([[RealNumber('0.7'),RealNumber('0.3')],[RealNumber('0.9'),RealNumber('0.1')]],
...                                   [[RealNumber('0.5'),RealNumber('.5')],[RealNumber('.1'),RealNumber('.9')]],
...                                   [RealNumber('0.3'),RealNumber('0.7')])
>>> T = M.transition_matrix(); T
[0.7 0.3]
[0.9 0.1]

The returned matrix is mutable, but changing it does not change the transition matrix for the model:

sage: T[0,0] = .1; T[0,1] = .9
sage: M.transition_matrix()
[0.7 0.3]
[0.9 0.1]

>>> from sage.all import *
>>> T[Integer(0),Integer(0)] = RealNumber('.1'); T[Integer(0),Integer(1)] = RealNumber('.9')
>>> M.transition_matrix()
[0.7 0.3]
[0.9 0.1]

Transition matrices for other types of models:

sage: M = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]],
....:                                   [(1,1), (-1,1)],
....:                                   [.5,.5])
sage: M.transition_matrix()
[0.1 0.9]
[0.5 0.5]
sage: M = hmm.GaussianMixtureHiddenMarkovModel(
....:         [[.9,.1],[.4,.6]],
....:         [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]],
....:         [.7,.3])
sage: M.transition_matrix()
[0.9 0.1]
[0.4 0.6]

>>> from sage.all import *
>>> M = hmm.GaussianHiddenMarkovModel([[RealNumber('.1'),RealNumber('.9')],[RealNumber('.5'),RealNumber('.5')]],
...                                   [(Integer(1),Integer(1)), (-Integer(1),Integer(1))],
...                                   [RealNumber('.5'),RealNumber('.5')])
>>> M.transition_matrix()
[0.1 0.9]
[0.5 0.5]
>>> M = hmm.GaussianMixtureHiddenMarkovModel(
...         [[RealNumber('.9'),RealNumber('.1')],[RealNumber('.4'),RealNumber('.6')]],
...         [[(RealNumber('.4'),(Integer(0),Integer(1))), (RealNumber('.6'),(Integer(1),RealNumber('0.1')))],[(Integer(1),(Integer(0),Integer(1)))]],
...         [RealNumber('.7'),RealNumber('.3')])
>>> M.transition_matrix()
[0.9 0.1]
[0.4 0.6]
sage.stats.hmm.hmm.unpickle_discrete_hmm_v0(A, B, pi, emission_symbols, name)[source]#
sage.stats.hmm.hmm.unpickle_discrete_hmm_v1(A, B, pi, n_out, emission_symbols, emission_symbols_dict)[source]#

Return a DiscreteHiddenMarkovModel, restored from the arguments.

This function is used internally for unpickling.