Rabiner 1989 - Problem 2 - State Posteriors and the Viterbi Algorithm

Summary

Two decoders answer Problem 2. Posterior (marginal) decoding picks $q_t=\arg {\max }_i{\gamma }_t(i)$ with ${\gamma }_t(i)={\alpha }_t(i){\beta }_t(i)/P(O\mid \lambda )$, maximising the expected number of correct states but possibly producing impossible sequences. Viterbi dynamic programming finds the single best path maximising $P(Q,O\mid \lambda )$ with ${\delta }_{t+1}(j)=[{\max }_i{\delta }_t(i)a_{ij}]b_j(O_{t+1})$ plus backpointers.

Overview

“There are several possible ways of solving Problem 2”: the difficulty is the optimality criterion. Different criteria give different state sequences.

Main Content

State posterior $\gamma$ (eqs 26-28)

${\gamma }_t(i)=P(q_t=S_i\mid O,\lambda )=\frac{{\alpha }_t(i){\beta }_t(i)}{P(O\mid \lambda )}=\frac{{\alpha }_t(i){\beta }_t(i)}{{\sum }_{i=1}^N{\alpha }_t(i){\beta }_t(i)}$ with ${\sum }_i{\gamma }_t(i)=1$. ${\alpha }_t$ accounts for $O_1..O_t$ and state $S_i$ at $t$; ${\beta }_t$ for $O_{t+1}..O_T$ given state $S_i$.

Individually most likely states (eq. 29)

$q_t=\arg {\max }_{1\le i\le N}{\gamma }_t(i)$ maximises the expected number of correctly decoded individual states. Caveat (paper): it ignores the probability of sequences of states; if some $a_{ij}=0$ the resulting sequence may not even be a valid path. Alternatives: maximise expected correct pairs $(q_t,q_{t+1})$ or triples.

Viterbi algorithm (eqs 30-35)

Define ${\delta }_t(i)={\max }_{q_1..q_{t-1}}P[q_1..q_{t-1},q_t=i,O_1..O_t\mid \lambda ]$ (eq. 30), the best score along a single path ending in state $i$. Induction: ${\delta }_{t+1}(j)=[{\max }_i{\delta }_t(i)a_{ij}]b_j(O_{t+1})$ (eq. 31).

1. Initialisation: ${\delta }_1(i)={\pi }_i{b}_i(O_1)$, ${\psi }_1(i)=0$.
2. Recursion ($2\le t\le T$, $1\le j\le N$): ${\delta }_t(j)={\max }_i[{\delta }_{t-1}(i)a_{ij}]b_j(O_t)$, ${\psi }_t(j)=\arg {\max }_i[{\delta }_{t-1}(i)a_{ij}]$.
3. Termination: $P^{\ast }={\max }_i{\delta }_T(i)$, $q_T^{\ast }=\arg {\max }_i{\delta }_T(i)$.
4. Backtracking: $q_t^{\ast }={\psi }_{t+1}(q_{t+1}^{\ast })$, $t=T-1,\dots ,1$. It equals the forward calculation with $\max$ replacing $\sum$ plus backtracking; same lattice. In logs (eqs 104-105): ${\varphi }_t(j)={\max }_i[{\varphi }_{t-1}(i)+\log a_{ij}]+\log b_j(O_t)$, no scaling needed, and $\log a_{ij}$, $\log b_j(O_t)$ can be precomputed.

Why the two decoders differ (illustration)

Posterior decoding can return a sequence containing a zero-probability transition (the paper’s caveat), whereas Viterbi never does, because Viterbi maximises the joint path probability. Constructing such a case needs at least three states: with two states and $a_{12}=0$, state 1 can only be followed by state 1, so ${\gamma }_{t+1}(1)\ge {\gamma }_t(1)$ and the decoded pair $(S_1,S_2)$ cannot occur. In the simulated 2-state Poisson HMM (200 steps, rates 1 and 6, $A=\left(\begin{array}{cc}.9& .1\\ .2& .8\end{array}\right)$) Viterbi was 99.5% accurate and posterior decoding 99.0%: close here because states are well separated.

import numpy as np
def viterbi(logB, logA, logpi):
    T, N = logB.shape; d = np.zeros((T, N)); psi = np.zeros((T, N), int); d[0] = logpi + logB[0]
    for t in range(1, T):
        s = d[t-1][:, None] + logA; psi[t] = s.argmax(0); d[t] = s.max(0) + logB[t]   # eqs 33a-b in logs
    q = np.zeros(T, int); q[-1] = d[-1].argmax()
    for t in range(T-2, -1, -1): q[t] = psi[t+1][q[t+1]]                              # eq. 35
    return q, d[-1].max()
# posterior decoding: gamma = exp(la + lb - logsumexp(la[-1])); q = gamma.argmax(1)

Examples (PyMC-style)

E24 runs forward-backward and Viterbi per posterior draw to get regime probabilities and a decoded path with uncertainty (the Viterbi path varies across draws of $(A,B)$; posterior state probabilities average over them). Untested sketch: for each draw $s$ in idata.posterior, compute la, lb from the sampled $A^{(s)},{\lambda }^{(s)}$ as above, then average ${\gamma }^{(s)}$ across draws.

Connections

Used in Problem 3 as the segmentation step in segmental $k$-means (Rabiner 1989 - Implementation Issues - Scaling, Multiple Sequences and Initialisation). Movement-state decoding: Model Building with Latent Variables - Animal Movement.

See Also

Rabiner 1989 - Problem 1 - Forward and Backward Algorithms, Rabiner 1989 - Problem 3 - Baum-Welch Reestimation