Hidden Markov models: the Viterbi algorithm and the dishonest casino

Let N, M ≥ 1 be integers, S = {1, …, N} the set of hidden states and V = {1, …, M} the set of observable symbols. Further let

• π = (π₁, …, π_N) ∈ ℝ^1×N be a probability distribution on S (π_i ≥ 0, ∑_i π_i = 1) — the initial distribution;
• A = (a_ij) ∈ ℝ^N×N a stochastic matrix — the transition matrix;
• B = (b_i(o))_{i ∈ S, o ∈ V} ∈ ℝ^N×M a stochastic matrix — the emission probabilities.

A hidden Markov model (HMM) with parameters λ = (π, A, B) is a pair of sequences of random variables X₁, X₂, … : Ω → S (hidden states) and O₁, O₂, … : Ω → V (observations) on a probability space (Ω, 𝓕, P) such that:

1. (X_t)_{t ≥ 1} is a homogeneous Markov chain with initial distribution P(X₁ = i) = π_i and transition matrix A, i.e. a_ij = P(X_t+1 = j | X_t = i);
2. for every T ≥ 1, t ∈ {1, …, T}, o ∈ V and all values x₁, …, x_T ∈ S and o_s ∈ V (s ≠ t) for which the condition has positive probability, P(O_t = o | X₁ = x₁, …, X_T = x_T, O_s = o_s for s ∈ {1, …, T} ∖ {t}) = b_xt(o).

In words: the hidden states change like an ordinary Markov chain, and the observation at time t depends only on the state X_t — not on the other states or on the other observations. We observe only O₁, O₂, …, not the states X_t.

Our sequence begins with the rolls 1, 6, 3. Let F denote the fair die and L the loaded one. The history “loaded three times”:

P(LLL, 1, 6, 3) = 0.5 · 0.1 · 0.9 · 0.5 · 0.9 · 0.1 = 0.002025,

the history “fair three times”:

P(FFF, 1, 6, 3) = 0.5 · 1/6 · 0.95 · 1/6 · 0.95 · 1/6 = 0.002089

(colours as in the formula above: initial distribution, emissions, transitions). The fair history comes out a whisker more probable: the six on the loaded die does get 0.5 instead of 1/6, but the one and the three get only 0.1 instead of 1/6. (In the simulation, all three rolls were in fact made with the loaded die — three rolls are not enough to tell.)

The probability of the observations is the sum of the joint probability over all possible histories: P(O | λ) = ∑_x P(x, O | λ). For the rolls 1, 6, 3 there are 2³ = 8 of them:

Each history costs 2T − 1 = 5 multiplications, 8 · 5 = 40 in total. For T rolls it is (2T − 1) · 2^T — for 300 rolls an unimaginable number.

Input: a hidden Markov model λ = (π, A, B) with states S = {1, …, N} and symbols V; an observed sequence O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Output: P(O | λ) = P(O₁, …, O_T) and the forward variables

α_t(j) = P(O₁, …, O_t, X_t = j | λ),   t = 1, …, T,   j ∈ S.

1. Initialization: for every j ∈ S set α₁(j) = π_j b_j(O₁).
2. Induction: for t = 1, …, T − 1 and every j ∈ S set α_t+1(j) = [∑_i=1^N α_t(i) a_ij] · b_j(O_t+1).
3. Termination: after processing time T, return P(O | λ) = ∑_j=1^N α_T(j).

The algorithm stops after T − 1 induction steps; each costs O(N²) operations, so O(N² T) in total instead of O(T · N^T) when summing over all histories. To keep the numbers from underflowing, at every time t we compute c_t = ∑_j α_t(j) (with the already normalized α_t−1 in the induction), replace α_t(j) by α_t(j)/c_t and return log P(O | λ) = ∑_t=1^T log c_t.

Start: α₁(F) = 0.5 · 1/6 = 0.0833, α₁(L) = 0.5 · 0.1 = 0.0500. The second roll is a six (colours as in the formula above):

α₂(F) = (0.0833 · 0.95 + 0.0500 · 0.1) · 1/6 = 0.0140,

α₂(L) = (0.0833 · 0.05 + 0.0500 · 0.9) · 0.5 = 0.0246.

