5  Inference: filtering, smoothing, decoding

Given a fully specified HMM and a stretch of data, three questions arise. What is the state now, given everything seen so far? What was the state at time \(t\), given the whole series? And which single sequence of states best explains the data? The three have different answers, computed by different algorithms, and confusing them is a common source of error.

All three algorithms are dynamic programming on a chain, all three cost \(\mathcal{O}(K^2T)\), and all three come out of the factorisation in Theorem 3.1. The treatment follows Zucchini et al. (2016, ch. 4) and Rabiner (1989), with the derivations written out and the numerical issues treated properly rather than deferred.

By the end of this chapter you will be able to

• derive the forward and backward recursions from the model assumptions, and prove they compute what they claim;
• implement either of them without underflow, and explain what the scale factors mean;
• say precisely how Viterbi decoding differs from taking the marginal mode at each time, and construct a case where they disagree;
• explain why the choice of initial distribution stops mattering after a few dozen observations.

NoteThe three problems

	Question	Object	Algorithm
Filtering	state now, given the past	\(\mathbb{P}(S_t = i \mid Y_{1:t})\)	forward
Smoothing	state then, given everything	\(\mathbb{P}(S_t = i \mid Y_{1:T})\)	forward–backward
Decoding	best whole path	\(\arg\max_{s_{1:T}} \mathbb{P}(S_{1:T} = s_{1:T} \mid Y_{1:T})\)	Viterbi

Throughout, \(\theta\) is fixed and known, and we suppress it from the notation. Observations \(y_{1:T}\) are fixed.

5.1 Filtering: the forward algorithm

5.1.1 Forward variables

Definition 5.1 (Forward variables) For \(t = 1,\dots,T\) and \(i \in E\), the unnormalised forward variable is \[ \tilde\alpha_t(i) \;=\; p(S_t = i, \, Y_{1:t} = y_{1:t}) , \] and the filtering probability is \(\alpha_t(i) = \mathbb{P}(S_t = i \mid Y_{1:t} = y_{1:t}) = \tilde\alpha_t(i) / \sum_j \tilde\alpha_t(j)\).

Working with the unnormalised version first is deliberate: it satisfies an exact linear recursion, whereas the normalised version satisfies a ratio of linear recursions. The normalising constants turn out to be useful in their own right.

Theorem 5.1 (Forward recursion) The variables of Definition 5.1 satisfy \[ \tilde\alpha_1(i) = \delta_i f_i(y_1), \qquad \tilde\alpha_{t+1}(j) = f_j(y_{t+1}) \sum_{i=1}^{K} \tilde\alpha_t(i) \, \gamma_{ij} \quad (1 \le t \le T-1), \] and \(L_T(\theta) = \sum_{i=1}^K \tilde\alpha_T(i)\).

Proof. The initial case is the definition of a joint density: \(p(S_1 = i, Y_1 = y_1) = \mathbb{P}(S_1 = i) \, p(y_1 \mid S_1 = i) = \delta_i f_i(y_1)\).

For the induction, marginalise over the state at time \(t\) and factorise: \[ \begin{aligned} \tilde\alpha_{t+1}(j) &= \sum_{i=1}^K p(S_t = i, \, S_{t+1} = j, \, Y_{1:t+1} = y_{1:t+1}) \\ &= \sum_{i=1}^K p(S_t = i, Y_{1:t} = y_{1:t}) \; \mathbb{P}(S_{t+1} = j \mid S_t = i, Y_{1:t} = y_{1:t}) \; p(y_{t+1} \mid S_{t+1} = j, S_t = i, Y_{1:t} = y_{1:t}) . \end{aligned} \] The second factor equals \(\gamma_{ij}\) by assumption (A1) in the form stated in Definition 3.2, and the third equals \(f_j(y_{t+1})\) by assumption (A2). Since \(f_j(y_{t+1})\) does not depend on \(i\) it comes outside the sum, giving the stated recursion.

