9  Bayesian inference

Maximum likelihood returns a point and, with work, a standard error. A Bayesian analysis returns a distribution, and for HMMs it also solves several problems that are awkward otherwise: unbounded likelihoods, states that empty out, boundary parameters, and — the subject of Chapter 11 — an unknown number of states.

The computational story is unusually clean. The latent states are the natural augmentation variable, and conditional on them the model becomes a collection of independent standard problems. The one algorithm that makes this work, forward filtering with backward sampling, is a small modification of the recursions already derived in Chapter 5.

This chapter sets up priors, derives the sampler, and confronts label switching, which is more consequential in the Bayesian setting than the classical one. Frühwirth-Schnatter (2006, chs. 10–11) is the standard reference throughout.

By the end of this chapter you will be able to

  • specify conjugate priors for an HMM and explain what the Dirichlet hyperparameters encode;
  • derive forward filtering with backward sampling, and say why the states must be drawn as a block;
  • write out a full Gibbs sweep and compare each step with the corresponding EM update;
  • recognise label switching in MCMC output and choose an appropriate response.

9.1 Priors

Write \(\theta = (\boldsymbol{\delta}, \boldsymbol{\Gamma}, \phi_{1:K})\) as before. A prior must be placed on each block.

9.1.1 Transition matrix

Rows of \(\boldsymbol{\Gamma}\) live on the simplex, and the Dirichlet distribution is conjugate to the multinomial counts that arise from an observed state sequence.

Definition 9.1 (Independent Dirichlet rows) Take the rows independently, \[ \boldsymbol{\Gamma}_{j,\cdot} \sim \mathrm{Dir}(e_{j1}, \dots, e_{jK}), \qquad j = 1,\dots,K , \] with hyperparameters \(e_{jk} > 0\). The symmetric choice \(e_{jk} = e_0\) for all \(j,k\) is the default.

The hyperparameters act as prior transition counts, and their size controls how strongly the posterior is pulled towards uniformity. Two choices are worth distinguishing. Setting \(e_0 = 1\) gives a uniform prior on each row, which sounds neutral but in fact favours transition matrices with little persistence — a problem when the states are meant to represent regimes that last. Putting extra mass on the diagonal, \(e_{jj} = e_0 + \varkappa\) with \(\varkappa > 0\), encodes persistence directly; this is the “sticky” prior that Fox et al. (2011) introduced for the nonparametric case and it applies equally to finite models.

Small values \(e_0 < 1\) favour sparse rows, pushing mass towards transition matrices with near-zero entries. That behaviour is exploited deliberately in Chapter 11.

The initial distribution takes a Dirichlet prior in the same way, \(\boldsymbol{\delta} \sim \mathrm{Dir}(d_1,\dots,d_K)\), unless stationarity is imposed, in which case \(\boldsymbol{\delta}\) is determined by \(\boldsymbol{\Gamma}\) and carries no prior of its own.

9.1.2 Emission parameters

The natural choice is the conjugate prior for the emission family, applied independently across states: \[ \phi_i \stackrel{\text{iid}}{\sim} G_0, \qquad i = 1,\dots,K . \] For Gaussian emissions with unknown mean and variance, \(G_0\) is normal–inverse-gamma; for Poisson emissions, gamma; for categorical emissions, Dirichlet. Conjugacy is not required, but it makes the corresponding sampler step a direct draw rather than a Metropolis update.

TipThe prior fixes the unbounded likelihood

Chapter 4 noted that a Gaussian HMM has an unbounded likelihood: sending \(\sigma_i \to 0\) with \(\mu_i\) at an observation drives \(\ell_T \to \infty\). A proper inverse-gamma prior on \(\sigma_i^2\) has a density vanishing at zero, so the posterior remains proper and the degeneracy disappears. This is a real advantage of the Bayesian formulation, not a technicality — it is why Bayesian HMM fits do not need variance floors.

