yahmm

Yet Another Hidden Markov Model library

<b>NOTE: While yahmm is still fully functional, active development has moved over to pomegranate. Please switch over at your convenience.</b>

This module implements Hidden Markov Models (HMMs) with a compositional, graph- based interface. Models can be constructed node by node and edge by edge, built up from smaller models, loaded from files, baked (into a form that can be used to calculate probabilities efficiently), trained on data, and saved.

Implements the forwards, backwards, forward-backward, and Viterbi algorithms, and training by both Baum-Welch and Viterbi algorithms.

Silent states are accounted for, but loops containing all silent states are prohibited. Tied states are also implemented, and handled appropriately in the training of models.

Installation

Since yahmm is on PyPi, installation is as easy as running

pip install yahmm

Contributing

If you would like to contribute a feature then fork the master branch (fork the release if you are fixing a bug). Be sure to run the tests before changing any code. You'll need to have nosetests installed. The following command will run all the tests:

nosetests -w tests/

Let us know what you want to do just in case we're already working on an implementation of something similar. This way we can avoid any needless duplication of effort. Also, please don't forget to add tests for any new functions.

Documentation

See the wiki for documentation of yahmm's functions and design. For real-world usage check out the examples.

Tutorial

For our examples here we're going to make the random number generator deterministic:

>>> random.seed(0)

To use this module, first create a Model, which is the main HMM class:

>>> model = Model(name="ExampleModel")

You then need to populate the Model with State objects. States are constructed from emission distributions; right now a few continuous distributions over floats are available, but new Distribution classes are simple to write. For our example, we will use the UniformDistribution:

>>> distribution = UniformDistribution(0.0, 1.0)

And then construct a state that emits from the distribution:

>>> state = State(distribution, name="uniform")

And another state, emitting from a normal distribution with mean 0 and standard deviation 2:

>>> state2 = State(NormalDistribution(0, 2), name="normal")

If None is used as the distribution when creating a state, that state is a "silent state". Silent states don't emit anything, but are useful for wiring together complex HMMs. By default, a model has two special silent states: a start state Model.start, and an end state Model.end.

Topologies which include cycles of only silent states are prohibited; most HMM algorithms cannot process them.

>>> silent = State(None, name="silent")

We then add states to the HMM with the Model.add_state method:

>>> model.add_state(state)
>>> model.add_state(state2)

You can then add transitions between states, with associated probabilities. Out-edge probabilities are normalized to 1. for every state when the model is baked, not before.

>>> model.add_transition(state, state, 0.4)
>>> model.add_transition(state, state2, 0.4)
>>> model.add_transition(state2, state2, 0.4)
>>> model.add_transition(state2, state, 0.4)

Don't forget transitions in from the start state and out to the end state:

>>> model.add_transition(model.start, state, 0.5)
>>> model.add_transition(model.start, state2, 0.5)
>>> model.add_transition(state, model.end, 0.2)
>>> model.add_transition(state2, model.end, 0.2)

If you want to look at your model, try Model.draw(). Note that this unfortunately cannot plot self loops. If you want to do a better job of drawing the model, the underlying HMM graph is accessible as the graph attribute of the model object.

If you want to compose two Models together, use the Model.add_model() method. Note that you should not try to use the added model separately once you do this. You can also make use of the Model.concatenate_model() method, which will assume you simply want to connect model_a.end to model_b.start with a 1. probability edge.

Once we've finished building our model, we have to bake it. Internally, baking the model generates the transition log probability matrix, and imposes a numerical ordering on the states. If you add more states to the model, you will have to bake it again. Baking is also where edge normalization occurs to ensure that the out-edges for all nodes (except Model.end) sum to 1. Lastly, a simplification of the graph occurs here, merging any silent states which are connected simply by a 1.0 probability edge, as they cannot add value to the graph. You may toggle 'verbose=True' in the bake method to get a log of when either change occurs to your graph.

>>> model.bake()

Now that our model is complete, we can generate an example sequence from it:

>>> sequence = model.sample()
>>> sequence
[0.7579544029403025, 0.25891675029296335, 0.4049341374504143, \
0.30331272607892745, 0.5833820394550312]

And another:

>>> model.sample()
[0.28183784439970383, 0.6183689966753316, -2.411068768608379]