Summing \(\tilde\alpha_T(i)\) over \(i\) marginalises out \(S_T\) and leaves \(p(y_{1:T}) = L_T(\theta)\).

5.1.2 Matrix form

Collecting \(\tilde\alpha_t = (\tilde\alpha_t(1),\dots,\tilde\alpha_t(K))\) as a row vector and using \(\mathbf{Q}(y)\) from Definition 3.4, the recursion reads \[ \tilde{\boldsymbol{\alpha}}_1 = \boldsymbol{\delta}^\top \mathbf{Q}(y_1), \qquad \tilde{\boldsymbol{\alpha}}_{t+1} = \tilde{\boldsymbol{\alpha}}_t \, \boldsymbol{\Gamma} \, \mathbf{Q}(y_{t+1}), \] so that \(\tilde{\boldsymbol{\alpha}}_T \mathbf{1}\) reproduces Theorem 3.2. The forward algorithm is the left-to-right evaluation of that matrix product; nothing new has been introduced, only an order of operations. Each step costs one vector–matrix multiplication, \(\mathcal{O}(K^2)\), for \(\mathcal{O}(K^2T)\) overall, and only the current vector need be retained if the likelihood is all that is wanted.

5.1.3 Prediction

Filtering also delivers forecasts. The one-step-ahead predictive distribution is \[ p(y_{t+1} \mid y_{1:t}) \;=\; \sum_{j=1}^K \Bigl( \sum_{i=1}^K \alpha_t(i) \gamma_{ij} \Bigr) f_j(y_{t+1}) \;=\; \boldsymbol{\alpha}_t^\top \boldsymbol{\Gamma} \, \bigl( f_j(y_{t+1}) \bigr)_j , \] a mixture of the emission densities with weights given by the predicted state distribution. Iterating \(\boldsymbol{\Gamma}\) gives \(h\)-step forecasts, and by Theorem 2.1 those weights converge to \(\boldsymbol{\pi}\) as \(h\) grows: long-horizon forecasts from an ergodic HMM revert to the stationary mixture, regardless of the current state. Forecast pseudo-residuals built from these predictive distributions are the standard HMM diagnostic in Zucchini et al. (2016, ch. 6).

5.1.4 Numerical stability

The recursion is exact and, implemented naively, useless. Each step multiplies by an emission density and a transition probability, so \(\tilde\alpha_t(i)\) decays roughly geometrically in \(t\). With typical per-step factors of order \(10^{-2}\), double precision underflows to zero somewhere around \(t \approx 150\), and the computed likelihood becomes exactly \(0\) with logarithm \(-\infty\).

Two standard remedies, both exact in principle.

Scaling. Normalise at every step. Set \(c_t = \sum_i \tilde\alpha_t(i)\) computed from the current unnormalised vector, and propagate the normalised vector \(\hat\alpha_t = \tilde\alpha_t / c_t\) instead. Concretely, with \(\hat{\boldsymbol{\alpha}}_0^\top := \boldsymbol{\delta}^\top\), \[ \mathbf{v}_t = \hat{\boldsymbol{\alpha}}_{t-1}^\top \boldsymbol{\Gamma} \mathbf{Q}(y_t), \qquad c_t = \mathbf{v}_t \mathbf{1}, \qquad \hat{\boldsymbol{\alpha}}_t^\top = \mathbf{v}_t / c_t . \]

Proposition 5.1 (Scaling identity) With the above definitions, \(\hat\alpha_t(i) = \mathbb{P}(S_t = i \mid Y_{1:t} = y_{1:t})\), the scale factor is the one-step predictive density \(c_t = p(y_t \mid y_{1:t-1})\), and \[ L_T(\theta) = \prod_{t=1}^{T} c_t, \qquad \ell_T(\theta) = \sum_{t=1}^{T} \log c_t . \]