Independence across states in \(G_0\) says that knowing state \(1\)’s parameters tells you nothing about state \(2\)’s. When there are several related series, that assumption is usually wrong and wasteful; Chapter 10 replaces it with a hierarchical prior that lets the states borrow strength.

9.2 Posterior simulation

The posterior \(p(\theta \mid y_{1:T}) \propto p(\theta) L_T(\theta)\) is not tractable, because \(L_T\) is a sum over \(K^T\) sequences. Data augmentation makes it tractable by treating the states as unknowns to be sampled rather than summed out.

9.2.1 The augmented target

Target the joint posterior of parameters and states, \[ p(\theta, s_{1:T} \mid y_{1:T}) \;\propto\; p(\theta) \; \delta_{s_1} f_{s_1}(y_1;\phi_{s_1}) \prod_{t=2}^{T} \gamma_{s_{t-1}s_t} f_{s_t}(y_t;\phi_{s_t}) , \] and read off the two conditionals. Given \(\theta\), the states have the distribution that FFBS below can sample from. Given the states, the complete-data structure of Chapter 3 applies, and the parameter blocks separate exactly as they did in the M-step of Proposition 6.2 — with conjugate priors turning “maximise” into “sample”.

9.2.2 Forward filtering, backward sampling

The states must be drawn jointly, not one at a time. Sampling \(S_t\) from its full conditional given all other states, a single-site Gibbs move, is correct but mixes disastrously: consecutive states are strongly dependent, so the chain moves through state space one time point at a time and takes \(\mathcal{O}(T^2)\) sweeps to traverse it. Drawing the whole trajectory in one block removes the problem.

Theorem 9.1 (Forward filtering, backward sampling) Let \(\alpha_t(i) = \mathbb{P}_\theta(S_t = i \mid y_{1:t})\) be the filtering probabilities from Chapter 5. Then \[ p_\theta(s_{1:T} \mid y_{1:T}) \;=\; p_\theta(s_T \mid y_{1:T}) \prod_{t=T-1}^{1} p_\theta(s_t \mid s_{t+1}, y_{1:t}), \] where \[ \mathbb{P}_\theta(S_t = i \mid S_{t+1} = j, \, y_{1:t}) \;=\; \frac{\alpha_t(i) \, \gamma_{ij}}{\sum_{l=1}^{K} \alpha_t(l)\gamma_{lj}} . \]

Proof. For the factorisation, apply the chain rule backwards and use part 2 of Proposition 3.1: conditional on \(S_{t+1}\), the state \(S_t\) is independent of \(S_{t+2:T}\) and of \(y_{t+1:T}\), so the conditioning set collapses from \((s_{t+1:T}, y_{1:T})\) to \((s_{t+1}, y_{1:t})\).

For the conditional itself, Bayes’ theorem gives \[ \mathbb{P}(S_t = i \mid S_{t+1}=j, y_{1:t}) \propto \mathbb{P}(S_t = i \mid y_{1:t}) \, \mathbb{P}(S_{t+1} = j \mid S_t = i, y_{1:t}) = \alpha_t(i)\gamma_{ij}, \] the second factor reducing to \(\gamma_{ij}\) by the Markov property. Normalising over \(i\) gives the stated expression.

So one forward pass storing the filtered vectors \(\alpha_1,\dots,\alpha_T\), followed by a backward sweep drawing \(S_T \sim \alpha_T\) and then each \(S_t\) from a \(K\)-point distribution, produces an exact draw from \(p_\theta(s_{1:T} \mid y_{1:T})\). Cost is \(\mathcal{O}(K^2T)\) time and \(\mathcal{O}(KT)\) memory, the same as forward–backward. The algorithm is due to Frühwirth-Schnatter (1994) and Carter and Kohn (1994) for state-space models and Chib (1996) for the discrete case; no backward variables are needed, only backward sampling.

9.2.3 The Gibbs sampler