And another:

>>> model.sample()
[0.47214271545271336, -0.5804485412450214]

We can calculate the log probability of the sequence given the model (the log likelihood), summing over all possible paths, using both the forward and backward algorithms. Log probability is reported in nats (i.e. it is natural log probability). Both algorithms return the full table of size len( observations ) x len( states ). For the forward algorithm, the entry at position i, j represents the log probability of beginning at the start of the sequence, and summing over all paths to align observation i to hidden state j. This state can be recovered by pulling it from model.states[j].

>>> model.forward(sequence)
[[       -inf        -inf        -inf  0.        ]
 [-2.37704475 -0.69314718 -2.1322948         -inf]
 [-3.05961307 -1.43914762 -2.86809348        -inf]
 [-3.80752847 -2.1749463  -3.60588302        -inf]
 [-4.53632138 -2.91273584 -4.34219628        -inf]
 [-5.30367664 -3.6490491  -5.08355666        -inf]]

In order to get the log probability of the full sequence given the model, you can write the following:

>>> model.forward(sequence)[ len(sequence), model.end_index ]
-5.0835566645

Or, use a wrapper to get that value by default:

>>> model.log_probability(sequence)
-5.0835566645

The same paradigm is used for the backward algorithm. Indices i, j represent the probability of having aligned observation i to state j and continued aligning the remainder of the sequence till the end.

>>> model.backward(sequence)
[[-5.30670022 -5.30670022        -inf -5.08355666]
 [-4.56069977 -4.56069977        -inf -4.33755622]
 [-3.8249011  -3.8249011         -inf -3.60175755]
 [-3.08711156 -3.08711156        -inf -2.863968  ]
 [-2.3507983  -2.3507983         -inf -2.12765475]
 [-1.60943791 -1.60943791  0.                -inf]]

>>> model.backward(sequence)[ 0, model.start_index ]
-5.0835566645

The forward-backward algorithm is also implemented in a similar manner. It will return a tuple of the estimated transition probabilities given with that sequence and the table of log probabilities of the sum of all paths of the alignment of observation i with state j. Indices i, j represent having started at the beginning of the sequence, aligned observation i to state j, and then continued on to align the remainder of the sequence to the model.

>>> model.forward_backward(sequence)
(array([[-2.03205947, -0.39913252, -1.61932212,        -inf],
       [-2.03481952, -0.40209763, -1.60753724,        -inf],
       [       -inf,        -inf,        -inf,        -inf],
       [-1.85418786, -0.17029029,        -inf,        -inf]]), 
array([[-1.85418786, -0.17029029,  0.        ,  0.        ],
       [-1.80095751, -0.18049206,  0.        ,  0.        ],
       [-1.81108336, -0.17850119,  0.        ,  0.        ],
       [-1.80356301, -0.17997747,  0.        ,  0.        ],
       [-1.82955788, -0.17493035,  0.        ,  0.        ]]))

We can also find the most likely path, and the probability thereof, using the Viterbi algorithm. This returns a tuple of the likelihood under the ML path and the ML path itself. The ML path is in turn a list of tuples of State objects and the number of items in the sequence that had been generated by that point in the path (to account for the presence of silent states).

>>> model.viterbi(sequence)
(-5.9677480204906654, \
[(0, State(ExampleModel-start, None)), \
(1, State(uniform, UniformDistribution(0.0, 1.0))), \
(2, State(uniform, UniformDistribution(0.0, 1.0))), \
(3, State(uniform, UniformDistribution(0.0, 1.0))), \
(4, State(uniform, UniformDistribution(0.0, 1.0))), \
(5, State(uniform, UniformDistribution(0.0, 1.0))), \
(5, State(ExampleModel-end, None))])

Given a list of sequences, we can train our HMM by calling Model.train(). This returns the final log score: the log of the sum of the probabilities of all training sequences. It also prints the improvement in log score on each training iteration, and stops if the improvement gets too small or actually goes negative.

>>> sequences = [sequence]
>>> model.log_probability(sequence)
-5.0835566644993735

>>> log_score = model.train(sequences)
Training improvement: 5.81315226327
Training improvement: 0.156159401683
Training improvement: 0.0806734819188
Training improvement: 0.0506679952827
Training improvement: 0.142593661095
Training im