Proof. By induction, \(\mathbf{v}_t = \tilde{\boldsymbol{\alpha}}_t / \prod_{s<t} c_s\), so \(c_t = \bigl(\sum_i \tilde\alpha_t(i)\bigr) / \prod_{s<t} c_s = p(y_{1:t}) / p(y_{1:t-1})\), which is the predictive density, and the telescoping product of the \(c_t\) is \(p(y_{1:T})\). Normalising \(\mathbf{v}_t\) removes the common factor and leaves the conditional distribution.

The log-likelihood therefore accumulates as a sum of logarithms of quantities of order one. This is the implementation used throughout Zucchini et al. (2016), and it has the pleasant side effect that the quantities you must compute anyway are exactly the predictive densities needed for residual diagnostics.

Log domain. Alternatively, store \(a_t(i) = \log \tilde\alpha_t(i)\) and compute \[ a_{t+1}(j) = \log f_j(y_{t+1}) + \operatorname{logsumexp}_i \bigl\{ a_t(i) + \log \gamma_{ij} \bigr\}, \] using the stable form \(\operatorname{logsumexp}(z) = m + \log \sum_i e^{z_i - m}\) with \(m = \max_i z_i\). This handles a wider dynamic range than scaling and copes with structural zeros in \(\boldsymbol{\Gamma}\) without special cases, at the cost of \(K^2\) exponentials and logarithms per step. Scaling is faster; the log domain is more robust. Chapter 7 compares them in detail.

5.1.5 Filter stability

A practical question: how much does the answer depend on \(\boldsymbol{\delta}\), which is often the least well-determined part of the model?