Definition 9.2 (Data augmentation sampler for a finite HMM) Given the current \(\theta^{(m)}\), one sweep produces \(\theta^{(m+1)}\) as follows.

  1. States. Draw \(s_{1:T}^{(m+1)} \sim p_{\theta^{(m)}}(\cdot \mid y_{1:T})\) by FFBS (Theorem 9.1).
  2. Transitions. Form the counts \(N_{jk} = \#\{t \ge 2 : s_{t-1} = j, s_t = k\}\) and draw the rows independently, \[ \boldsymbol{\Gamma}^{(m+1)}_{j,\cdot} \sim \mathrm{Dir}\bigl( e_{j1} + N_{j1}, \dots, e_{jK} + N_{jK} \bigr). \]
  3. Initial distribution. Draw \(\boldsymbol{\delta}^{(m+1)} \sim \mathrm{Dir}(d_1 + \mathbf{1}\{s_1 = 1\}, \dots, d_K + \mathbf{1}\{s_1 = K\})\).
  4. Emissions. For each \(i\), draw \(\phi_i^{(m+1)}\) from its conjugate posterior given the observations assigned to state \(i\), namely \(\{ y_t : s_t = i \}\).

Each step is a standard conjugate update. Step 2 is Dirichlet–multinomial conjugacy applied row by row, which works because rows are a priori independent and the complete-data likelihood factorises across them. Step 4 is \(K\) separate conjugate updates, each involving only the observations currently assigned to that state — for Gaussian emissions, a normal–inverse-gamma update using the within-state sample mean and sum of squares.

Comparing with Proposition 6.2 is instructive. The M-step normalised expected counts; the Gibbs step draws from a Dirichlet centred near the normalised realised counts. EM propagates a soft assignment \(\gamma_t(i)\); the sampler propagates a hard assignment drawn from the same distribution. The two algorithms have the same skeleton, and Rydén (2008) compares their behaviour on identical problems in detail.

9.2.4 Mixing and diagnostics

Blocking the states is what makes the sampler work at all, but it does not make it fast in every regime. Mixing degrades when states overlap heavily, since the assignments then change little from sweep to sweep, and when the transition matrix is very persistent, since a long run in one state resists being reassigned. Standard remedies are to monitor the number of occupied states and the log-likelihood rather than individual parameters, to run several chains from dispersed starts, and — where the emission model permits — to integrate out \(\phi\) analytically and sample the states from the collapsed target, which usually improves mixing at the cost of losing the conditional independence that makes step 4 cheap.

9.3 Label switching

Definition 9.3 (The problem) By Proposition 3.2 the likelihood is invariant under relabelling. If the prior is also invariant — which it is under symmetric Dirichlet rows and exchangeable \(G_0\) — then the posterior has \(K!\) identical modes, and any state-specific posterior summary is meaningless: \(\mathbb{E}[\phi_1 \mid y]\), \(\mathbb{E}[\phi_2 \mid y], \dots\) are all equal to the same average over all states.

This is worse than the classical version of the problem. An optimiser converges to one mode and stays there, so the symmetry is invisible. A sampler that mixes properly must visit all \(K!\) modes, and then the marginal posterior of \(\phi_1\) is a \(K\)-fold symmetric mixture. Paradoxically, a sampler that appears to give clean, unimodal state-specific posteriors is often one that has failed to mix.

Three responses, in increasing order of respectability.

Identifiability constraints. Impose an ordering, such as \(\mu_1 < \mu_2 < \dots < \mu_K\), and permute each draw to satisfy it. Simple and often adequate, but the choice of ordering functional is arbitrary and a poor choice can cut through a mode rather than separating modes, producing artefacts. Jasra et al. (2005) give examples.

Relabelling algorithms. Post-process the output by choosing, for each draw, the permutation that minimises a loss relative to a reference — a Kullback–Leibler criterion in Stephens (2000), or a clustering of the draws in the space of allocation vectors. More work, but it adapts to the geometry of the posterior instead of assuming it.

Report only invariant quantities. Many questions do not require labelled states. The number of occupied states, the posterior predictive distribution, the co-clustering probabilities \(\mathbb{P}(S_t = S_{t'} \mid y_{1:T})\), and the log marginal likelihood are all permutation-invariant and need no relabelling at all. When the scientific question can be phrased this way, it should be.

Deliberately randomising the labels — applying a uniformly random permutation after each sweep, the random permutation sampler of Frühwirth-Schnatter (2006) — is a useful diagnostic. It forces the sampler to explore the full symmetric posterior, so the resulting output can be checked for the balance that a well-mixing chain should exhibit, and it provides a clean input to a relabelling step.

Check your understanding. Your sampler for a three-state model gives clean, well-separated posteriors for \(\mu_1, \mu_2, \mu_3\), with narrow credible intervals and no multimodality. Should you be pleased?

Probably not. The posterior is invariant under the \(3! = 6\) relabellings (Definition 9.3), so a chain that is exploring it properly should show each marginal as the same six-fold symmetric mixture. Three distinct unimodal posteriors mean the chain has stayed in one labelling for the whole run — it has not mixed over the symmetric modes.

That may still be fine for practical purposes, and it is what an identifiability constraint deliberately produces. But it is not evidence of good mixing, and it should not be reported as though the sampler had converged. Run the random permutation sampler as a diagnostic: if the chain then fails to explore all six modes, the mixing problem is real and affects more than the labels.

9.4 Marginal likelihood

Bayesian model comparison across different \(K\) needs the marginal likelihood \(p(y_{1:T} \mid K) = \int p_\theta(y_{1:T}) p(\theta \mid K) \, d\theta\). It is not available from the Gibbs output directly, and naive estimators such as the harmonic mean are notoriously unstable.

The standard approach is Chib (1996), which uses the basic marginal identity \[ \log p(y_{1:T}) = \log p_{\theta^*}(y_{1:T}) + \log p(\theta^*) - \log p(\theta^* \mid y_{1:T}) \] at a fixed high-density point \(\theta^*\). The first term comes from the forward recursion, the second from the prior, and the third from a sequence of reduced Gibbs runs that estimate the posterior ordinate blockwise. The method works well for HMMs precisely because the likelihood is available in closed form; its weakness is that the posterior ordinate estimate must account for label switching, or it will be wrong by a factor of \(K!\).

Where the goal is to choose \(K\) rather than to report evidence for a fixed \(K\), putting a prior on \(K\) and sampling it is usually better than computing marginal likelihoods for a grid of values. That is the subject of Chapter 11.

9.5 Summary

  • Conjugate priors — Dirichlet on the rows of \(\boldsymbol{\Gamma}\), a conjugate family for the emissions — make every parameter update a direct draw, and a proper prior on the variances removes the unbounded-likelihood pathology of Chapter 4.
  • Data augmentation targets the joint posterior of parameters and states. Conditional on the states the model decomposes exactly as the complete-data problem of Chapter 3.
  • States must be drawn as a block. Forward filtering with backward sampling (Theorem 9.1) does this exactly in \(\mathcal{O}(K^2T)\), using only the filtered probabilities.
  • Label switching is sharper here than in the classical setting: a well-mixing chain has \(K!\) symmetric modes. Report invariant quantities where possible, and relabel deliberately where not.

9.6 Further reading

Frühwirth-Schnatter (2006, chs. 10–11) is the reference treatment of Bayesian HMMs, including priors, samplers and label switching, and Frühwirth-Schnatter et al. (2019) collects the current state of the mixture-model literature it belongs to. Robert and Casella (2004) is the general reference for the Monte Carlo methods used here. Scott (2002) is a readable applied account of recursive Bayesian computation for HMMs. Chib (1996) introduces both the FFBS scheme in this setting and the marginal likelihood estimator. On label switching, Stephens (2000) and Jasra et al. (2005). Rydén (2008) compares EM and MCMC directly and is worth reading before choosing between them. Robert et al. (2000) handle unknown \(K\) by reversible jump, the approach that Chapter 11 is designed to avoid.