Infinite Hidden Markov Model

authors:
Jacob Schreiber [jmschreiber91@gmail.com]
Nicholas Farn [nicholasfarn@gmail.com]

This example shows how to use pomegranate to sample from an infinite HMM. The premise is that you have an HMM which does not have transitions to the end state, and so can continue on forever. This is done by not adding transitions to the end state. If you bake a model with no transitions to the end state, you get an infinite model, with no extra work! This change is passed on to all the algorithms.


In [1]:
from pomegranate import *
import itertools as it
import numpy as np

First we define the possible states in the model. In this case we make them all have normal distributions.


In [2]:
s1 = State( NormalDistribution( 5, 2 ), name="S1" )
s2 = State( NormalDistribution( 15, 2 ), name="S2" )
s3 = State( NormalDistribution( 25, 2 ), name="S3" )

We then create the HMM object, naming it, logically, "infinite".


In [3]:
model = HiddenMarkovModel( "infinite" )

We then add the possible transition, making sure not to add an end state. Thus with no end state, the model is infinite!


In [4]:
model.add_transition( model.start, s1, 0.7 )
model.add_transition( model.start, s2, 0.2 )
model.add_transition( model.start, s3, 0.1 )
model.add_transition( s1, s1, 0.6 )
model.add_transition( s1, s2, 0.1 )
model.add_transition( s1, s3, 0.3 )
model.add_transition( s2, s1, 0.4 )
model.add_transition( s2, s2, 0.4 )
model.add_transition( s2, s3, 0.2 )
model.add_transition( s3, s1, 0.05 )
model.add_transition( s3, s2, 0.15 )
model.add_transition( s3, s3, 0.8 )

Finally we "bake" the model, finalizing the model.


In [5]:
model.bake()

Now we can check whether or not our model is infinite.


In [6]:
# Not implemented: print model.is_infinite()

Now lets the possible states in the model.


In [7]:
print("States")
print("\n".join( state.name for state in model.states ))


States
S1
S2
S3
infinite-start
infinite-end

Now lets test out our model by feeding it a sequence of values. We feed our sequence of values first through a forward algorithm in our HMM.


In [8]:
sequence = [ 4.8, 5.6, 24.1, 25.8, 14.3, 26.5, 15.9, 5.5, 5.1 ]

print("Forward")
print(model.forward( sequence ))


Forward
[[        -inf         -inf         -inf   0.                 -inf]
 [ -1.97376066 -16.22652363 -54.91967081         -inf         -inf]
 [ -4.14167156 -16.93342888 -51.83481875         -inf         -inf]
 [-51.86583105 -18.40758124  -7.05897823         -inf         -inf]
 [-65.74670193 -25.1481525   -8.97420455         -inf         -inf]
 [-24.39327178 -12.54465999 -25.12068379         -inf         -inf]
 [-72.85427528 -31.60428336 -16.04740908         -inf         -inf]
 [-35.50647567 -19.65786431 -28.2338883          -inf         -inf]
 [-22.21746699 -33.46742001 -70.40988373         -inf         -inf]
 [-24.34161966 -38.38333577 -74.53476684         -inf         -inf]]

That looks good as well. Now lets feed our sequence into the model through a backwards algorithm.


In [9]:
print("Backward")
print(model.backward( sequence ))


Backward
[[-24.49576972 -24.90123357 -26.9806721  -24.34161886         -inf]
 [-22.36785863 -22.77331446 -24.85273373 -22.21370663         -inf]
 [-20.19994915 -20.60539952 -19.21912118 -21.29854671         -inf]
 [-18.26347263 -18.6689355  -17.28264358 -19.36208268         -inf]
 [-15.77270345 -14.38657639 -15.36741456 -15.07964315         -inf]
 [-11.39151868 -11.7969834  -10.41068947 -12.49013058         -inf]
 [ -8.69968233  -7.31340196  -8.29421024  -8.00654767         -inf]
 [ -4.27832043  -4.68377831  -6.76320251  -4.12416872         -inf]
 [ -2.12416054  -2.52962166  -4.60905363  -1.97000929         -inf]
 [  0.           0.           0.                 -inf         -inf]]

Continuing on we now feed the sequence in through a forward-backward algorithm.


In [10]:
print("Forward-Backward")
trans, emissions = model.forward_backward( sequence )
print(trans)
print(emissions)


Forward-Backward
[[1.99998844e+00 5.59973220e-06 1.00000592e+00 0.00000000e+00
  0.00000000e+00]
 [9.99976663e-01 9.85835713e-06 9.99979797e-01 0.00000000e+00
  0.00000000e+00]
 [3.44897138e-05 1.99995123e+00 1.00004800e+00 0.00000000e+00
  0.00000000e+00]
 [9.99999569e-01 4.30542744e-07 8.44725590e-25 0.00000000e+00
  0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00]]
[[-4.30542826e-07 -1.46582192e+01 -5.54307857e+01]
 [-1.85577414e-06 -1.31972095e+01 -4.67123211e+01]
 [-4.57876848e+01 -1.27348979e+01 -2.94648857e-06]
 [-5.71777865e+01 -1.51931100e+01 -2.52183224e-07]
 [-1.14431716e+01 -2.45377600e-05 -1.11897544e+01]
 [-5.72123388e+01 -1.45760665e+01 -4.67406625e-07]
 [-1.54431772e+01 -2.37681755e-05 -1.06554720e+01]
 [-8.67193612e-06 -1.16554228e+01 -5.06773185e+01]
 [-7.97553827e-07 -1.40417169e+01 -5.01931480e+01]]

Finally we feed the sequence through a Viterbi algorithm to find the most probable sequence of states.


In [11]:
print("Viterbi")
prob, states = model.viterbi( sequence )
print("Prob: {}".format( prob ))
print("\n".join( state[1].name for state in states ))
print()
print("MAP")
prob, states = model.maximum_a_posteriori( sequence )
print("Prob: {}".format( prob ))
print("\n".join( state[1].name for state in states ))


Viterbi
Prob: -24.341682585068853
infinite-start
S1
S1
S3
S3
S2
S3
S2
S1
S1

MAP
Prob: -6.372782082308959e-05
S1
S1
S3
S3
S2
S3
S2
S1
S1

Finally we try and reproduce the transition matrix from 100,000 samples.


In [ ]:
print("Should produce a matrix close to the following: ")
print(" [ [ 0.60, 0.10, 0.30 ] ")
print("   [ 0.40, 0.40, 0.20 ] ")
print("   [ 0.05, 0.15, 0.80 ] ] ")
print()
print("Transition Matrix From 100000 Samples:")
sample, path = model.sample( 100000, path=True )
trans = np.zeros((3,3))

for state, n_state in it.izip( path[1:-2], path[2:-1] ):
	state_name = int( state.name[1:] )-1
	n_state_name = int( n_state.name[1:] )-1
	trans[ state_name, n_state_name ] += 1

trans = (trans.T / trans.sum( axis=1 )).T
print(trans)