This example demonstrates a Hidden Markov Model (HMM) using discrete Bayes nets.
A Hidden Markov Model consists of:
A sequence of hidden states over time
Transition probabilities between states
Observation/measurement probabilities for each state
In this example:
We have 4 nodes (time steps) with 3 possible states each
We define a transition model that favors staying in the same state
We add measurements/observations at some (not all) time steps
import numpy as np
import gtsam
from gtsam import (DiscreteBayesNet, DiscreteFactorGraph, DiscreteMarginals,
DiscreteKeys, Ordering)Model Setup¶
nr_nodes = 4
nr_states = 3
# Define variables as well as ordering
ordering = Ordering()
keys = []
for k in range(nr_nodes):
key_i = (k, nr_states)
keys.append(key_i)
ordering.push_back(k)Create HMM as a DiscreteBayesNet¶
hmm = DiscreteBayesNet()
# Define backbone transition model
# "8/1/1 1/8/1 1/1/8" means the transition matrix favors staying in the same state
# For each parent state (0,1,2), it gives probabilities for child states (0,1,2)
transition = "8/1/1 1/8/1 1/1/8"
# Add transition conditionals for k=1,2,3 conditioned on k-1
for k in range(1, nr_nodes):
# Helper function to create DiscreteKeys for parents
parents = DiscreteKeys()
parents.push_back(keys[k - 1])
hmm.add(keys[k], parents, transition)# Add some measurements, not needed for all time steps!
# These are observation likelihoods
hmm.add(keys[0], "7/2/1") # Measurement at time 0
hmm.add(keys[1], "1/9/0") # Measurement at time 1
hmm.add(keys[3], "5/4/1") # Measurement at time 3 (last node)
# Print the HMM
hmm.print("HMM\n")HMM
size: 6
conditional 0: P( 1 | 0 ):
Choice(1)
0 Choice(0)
0 0 Leaf 0.8
0 1 Leaf 0.1
0 2 Leaf 0.1
1 Choice(0)
1 0 Leaf 0.1
1 1 Leaf 0.8
1 2 Leaf 0.1
2 Choice(0)
2 0 Leaf 0.1
2 1 Leaf 0.1
2 2 Leaf 0.8
conditional 1: P( 2 | 1 ):
Choice(2)
0 Choice(1)
0 0 Leaf 0.8
0 1 Leaf 0.1
0 2 Leaf 0.1
1 Choice(1)
1 0 Leaf 0.1
1 1 Leaf 0.8
1 2 Leaf 0.1
2 Choice(1)
2 0 Leaf 0.1
2 1 Leaf 0.1
2 2 Leaf 0.8
conditional 2: P( 3 | 2 ):
Choice(3)
0 Choice(2)
0 0 Leaf 0.8
0 1 Leaf 0.1
0 2 Leaf 0.1
1 Choice(2)
1 0 Leaf 0.1
1 1 Leaf 0.8
1 2 Leaf 0.1
2 Choice(2)
2 0 Leaf 0.1
2 1 Leaf 0.1
2 2 Leaf 0.8
conditional 3: P( 0 ):
Choice(0)
0 Leaf 0.7
1 Leaf 0.2
2 Leaf 0.1
conditional 4: P( 1 ):
Choice(1)
0 Leaf 0.1
1 Leaf 0.9
2 Leaf 0
conditional 5: P( 3 ):
Choice(3)
0 Leaf 0.5
1 Leaf 0.4
2 Leaf 0.1
Inference: Maximum Probable Explanation (MPE)¶
# Convert to factor graph
factor_graph = DiscreteFactorGraph(hmm)
# Do max-product to find the most probable explanation (MPE)
mpe = factor_graph.optimize()
print("mpe:", end="")
for k in range(nr_nodes):
print(f" ({k}, {mpe[k]})", end="")
print()
mpe: (0, 1) (1, 1) (2, 1) (3, 1)
Elimination, Sampling¶
# Create solver and eliminate
# This will create a DAG ordered with arrow of time reversed
chordal = factor_graph.eliminateSequential(ordering)
chordal.print("Eliminated\n")
# We can also sample from it
print("\n10 samples:")
for k in range(10):
sample = chordal.sample()
print("sample:", end="")
for i in range(nr_nodes):
print(f" ({i}, {sample[i]})", end="")
print()Eliminated
size: 4
conditional 0: P( 0 | 1 ):
Choice(1)
0 Choice(0)
0 0 Leaf 0.94915254
0 1 Leaf 0.033898305
0 2 Leaf 0.016949153
1 Choice(0)
1 0 Leaf 0.29166667
1 1 Leaf 0.66666667
1 2 Leaf 0.041666667
2 Choice(0)
2 0 Leaf 0.41176471
2 1 Leaf 0.11764706
2 2 Leaf 0.47058824
conditional 1: P( 1 | 2 ):
Choice(2)
0 Choice(1)
0 0 Leaf 0.68604651
0 1 Leaf 0.31395349
0 2 Leaf 0
1 Choice(1)
1 0 Leaf 0.033016228
1 1 Leaf 0.96698377
1 2 Leaf 0
2 Choice(1)
2 0 Leaf 0.21454545
2 1 Leaf 0.78545455
2 2 Leaf 0
conditional 2: P( 2 | 3 ):
Choice(3)
0 Choice(2)
0 0 Leaf 0.72746497
0 1 Leaf 0.23618821
0 2 Leaf 0.036346815
1 Choice(2)
1 0 Leaf 0.045088145
1 1 Leaf 0.9368897
1 2 Leaf 0.018022151
2 Choice(2)
2 0 Leaf 0.14716578
2 1 Leaf 0.38224599
2 2 Leaf 0.47058824
conditional 3: P( 3 ):
Choice(3)
0 Leaf 0.36536251
1 Leaf 0.58948629
2 Leaf 0.045151196
10 samples:
sample: (0, 0) (1, 1) (2, 1) (3, 1)
sample: (0, 1) (1, 1) (2, 1) (3, 1)
sample: (0, 0) (1, 0) (2, 0) (3, 0)
sample: (0, 1) (1, 1) (2, 1) (3, 1)
sample: (0, 0) (1, 1) (2, 1) (3, 1)
sample: (0, 0) (1, 1) (2, 1) (3, 1)
sample: (0, 0) (1, 0) (2, 0) (3, 0)
sample: (0, 1) (1, 1) (2, 1) (3, 1)
sample: (0, 1) (1, 1) (2, 1) (3, 1)
sample: (0, 0) (1, 1) (2, 1) (3, 1)
# Or compute the marginals. This re-eliminates the FG into a Bayes tree
print("\nComputing Node Marginals ..")
marginals = DiscreteMarginals(factor_graph)
for k in range(nr_nodes):
marg_probs = marginals.marginalProbabilities((k, nr_states))
print(f"marginal {k}{marg_probs}")
Computing Node Marginals ..
marginal 0[0.44714654 0.5170319 0.03582156]
marginal 1[0.23647637 0.76352363 0. ]
marginal 2[0.29901199 0.65583682 0.0451512 ]
marginal 3[0.36536251 0.58948629 0.0451512 ]