The sum in the last row is exactly the sum of the eight histories from the previous example, only instead of 40 multiplications 2 + 2 · 6 = 14 were enough. The last column c_t = P(O_t | O₁, …, O_t−1) is the normalizing constant from the scaling: how probable roll t was given the previous ones. Their product is P(O | λ), so log P(O | λ) = ∑_t log c_t.

Why logarithms? For 3,000 rolls P(O | λ) comes out around 10^−2,278, and the smallest positive number a computer can store in ordinary double precision (IEEE 754) is about 4.9 · 10⁻³²⁴. Without scaling, α would underflow to zero.

Let (Ω, 𝓕, P) be a probability space, M₁, M₂ ∈ 𝓕 the events “model 1 holds” and “model 2 holds” with P(M₁) > 0, P(M₂) > 0, and O ∈ 𝓕 the event “we observed the data” with P(O | M₂) > 0. The Bayes factor of model M₁ against M₂ is the number

K = P(O | M₁)/P(O | M₂),

i.e. the ratio of the likelihoods of the data under the two models. K > 1 means that the data are more probable under the model M₁. If moreover P(O) > 0, then Bayes’ theorem written for both models and divided (P(O) cancels out) gives:

Calculation. log₁₀ K = −230.03 − (−233.45) = 3.42 (3.413 from the unrounded values), so K ≈ 10^3.413 ≈ 2,590, roughly 2,600. If you considered cheating and a fair game equally likely beforehand (prior odds 1), after 300 rolls the odds are 2,600 · 1 = 2,600 : 1 in favour of cheating. On the widely used scale of Kass and Raftery (1995), a Bayes factor above 150 is “very strong” evidence.

Input: a hidden Markov model λ = (π, A, B) with states S = {1, …, N}; observations O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Output: the most probable sequence of hidden states and its probability

x* = argmax_{x ∈ S^T} P(X₁ = x₁, …, X_T = x_T, O | λ),   P* = P(x*, O | λ).

For t = 1, …, T and j ∈ S let δ_t(j) be the probability of the best path that ends in state j at time t:

δ_t(j) = max_{x1, …, xt−1 ∈ S} P(X₁ = x₁, …, X_t−1 = x_t−1, X_t = j, O₁, …, O_t | λ).

1. Initialization: for every j ∈ S set δ₁(j) = π_j b_j(O₁).
2. Induction: for t = 1, …, T − 1 and every j ∈ S set δ_t+1(j) = max_{i ∈ S} [δ_t(i) a_ij] · b_j(O_t+1),   ψ_t+1(j) = argmax_{i ∈ S} δ_t(i) a_ij (ψ_t+1(j) is the best predecessor; on a tie choose any).
3. Termination: P* = max_{j ∈ S} δ_T(j) and x*_T = argmax_{j ∈ S} δ_T(j).
4. Backtracking: for t = T − 1, T − 2, …, 1 set x*_t = ψ_t+1(x*_t+1); after t = 1 return x* = (x*₁, …, x*_T) and P*.

The algorithm stops after T − 1 induction steps and T − 1 backtracking steps; the time is again O(N² T). In practice logarithms are added instead of multiplying products, so that the numbers do not underflow.

The start is the same as for α: δ₁(F) = 0.0833, δ₁(L) = 0.0500. For the six on the second roll:

δ₂(F) = max(0.0833 · 0.95, 0.0500 · 0.1) · 1/6 = 0.0792 · 1/6 = 0.0132   (F wins),

δ₂(L) = max(0.0833 · 0.05, 0.0500 · 0.9) · 0.5 = 0.0450 · 0.5 = 0.0225   (L wins)

(colours as in the formula above).

At the end F wins (0.002089 > 0.002025), and following the pointers ψ backwards gives the path FFF — exactly the history from the first example that was the most probable of all eight. The forward algorithm, however, gives for the third roll P(L | O₁, O₂, O₃) = α₃(L)/α₃(F) + α₃(L) = 0.002283/0.004913 = 0.465. The best path and the most probable state agree here (F), but the “victory” is narrow and three rolls are not enough to be sure. In general they need not agree — we will come back to that with posterior decoding.