Theorem 5.2 (Forgetting of the initial distribution) Let \(\alpha_t^{\boldsymbol{\delta}}\) and \(\alpha_t^{\boldsymbol{\delta}'}\) denote the filtering distributions at time \(t\) obtained from the same observations \(y_{1:t}\) but different initial distributions. If the transition matrix satisfies the Doeblin condition \(\gamma_{ij} \ge \epsilon \nu_j\) for some \(\epsilon > 0\) and probability vector \(\boldsymbol{\nu}\), then \[ \bigl\lVert \alpha_t^{\boldsymbol{\delta}} - \alpha_t^{\boldsymbol{\delta}'} \bigr\rVert_{\mathrm{TV}} \;\le\; (1-\epsilon)^{t-1} \, \bigl\lVert \boldsymbol{\delta} - \boldsymbol{\delta}' \bigr\rVert_{\mathrm{TV}} . \]

The proof runs exactly as Proposition 2.7. The filter update is the map \(\mu \mapsto \mu^\top \boldsymbol{\Gamma} \mathbf{Q}(y_t) / (\text{normalisation})\), and the diagonal reweighting by \(\mathbf{Q}(y_t)\) does not increase the Dobrushin coefficient, so the composite map contracts total variation by at least \((1-\epsilon)\) per step. Cappé et al. (2005, ch. 4) and Del Moral and Guionnet (2001) give the general treatment.

The practical reading: two analysts who disagree about \(\boldsymbol{\delta}\) but agree on everything else will produce filtering distributions that are indistinguishable after a few dozen observations, provided the chain mixes. It also explains why the same misspecification barely moves \(\ell_T\) for large \(T\), and hence why \(\boldsymbol{\delta}\) is estimated poorly from a single long series.

5.2 Smoothing: the forward–backward algorithm

Filtering uses only the past. For anything retrospective — state estimates for interpretation, or the expected counts needed by EM — the whole series should be used.

5.2.1 Backward variables

Definition 5.2 (Backward variables) For \(t = 1,\dots,T\) and \(i \in E\), \[ \beta_t(i) \;=\; p\bigl( Y_{t+1:T} = y_{t+1:T} \mid S_t = i \bigr), \] with the convention \(\beta_T(i) = 1\) (an empty product).

Note the asymmetry with Definition 5.1: \(\tilde\alpha_t\) is a joint density and \(\beta_t\) is a conditional one. The asymmetry is what makes their product come out right.

Theorem 5.3 (Backward recursion) For \(t = T-1, T-2, \dots, 1\), \[ \beta_t(i) \;=\; \sum_{j=1}^{K} \gamma_{ij} \, f_j(y_{t+1}) \, \beta_{t+1}(j), \qquad\text{i.e.}\qquad \boldsymbol{\beta}_t = \boldsymbol{\Gamma}\mathbf{Q}(y_{t+1}) \boldsymbol{\beta}_{t+1} \] with \(\boldsymbol{\beta}_t\) a column vector.

Proof. Condition on the next state and split the future observations at time \(t+1\): \[ \beta_t(i) = \sum_j \mathbb{P}(S_{t+1} = j \mid S_t = i) \; p(y_{t+1} \mid S_{t+1} = j) \; p(y_{t+2:T} \mid S_{t+1} = j, S_t = i, y_{t+1}) . \] The first two factors are \(\gamma_{ij}\) and \(f_j(y_{t+1})\). In the third, assumptions (A1) and (A2) let us drop \(S_t\) and \(y_{t+1}\): the future after \(t+1\) depends on the past only through \(S_{t+1}\). What remains is \(\beta_{t+1}(j)\).

5.2.2 Combining the two passes

Theorem 5.4 (Two-filter formula) For every \(t\) and \(i\), \[ p(S_t = i, \, Y_{1:T} = y_{1:T}) \;=\; \tilde\alpha_t(i) \, \beta_t(i), \] and hence the smoothing probability is \[ \gamma_t(i) \;:=\; \mathbb{P}(S_t = i \mid Y_{1:T} = y_{1:T}) \;=\; \frac{\tilde\alpha_t(i)\beta_t(i)}{L_T(\theta)}, \qquad L_T(\theta) = \sum_{j} \tilde\alpha_t(j)\beta_t(j) \ \text{ for any } t . \]

Proof. Factorise the joint density by splitting the observations at \(t\): \[ p(S_t = i, y_{1:T}) = p(S_t = i, y_{1:t}) \; p(y_{t+1:T} \mid S_t = i, y_{1:t}) . \] The first factor is \(\tilde\alpha_t(i)\). In the second, part 1 of Proposition 3.1 says future observations are conditionally independent of past observations given \(S_t\), so it equals \(p(y_{t+1:T} \mid S_t = i) = \beta_t(i)\).

That the likelihood can be recovered at any \(t\) is a useful implementation check: if \(\sum_j \tilde\alpha_t(j)\beta_t(j)\) drifts with \(t\), there is a bug or a scaling inconsistency.

5.2.3 Pairwise probabilities

EM needs expected transition counts, which requires the joint smoothed distribution of consecutive states.

Proposition 5.2 (Pairwise smoothing probabilities) For \(t = 1,\dots,T-1\), \[ \xi_t(i,j) \;:=\; \mathbb{P}(S_t = i, S_{t+1} = j \mid Y_{1:T} = y_{1:T}) \;=\; \frac{\tilde\alpha_t(i) \, \gamma_{ij} \, f_j(y_{t+1}) \, \beta_{t+1}(j)}{L_T(\theta)} . \] Moreover \(\sum_j \xi_t(i,j) = \gamma_t(i)\) and \(\sum_i \xi_t(i,j) = \gamma_{t+1}(j)\).

The derivation is the two-filter argument with the split made between \(t\) and \(t+1\), and the four factors have a clean reading: everything up to and including \(t\), the transition, the emission at \(t+1\), and everything after \(t+1\). The marginal consistency conditions are the second standard implementation check.

5.2.4 Scaled smoothing

Combining scaled forward variables with backward variables requires care, since the two passes must be normalised compatibly. The cleanest route is to scale the backward variables by the same constants used in the forward pass: set \(\hat\beta_T = \mathbf{1}\) and \[ \hat{\boldsymbol{\beta}}_t = \boldsymbol{\Gamma} \mathbf{Q}(y_{t+1}) \hat{\boldsymbol{\beta}}_{t+1} \big/ c_{t+1} . \] Then \(\tilde\alpha_t(i)\beta_t(i) / L_T = \hat\alpha_t(i)\hat\beta_t(i)\) exactly, so \[ \gamma_t(i) = \hat\alpha_t(i)\hat\beta_t(i), \qquad \xi_t(i,j) = \hat\alpha_t(i) \gamma_{ij} f_j(y_{t+1}) \hat\beta_{t+1}(j) / c_{t+1} , \] with no further normalisation needed. The alternative is to normalise \(\hat\beta_t\) to sum to one at each step and renormalise \(\gamma_t\) and \(\xi_t\) afterwards; that is simpler to reason about and equally stable, but loses the exact identity above.

5.3 Decoding: the Viterbi algorithm

Smoothing gives a marginal distribution at each time separately. If a single trajectory is wanted, there are two natural definitions and they do not agree.

The algorithm below was introduced by Viterbi (1967) for decoding convolutional codes, and Forney (1973) identified it as dynamic programming on a trellis — the formulation used here.

Definition 5.3 (Local and global decoding) Local (posterior) decoding takes \(\hat s_t = \arg\max_i \gamma_t(i)\) independently at each \(t\).

Global (MAP) decoding takes \[ \hat s_{1:T} \in \arg\max_{s_{1:T} \in E^T} \mathbb{P}(S_{1:T} = s_{1:T} \mid Y_{1:T} = y_{1:T}) = \arg\max_{s_{1:T}} p(s_{1:T}, y_{1:T}), \] the denominator being constant in \(s_{1:T}\).

Local decoding maximises the expected number of correctly assigned states. Global decoding maximises the probability of getting the entire sequence right. Neither dominates, and they answer different questions.

Example 5.1 (They can differ, and the local path can be impossible) Take \(K = 3\) and suppose \(\gamma_{13} = 0\). Local decoding might return \(\hat s_t = 1\) and \(\hat s_{t+1} = 3\) if state \(1\) and state \(3\) each happen to be the marginal mode at their respective times. That path has probability zero under the model. Global decoding cannot produce it, because the transition factor \(\gamma_{13} = 0\) annihilates any path containing that step.

The example is not pathological: whenever \(\boldsymbol{\Gamma}\) has structural zeros — common in profile HMMs and in duration-constrained models — local decoding can and does return infeasible paths. When the states are to be read as a trajectory, use Viterbi.

5.3.1 The recursion

Theorem 5.5 (Viterbi algorithm) Define \[ m_t(j) = \max_{s_{1:t-1} \in E^{t-1}} \; p(S_{1:t-1} = s_{1:t-1}, \, S_t = j, \, Y_{1:t} = y_{1:t}) . \] Then \[ m_1(j) = \delta_j f_j(y_1), \qquad m_t(j) = f_j(y_t) \max_{i \in E} \bigl\{ m_{t-1}(i) \gamma_{ij} \bigr\} \quad (t \ge 2), \] and with backpointers \(\psi_t(j) \in \arg\max_i \{ m_{t-1}(i)\gamma_{ij}\}\), the path defined by \[ \hat s_T \in \arg\max_j m_T(j), \qquad \hat s_t = \psi_{t+1}(\hat s_{t+1}) \ \ (t = T-1,\dots,1) \] is a global maximiser in the sense of Definition 5.3.

Proof. The recursion first. Splitting the maximisation over \(s_{1:t-1}\) into the value of \(s_{t-1}\) and the rest, \[ m_t(j) = \max_{i} \max_{s_{1:t-2}} p(s_{1:t-2}, S_{t-1}=i, S_t=j, y_{1:t}) = \max_i \Bigl[ \max_{s_{1:t-2}} p(s_{1:t-2}, S_{t-1}=i, y_{1:t-1}) \Bigr] \gamma_{ij} f_j(y_t), \] where the factorisation of Theorem 3.1 lets the last two factors be pulled out of the inner maximisation because they do not involve \(s_{1:t-2}\). The bracket is \(m_{t-1}(i)\).

For the backtracking, note that \(\max_{s_{1:T}} p(s_{1:T}, y_{1:T}) = \max_j m_T(j)\), so \(\hat s_T\) is the last state of some optimal path. Suppose inductively that \(\hat s_{t+1}\) is the state at time \(t+1\) of an optimal path. By the recursion, the maximum over paths ending in \(\hat s_{t+1}\) at time \(t+1\) is attained by extending a path that attains \(m_t(\psi_{t+1}(\hat s_{t+1}))\), so setting \(\hat s_t = \psi_{t+1}(\hat s_{t+1})\) preserves optimality. This is the Bellman principle of optimality: every suffix-optimal path has prefix-optimal prefixes.

In practice one runs the recursion in logarithms, turning products into sums: \[ \log m_t(j) = \log f_j(y_t) + \max_i \bigl\{ \log m_{t-1}(i) + \log \gamma_{ij} \bigr\} . \] No log-sum-exp is needed, because \(\max\) commutes with \(\log\). Viterbi is therefore numerically better behaved than the forward algorithm, and structural zeros are handled by \(-\infty\) entries without special cases.

5.3.2 The semiring view

Compare the two recursions: \[ \tilde\alpha_t(j) = f_j(y_t) \; \textstyle\sum_i \tilde\alpha_{t-1}(i)\gamma_{ij}, \qquad m_t(j) = f_j(y_t) \; \textstyle\max_i \, m_{t-1}(i)\gamma_{ij} . \] They are the same computation over different semirings: \((+,\times)\) for the forward algorithm and \((\max,\times)\) for Viterbi. Taking logarithms turns the latter into \((\max,+)\), and negating gives \((\min,+)\), the shortest-path algebra. So Viterbi decoding is literally a shortest-path problem on the trellis graph whose nodes are \((t,i)\) pairs and whose edge weights are \(-\log \gamma_{ij} - \log f_j(y_t)\).

The abstraction earns its keep: the same code, with the semiring swapped, computes the likelihood, the MAP path, the \(n\) best paths, expected counts, and the entropy of the posterior over paths. Murphy (2012, ch. 17) develops the general message-passing formulation.

5.3.3 Cost

Time is \(\mathcal{O}(K^2T)\), matching forward–backward. Memory is \(\mathcal{O}(KT)\), because the backpointers \(\psi_t\) must all be retained for the backward pass — worse than the \(\mathcal{O}(K)\) of a likelihood-only forward pass. For long sequences, checkpointing recomputes segments of the trellis on demand and trades a constant factor of time for \(\mathcal{O}(K\sqrt{T})\) memory.

5.4 Other quantities

Once forward and backward variables are available, most other quantities of interest are one line away.

• State occupancy. \(\mathbb{E}[\#\{t : S_t = i\} \mid y_{1:T}] = \sum_{t=1}^T \gamma_t(i)\).
• Expected transitions. \(\mathbb{E}[\#\{t : S_{t-1}=i, S_t=j\} \mid y_{1:T}] = \sum_{t=2}^{T} \xi_{t-1}(i,j)\). These two are exactly the expected sufficient statistics that EM requires in Chapter 6.
• Conditional simulation. The state sequence can be drawn from \(p(s_{1:T} \mid y_{1:T})\) by sampling \(S_T \sim \gamma_T(\cdot)\) and then, backwards, \(\mathbb{P}(S_t = i \mid S_{t+1} = j, y_{1:T}) \propto \tilde\alpha_t(i)\gamma_{ij}\). This is the forward filtering–backward sampling algorithm, developed in Chapter 9, where it is the workhorse of every MCMC scheme for HMMs.
• Pseudo-residuals. With \(u_t = F(y_t \mid y_{1:t-1})\) the predictive CDF evaluated at the observation, the \(\Phi^{-1}(u_t)\) are approximately standard normal under correct specification, giving the diagnostic plots of Zucchini et al. (2016, ch. 6).

Check your understanding. You compute \(\sum_i \tilde\alpha_t(i)\beta_t(i)\) at \(t = 10\) and get \(3.2\times10^{-41}\), and at \(t = 400\) you get \(0\). What has happened, and is the model wrong?

TipShow answer

Nothing is wrong with the model; the implementation has underflowed. Theorem 5.4 says that sum equals \(L_T(\theta)\) for every \(t\), so a value that changes with \(t\) — let alone collapses to zero — is a numerical failure, not a statistical one.

The forward variables decay geometrically, so by \(t = 400\) they are below the smallest representable double. Fix it with the scaling of Proposition 5.1 or by working in the log domain; Chapter 7 compares the two. That this invariant catches the bug immediately is why it is the first implementation check listed there.

5.5 Summary

• Filtering, smoothing and decoding answer three different questions and are computed by three different recursions, all costing \(\mathcal{O}(K^2T)\).
• The forward recursion (Theorem 5.1) is the left-to-right evaluation of the matrix product for the likelihood; its scale factors are the one-step predictive densities and their logarithms sum to \(\ell_T\) (Proposition 5.1).
• Forward and backward variables multiply to give smoothed marginals (Theorem 5.4) and, with one transition and one emission factor inserted, pairwise probabilities (Proposition 5.2).
• Viterbi maximises instead of summing (Theorem 5.5). Its output can differ from the sequence of marginal modes, and unlike that sequence it is always a feasible path (Example 5.1).
• Under a Doeblin condition, filtering forgets its initial distribution geometrically fast (Theorem 5.2), which is why \(\boldsymbol{\delta}\) matters so little in long series.

5.6 Further reading

Zucchini et al. (2016, chs. 4–5) give the same algorithms with scaling handled exactly as here, plus R implementations. Rabiner (1989) remains the most readable algorithmic account and fixed the notation most of the field still uses. Cappé et al. (2005, ch. 3) treats the recursions as operator equations on measures, which is the right level of generality for Chapter 12. Murphy (2012, ch. 17) presents them as message passing on a factor graph and makes the semiring structure explicit. Forney (1973) remains the clearest account of the trellis view of decoding. For filter stability, Del Moral and Guionnet (2001) and Cappé et al. (2005, ch. 4).

Cappé, Olivier, Éric Moulines, and Tobias Rydén. 2005. Inference in Hidden Markov Models. Springer Series in Statistics. Springer. https://doi.org/10.1007/0-387-28982-8.

Del Moral, Pierre, and Alice Guionnet. 2001. “On the Stability of Interacting Processes with Applications to Filtering and Genetic Algorithms.” Annales de l’Institut Henri Poincaré (B) Probability and Statistics 37 (2): 155–94. https://doi.org/10.1016/S0246-0203(00)01064-5.

Forney, G. David. 1973. “The Viterbi Algorithm.” Proceedings of the IEEE 61 (3): 268–78. https://doi.org/10.1109/PROC.1973.9030.

Murphy, Kevin P. 2012. Machine Learning: A Probabilistic Perspective. MIT Press.

Rabiner, Lawrence R. 1989. “A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition.” Proceedings of the IEEE 77 (2): 257–86. https://doi.org/10.1109/5.18626.

Viterbi, Andrew J. 1967. “Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm.” IEEE Transactions on Information Theory 13 (2): 260–69. https://doi.org/10.1109/TIT.1967.1054010.

Zucchini, Walter, Iain L. MacDonald, and Roland Langrock. 2016. Hidden Markov Models for Time Series: An Introduction Using R. 2nd ed. Chapman; Hall/CRC. https://doi.org/10.1201/b20790.