Hidden Markov Models: Forward, Viterbi and Baum-Welch, Built From Scratch

Every model up to this point mapped a row of features to an answer. Logistic regression read a feature vector and returned a class. K-means grouped independent points. Even a Gaussian mixture treated each sample as a fresh draw with no memory of the last one.

A Hidden Markov Model assumes the opposite: your observations were produced by a hidden sequence of states that you cannot see, that evolves in time, and where what happens next depends on where you are now.

A casino quietly swaps a fair die for a loaded one and swaps back — you only ever see the rolls. A phone hears an acoustic waveform — the words are hidden. A sequencer reads bases — whether you are inside a coding region is hidden.

Three tables and two assumptions

An HMM is λ = (π, A, B):

• π[i] — probability the chain starts in state i
• A[i][j] — probability of moving from state i to state j
• B[i][k] — probability that state i emits symbol k

Every row sums to 1. Two assumptions make everything computable. The Markov property: the next state depends only on the current one, not the whole past. Output independence: an observation depends only on the state that emitted it.

model = {
  "pi": [0.85, 0.15],
  "A":  [[0.95, 0.05],          # FAIR   -> FAIR / LOADED
         [0.10, 0.90]],         # LOADED -> FAIR / LOADED
  "B":  [[1/6]*6,               # a fair die
         [.1,.1,.1,.1,.1,.5]]   # a die that loves sixes
}

You observe the rolls. You never observe FAIR vs LOADED.

The wall

P(O | λ) has a clear definition: sum the joint probability over every hidden path. Each path is a simple product. The problem is how many there are — N^T.

T=9,   N=2  ->            512 paths
T=20,  N=2  ->      1,048,576 paths
T=100, N=2  ->  1.3 × 10^30 paths

At a billion paths per second, that last sum finishes long after the sun burns out. What rescues HMMs is not an approximation. It is the observation that those paths are enormously redundant: vast numbers of them share the same prefix, and a shared prefix contributes the same factor to every path that shares it.

Forward — one column at a time

Define αₜ(i) = P(o₀…oₜ, qₜ = i): the total probability of everything seen so far and being in state i right now. The Markov property means this is enough to keep going — you need the accumulated mass in each state, not the paths that produced it.

alpha = [[pi[i] * B[i][o[0]] for i in range(N)]]
for t in range(1, T):
    col = []
    for j in range(N):
        hop = sum(alpha[t-1][i] * A[i][j] for i in range(N))
        col.append(hop * B[j][o[t]])      # gather, then emit
    alpha.append(col)
P_O = sum(alpha[T-1])

T·N² multiply-adds instead of N^T products. Exact, not approximate. The exponential did not get approximated away — it got factored away.

The bug that only shows up on long sequences

Written exactly like that, the forward algorithm is mathematically right and practically broken. Every step multiplies by two probabilities below 1, so α shrinks geometrically and underflows to exactly 0 in float64 after a few hundred timesteps. Everything downstream becomes NaN.

The fix is three lines: normalise each column, keep the normaliser.

c = []
for t in range(T):
    ...
    ct = sum(col); c.append(ct)
    col = [v / ct for v in col]
logP = sum(math.log(ct) for ct in c)      # never underflows

The normalised α̂ₜ(i) is worth having on its own — it is P(qₜ = i | o₀…oₜ), the filtered belief a real-time tracker reports. This is the single most common source of silent HMM bugs.

Viterbi — swap the sum for a max

Forward asks how much total probability reaches each cell. Viterbi asks what the best single path reaching it is. Identical recursion, Σ becomes max, and you record which predecessor won so you can walk the winners backwards.

d   = [[log(pi[i]) + log(B[i][o[0]]) for i in range(N)]]
psi = [[0]*N for _ in range(T)]
for t in range(1, T):
    col = []
    for j in range(N):
        best, arg = max((d[t-1][i] + log(A[i][j]), i) for i in range(N))
        col.append(best + log(B[j][o[t]])); psi[t][j] = arg
    d.append(col)

path = [max(range(N), key=lambda i: d[T-1][i])]
for t in range(T-1, 0, -1): path.append(psi[t][path[-1]])
path.reverse()

Do it in log space — sums of logs never underflow, and since log is monotone the argmax is untouched.

Viterbi returns a globally consistent path: every transition in it is one the model permits. The obvious alternative — picking the most likely state independently at each timestep — can stitch together a sequence containing a zero-probability transition, a path the model considers literally impossible.

Backward, and the posteriors

The filtered belief only knows the past. Once the whole sequence is in hand, evidence after time t also tells you about the state at time t. βₜ(i) captures that, computed by the same recursion run right to left. Multiply the two and you get γₜ(i) = P(qₜ = i | the entire sequence) — the smoothed posterior, always at least as sharp as the filtered one.