Let λ = (π, A, B) be a hidden Markov model with states S = {1, …, N} and O = (O₁, …, O_T) observations with P(O | λ) > 0. The posterior probability of the state i ∈ S at time t ∈ {1, …, T} is

γ_t(i) = P(X_t = i | O₁, …, O_T, λ).

It is computed with the backward variables (for t < T and P(X_t = i) > 0)

β_t(i) = P(O_t+1, …, O_T | X_t = i, λ),   β_T(i) = 1.

The forward–backward algorithm computes α_t(i) by the forward algorithm, β_t(i) from the end (for t = T − 1, …, 1) and from both γ_t(i) by the formulas below. In words: γ_t(i) is the probability that the model was in state i at time t, given the whole observed sequence — past and future.

Let us stay with the rolls 1, 6, 3 and ask about the second roll (the six). The backward step from β₃ = 1 with the third roll (the three):

β₂(F) = 0.95 · 1/6 + 0.05 · 0.1 = 0.1633,   β₂(L) = 0.1 · 1/6 + 0.9 · 0.1 = 0.1067.

(In each product the first factor is the transition a_ij, the second the emission of the three b_j(3), and β₃(j) = 1 is not written.) By the formula above (colours match its terms):

γ₂(L) = 0.024583 · 0.106667/0.004913 = 0.534.

If we knew only the first two rolls, we would have P(L | O₁, O₂) = α₂(L)/α₂(F) + α₂(L) = 0.637. The three that came up afterwards is less probable on the loaded die, so it retroactively lowers the suspicion about the second roll to 0.534.

Let λ be a hidden Markov model with states S, O = (O₁, …, O_T) observations and γ_t(i) the posterior probabilities of the states from the previous definition. Posterior decoding is the sequence x̂ = (x̂₁, …, x̂_T) ∈ S^T, where

x̂_t = argmax_{i ∈ S} γ_t(i)   for every t = 1, …, T

(the state is chosen at each time separately). Further let x ∈ S^T be the true sequence of states and y ∈ S^T any decoding (posterior or Viterbi). The accuracy of the decoding y is the fraction of times at which it agrees with reality:

accuracy(y, x) = |{t ∈ {1, …, T} : y_t = x_t}|/T ∈ [0, 1].

Posterior decoding maximizes the expected number of correctly identified rolls, Viterbi the probability that the whole path is correct at once. That is why a posterior sequence may even contain a transition that has probability zero in the model — the casino has no such transition, but in general it can happen.

It is an EM algorithm (Baum et al., around 1970).

Input: the number of hidden states N ≥ 1, the set of symbols V, observations O = (O₁, …, O_T) ∈ V^T with T ≥ 2, an initial estimate of the parameters λ⁽⁰⁾ = (π, A, B) with P(O | λ⁽⁰⁾) > 0 and a number of iterations K ≥ 1 (or a tolerance ε > 0).
Output: parameters λ^(K) with likelihood P(O | λ^(K)) ≥ P(O | λ⁽⁰⁾).

For k = 0, 1, …, K − 1 repeat with the parameters λ^(k):

1. E step: use the forward–backward algorithm to compute α_t(i), β_t(i), P(O | λ^(k)), the posterior probabilities γ_t(i) for t = 1, …, T and the transition probabilities ξ_t(i, j) = P(X_t = i, X_t+1 = j | O, λ^(k)),   t = 1, …, T − 1,   i, j ∈ S.
2. M step: the new parameters λ^(k+1) are, for all i, j ∈ S and o ∈ V, π_i = γ₁(i),   a_ij = ∑_t=1^T−1 ξ_t(i, j)/∑_t=1^T−1 γ_t(i),   b_i(o) = ∑_{t : Ot = o} γ_t(i)/∑_t=1^T γ_t(i).
3. Stopping: stop after K iterations, or earlier if log P(O | λ^(k+1)) − log P(O | λ^(k)) < ε.

P(O | λ^(k+1)) ≥ P(O | λ^(k)) holds — the likelihood never decreases in any iteration. The algorithm converges to a stationary point of the likelihood, in practice usually a local maximum that need not be global; the result therefore depends on the initial estimate.