Inferring hidden states from what you can observe.
On this page
A Markov chain is a sequence of states where each step picks the next state by fixed chances. Hide the state, show only what it emits, and you have an HMM, the model behind decades of speech recognition.
A hidden Markov model (HMM) is a Markov chain whose state you never see, only what each state emits. Given a sequence it answers three questions, as Rabiner (1989) organizes them: how likely the sequence is (the forward pass), which hidden state made each step (Viterbi), and which chances make the sequence most likely (Baum-Welch, which is EM). Baum-Welch fits once, offline, over the training sequences; after that, scoring or decoding a new sequence is one pass over it.
Take a toy weather model with made-up numbers. A sunny day stays sunny 80% of the time; a rainy day stays rainy 60% of the time. Nothing before today matters, which is the Markov assumption (the independence idea from Probability and Bayes Rule).
Figure 3. Stickier weather, longer runs; the rainy share stays a third.
This snippet simulates the chain from A alone and prints the share of rainy days; the rainy share comes out near one third.
import numpy as np
SUNNY, RAINY = 0, 1
A = np.array([[0.8, 0.2], [0.4, 0.6]])
rng = np.random.default_rng(1)
days = [SUNNY]
for _ in range(999):
days.append(rng.choice(2, p=A[days[-1]]))
print(np.mean(np.array(days) == RAINY)) # about 0.33Hide which die is in use, and you get an HMM. A dealer mostly rolls a fair die but sometimes swaps in a loaded one that rolls six half the time (the occasionally dishonest casino, Durbin et al. 1998). You see the rolls, never the die.
Adding up every possible die sequence is hopeless: 40 rolls have 240, about 1012, of them.
CTC, the speech recognition loss that adds up every way to line a transcript up with the audio, fills the same kind of table over label positions. The totals shrink by roughly a factor of 6 per roll, so real code works in logs (see Maximum Likelihood and EM) or rescales each column.
This snippet runs the forward pass in logs on nine rolls and prints the log chance of all of them.
import numpy as np
from scipy.special import logsumexp
logA = np.log([[0.95, 0.05], [0.10, 0.90]])
logB = np.log([[1/6] * 6, [0.1] * 5 + [0.5]])
pi = np.array([0.5, 0.5])
o = np.array([3, 1, 4, 2, 5, 3, 6, 6, 6]) - 1
log_a = np.log(pi) + logB[:, o[0]]
for t in range(1, len(o)):
log_a = logsumexp(log_a[:, None] + logA, axis=0) + logB[:, o[t]]
print(logsumexp(log_a)) # log P(rolls), about -15.41The Viterbi algorithm (Viterbi 1967) answers the second question, which die made each roll. It fills the same table as the forward pass but keeps only the best path into each cell, a max in place of the sum, and remembers which cell that path came from. At the last roll it walks those remembered choices back and reads off one die per roll.
Figure 7. Darker cell, better die so far; three sixes flip the path at roll 7.
This snippet runs Viterbi on the same nine rolls and prints the best die per roll.
import numpy as np
logA = np.log([[0.95, 0.05], [0.10, 0.90]])
logB = np.log([[1/6] * 6, [0.1] * 5 + [0.5]])
o = np.array([3, 1, 4, 2, 5, 3, 6, 6, 6]) - 1
log_delta = np.log([0.5, 0.5]) + logB[:, o[0]]
back = []
for t in range(1, len(o)):
scores = log_delta[:, None] + logA
back.append(scores.argmax(axis=0))
log_delta = scores.max(axis=0) + logB[:, o[t]]
path = [int(log_delta.argmax())]
for ptr in reversed(back):
path.append(int(ptr[path[-1]]))
print(path[::-1]) # [0, 0, 0, 0, 0, 0, 1, 1, 1], 0 = fairViterbi gives one best path. Sometimes you want a chance for every roll instead: how likely is it that the loaded die made roll t? Run the table backward as well and combine the two.
On 30 rolls, the model’s switch chance matters. Fair goes to loaded with chance q and back with 2q; the dealer’s real q is 0.05.
Figure 8. A low q almost never switches dice. At q = 0.3 the path switches 5 to 12 times in 30 rolls, jumping to loaded at lone sixes.
If A and B are unknown, learn them from the rolls with EM, Expectation-Maximization, here called Baum-Welch (Baum et al. 1970). The E step runs forward and backward to get soft counts; the M step divides counts. The updates are essentially EM’s, so the likelihood never drops and the fit can stop at a local peak.
The official hmmlearn casino example fits 50 times from random starts and keeps the best fit.
With hmmlearn the three questions are three calls on rolls stored as one column, and the fit runs from several starts as Section 10 advises; several sequences are stacked into one array, and lengths says where each one ends, so no step is counted across the join.
import numpy as np
from hmmlearn import hmm
casino = hmm.CategoricalHMM(n_components=2)
casino.startprob_ = np.array([0.5, 0.5])
casino.transmat_ = np.array([[0.95, 0.05], [0.10, 0.90]])
casino.emissionprob_ = np.array([[1/6] * 6, [0.1] * 5 + [0.5]])
rs = np.random.RandomState(1)
X1, _ = casino.sample(300, random_state=rs)
X2, _ = casino.sample(200, random_state=rs)
X, lengths = np.concatenate([X1, X2]), [len(X1), len(X2)]
best_score, model = -np.inf, None
for seed in range(10):
candidate = hmm.CategoricalHMM(n_components=2, n_iter=50, random_state=seed)
candidate.fit(X, lengths)
score = candidate.score(X, lengths)
if score > best_score:
best_score, model = score, candidate
log_prob, path = model.decode(X2)
gamma = model.predict_proba(X2)
print(model.transmat_.round(2), log_prob)rs and records their lengths.| step | when, on what, and its cost |
|---|---|
fit | Once, in Baum-Welch rounds; each round runs the forward and backward pass over every sequence. |
decode, predict_proba, score | Per new sequence, one Viterbi pass, or one forward and backward pass. |
| state | π, A and B (startprob_, transmat_, emissionprob_). |
| setting | raise it, lower it | default |
|---|---|---|
n_components, the hidden states j in A and B | More states fit finer patterns and soon overfit; 1 state means no hidden switching at all, one die for every roll. | 1 |
n_iter, tol, the Baum-Welch rounds | Fit stops after n_iter rounds or when a round gains less than tol in total log-likelihood. A smaller tol runs more rounds and lands closer to the peak; too few rounds stop short of it. | 10, 0.01 |
init_params, the starting π, A and B | Says which of π, A and B (letters s, t, e) fit sets fresh before the first Baum-Welch round, so it picks where the climb starts. With all three letters the random_state seed decides the start, so 10 seeds give 10 different starts. Fewer letters keep the values you set yourself; an empty string starts every seed from your own π, A and B. | "ste" |
lengths | Without it, separate sequences are read as one chain with a false step at every join. | one sequence |
Defaults are from the hmmlearn API reference. Gotchas are in Section 10.
Classic speech recognizers model each word or sound unit as a chain of left-to-right states that can stay or move on, never back. A Gaussian mixture under each state scores the sound frames, often MFCCs (see MFCC and Other Audio Features), and Viterbi finds the best path. Hinton et al. (2012) showed deep neural networks (see Multilayer Perceptrons) beating those mixtures as the frame scorers on many benchmarks, and CTC later dropped the explicit HMM while keeping a forward table.
Mixtures and HMMs look for hidden groups. Principal Component Analysis looks for the few directions along which the data stretches most.
hmmlearn (a Python library for HMMs) documentation, Dishonest Casino Example and API reference. The page’s switch and loaded-die numbers are the ones in this example and in PMTK, both citing Durbin et al. ch. 3; the even start is PMTK’s.