beta[T-1] = [1/c[T-1]] * N                     # reuse forward's scale factors
for t in range(T-2, -1, -1):
    for i in range(N):
        beta[t][i] = sum(A[i][j]*B[j][o[t+1]]*beta[t+1][j]
                         for j in range(N)) / c[t]

gamma = [[alpha[t][i]*beta[t][i]*c[t] for i in range(N)] for t in range(T)]

Filtering is what you can do online. Smoothing is what you can do afterwards, and it is strictly better informed.

Baum-Welch — learning with no labels at all

If the hidden states were visible, fitting the model would be counting: how often did state i start, how often did i lead to j, how often did i emit k. You cannot see them, so use expected counts under the current model. That is EM applied to sequences.

The E-step is the forward-backward pass you already wrote — γ and ξ are those expected counts. The M-step normalises them into new tables.

for iteration in range(200):
    gamma, xi = forward_backward(o, pi, A, B)          # E-step
    pi = gamma[0]                                      # M-step
    for i in range(N):
        denA = sum(gamma[t][i] for t in range(T-1))
        for j in range(N):
            A[i][j] = sum(xi[t][i][j] for t in range(T-1)) / denA
        denB = sum(gamma[t][i] for t in range(T))
        for k in range(M):
            B[i][k] = sum(gamma[t][i] for t in range(T) if o[t]==k) / denB

Every parameter is learned from observations alone. No labelled states anywhere.

The guarantee, and what it does not promise

Baum-Welch comes with a real theorem: each iteration produces parameters whose likelihood is greater than or equal to the previous ones. Never worse. No learning rate, no divergence — and a decreasing likelihood in your log is proof of a bug, which makes it unusually easy to debug.

What it does not promise is the global maximum, or the true parameters. EM climbs to whatever local optimum sits under its starting point, so real code runs several random restarts and keeps the best likelihood. It also converges to the maximum-likelihood fit of the sample you gave it, which on finite data is deliberately not the process that generated it — you will routinely see EM finish slightly above the true model's score. And because it optimises fit rather than labels, state numbering is arbitrary: a refit's state 0 may be your state 1.

Test it against the definition

Dynamic programming is easy to write and easy to get subtly wrong. An off-by-one in the backward recursion or a misplaced scale factor still produces plausible-looking numbers. The cure is an independent baseline — and for small T the brute-force enumeration is that baseline, since it is the definition coded literally.

T = 9                                  # 2^9 = 512 paths, enumerable
P  = posteriors(obs, pi, A, B)
V  = viterbi(obs, pi, A, B)
BF = brute(obs, pi, A, B)

assert abs(P.logLik  - BF.logLik)      < 1e-10   # forward == sum over paths
assert abs(V.logProb - BF.bestLogProb) < 1e-10   # viterbi == max over paths
assert V.path == BF.bestPath                     # and the SAME path
assert max_abs(P.gamma, BF.gamma)      < 1e-12   # gamma == exact marginals

Run against that baseline on three different HMMs, the implementation behind this walkthrough passes 38/38 assertions — including EM's monotonicity checked over six random initialisations, and parameter recovery within 0.03 of the truth once the sequence is long enough.

Where it belongs today

from hmmlearn.hmm import CategoricalHMM

best = None
for seed in range(10):                       # EM is LOCAL — restart
    m = CategoricalHMM(n_components=2, n_iter=200, random_state=seed)
    m.fit(X, lengths)
    if best is None or m.score(X, lengths) > best.score(X, lengths):
        best = m

best.score(X, lengths)          # the forward algorithm
best.predict(X, lengths)        # the Viterbi path
best.predict_proba(X, lengths)  # the gamma posteriors

Reach for an HMM when your data is a sequence, the quantity you care about is latent, and you want the states to mean something: market regimes, gene annotation, sleep stages, machine health, activity recognition. It shines with little data because it has very few parameters, it trains with no labels, it gives calibrated probabilities rather than a point guess, and you can read its tables and argue with them.

The limits are just as clear. The Markov assumption means the current state must summarise all relevant history, so genuinely long-range structure needs something else. The number of states is a choice you make by BIC or held-out likelihood, never by "it fit better". And O(T·N²) gets painful with many states. That is the trade against an RNN or a Transformer, which drop the assumption and win on scale but need labels or a mountain of data, and hand you an opaque vector instead of a table you can read.

Watch Viterbi recover a hidden FAIR/LOADED track from nothing but dice rolls, step through the forward lattice while the page enumerates all N^T paths beside it, and start Baum-Welch from random parameters: https://dev48v.infy.uk/ml/day59-hidden-markov-models.html