3 The hidden Markov model

This is the chapter the rest of the notes rests on. It states the model, writes down its two conditional independence assumptions, derives the factorisation those assumptions force, and turns that factorisation into a likelihood.

The last step is where the subject becomes computational. The likelihood is a sum over \(K^T\) state sequences, which is hopeless as written, and it is identical to a product of \(T\) small matrices, which is not. Everything in Parts III and IV follows from noticing that.

By the end of this chapter you will be able to

state the two conditional independence assumptions that define an HMM, and say what each one rules out;
derive the joint factorisation of states and observations, and read off which independences hold;
write the observed-data likelihood both as a sum over sequences and as a matrix product, and explain why the second costs \(\mathcal{O}(K^2T)\);
explain why the log-likelihood has at least \(K!\) global maxima, and what that does and does not matter for.

Notation follows Zucchini et al. (2016) throughout, so the formulas here match theirs term for term.

3.1 The two assumptions

An HMM over \(T\) time points has hidden states \(S_1,\dots,S_T\) taking values in \(E = \{1,\dots,K\}\) and observations \(Y_1,\dots,Y_T\) taking values in a measurable space \(\mathcal{Y}\). As a directed graphical model it has exactly two kinds of edge: state to next state, and state to its own observation.

flowchart LR
  S1((S₁)) --> S2((S₂)) --> S3((S₃)) --> S4((S₄))
  S1 --> Y1[Y₁]
  S2 --> Y2[Y₂]
  S3 --> Y3[Y₃]
  S4 --> Y4[Y₄]
  style Y1 fill:#dbe6f0,stroke:#14507a
  style Y2 fill:#dbe6f0,stroke:#14507a
  style Y3 fill:#dbe6f0,stroke:#14507a
  style Y4 fill:#dbe6f0,stroke:#14507a

The HMM as a directed graphical model. Shaded nodes are observed.

Definition 3.1 (The two HMM assumptions) (A1) Markov states. For every \(t\), \(\ S_{t+1} \perp\!\!\!\perp S_{1:t-1} \mid S_t\).

(A2) Conditionally independent emissions. For every \(t\), \[ Y_t \perp\!\!\!\perp \bigl( S_{1:t-1}, S_{t+1:T}, \; Y_{1:t-1}, Y_{t+1:T} \bigr) \; \bigm| \; S_t . \]

They fail in different ways and are relaxed by different extensions, so it is worth keeping them apart.

(A1) is a statement about memory in the latent process. It forces the time spent in a state to be geometrically distributed, since leaving is an independent coin flip at each step. Hidden semi-Markov models replace that with an arbitrary dwell-time distribution; higher-order chains let the next state depend on several previous ones.

(A2) says the current state screens off the current observation from everything else. Autoregressive HMMs weaken it by letting \(Y_t\) depend on \(Y_{t-1}\) as well as \(S_t\), which is what Chapter 14 needs for Markov-switching autoregressions and for modelling dinucleotide composition in genomes.

3.2 The model

Definition 3.2 (Finite-state hidden Markov model) A finite-state hidden Markov model consists of:

a finite state space \(E = \{1,\dots,K\}\);
an observation space \((\mathcal{Y}, \mathcal{B}_{\mathcal{Y}})\) with a \(\sigma\)-finite dominating measure \(\nu\);
an initial distribution \(\boldsymbol{\delta} = (\delta_1,\dots,\delta_K)^\top \in \Delta^{K-1}\);
a transition matrix \(\boldsymbol{\Gamma} = (\gamma_{ij})_{i,j=1}^K\), row-stochastic;
emission densities \(f_i = dF_i / d\nu\) for \(i \in E\).

The associated process \((S_t, Y_t)_{t \ge 1}\) is defined by: \(S_1 \sim \boldsymbol{\delta}\); \(\ \mathbb{P}(S_{t+1} = j \mid S_{1:t}, Y_{1:t}) = \gamma_{S_t j}\); and conditionally on \(S_{1:T}\) the observations \(Y_1,\dots,Y_T\) are independent with \(Y_t \sim F_{S_t}\).

We write \(\theta = (\boldsymbol{\delta}, \boldsymbol{\Gamma}, \phi_1, \dots, \phi_K)\) for the full parameter, where \(\phi_i\) parameterises \(f_i\).

Three remarks on the definition.

The dominating measure absorbs the discrete–continuous distinction. Taking \(\nu\) to be counting measure on a finite \(\mathcal{Y}\) makes \(f_i\) a probability mass function; taking Lebesgue measure on \(\mathbb{R}^d\) makes it a density. Every formula below holds in both cases unchanged, which is why these notes never split into separate discrete and continuous treatments.

The requirement \(\mathbb{P}(S_{t+1} = j \mid S_{1:t}, Y_{1:t}) = \gamma_{S_t j}\) conditions on past observations as well as past states. That is stronger than the bare Markov property of Definition 2.1 and is the precise form of (A1): the hidden chain evolves without reference to what has been seen. It is what lets the proof below drop \(y_{1:t-1}\) from the transition factor.

Nothing so far requires \(\boldsymbol{\delta}\) to be stationary. Setting \(\boldsymbol{\delta} = \boldsymbol{\pi}\) makes \((S_t, Y_t)\) strictly stationary and removes \(K-1\) free parameters. That choice is standard in the asymptotic theory of Chapter 8, where stationarity is needed to define the limiting log-likelihood, and it is often right for a long series. For a short one it forces an assumption the data may not support.

3.3 The factorisation

Theorem 3.1 (Joint factorisation) Under (A1) and (A2), for every \(T \ge 1\) the joint density of \((S_{1:T}, Y_{1:T})\) with respect to counting measure on \(E^T\) times \(\nu^{\otimes T}\) on \(\mathcal{Y}^T\) is \[ p(s_{1:T}, y_{1:T}) \;=\; \delta_{s_1} f_{s_1}(y_1) \prod_{t=2}^{T} \gamma_{s_{t-1} s_t} f_{s_t}(y_t) . \]

Proof. Apply the chain rule in the order \(S_1, Y_1, S_2, Y_2, \dots\): \[ p(s_{1:T}, y_{1:T}) = p(s_1) \, p(y_1 \mid s_1) \prod_{t=2}^T p(s_t \mid s_{1:t-1}, y_{1:t-1}) \, p(y_t \mid s_{1:t}, y_{1:t-1}). \] Assumption (A1), in the form stated in Definition 3.2, reduces \(p(s_t \mid s_{1:t-1}, y_{1:t-1})\) to \(\gamma_{s_{t-1}s_t}\). Assumption (A2) reduces \(p(y_t \mid s_{1:t}, y_{1:t-1})\) to \(f_{s_t}(y_t)\).

One line of algebra, and it is the source of every algorithm in these notes. The forward recursion sums it left to right, the backward recursion right to left, Viterbi maximises instead of summing, and EM takes expectations of its logarithm.

3.3.1 What follows, and what does not

Proposition 3.1 (Some consequences of the structure) Under (A1) and (A2):

\(Y_{1:t-1} \perp\!\!\!\perp Y_{t+1:T} \mid S_t\) — the current state separates past from future observations.
\(S_{t+1:T} \perp\!\!\!\perp Y_{1:t} \mid S_t\) — the current state is a sufficient summary of the past for predicting the future.
Given \(S_{1:T}\), the observations \(Y_1,\dots,Y_T\) are mutually independent.
Given \(Y_{1:T}\), the hidden states \(S_1,\dots,S_T\) form a Markov chain in \(t\), generally time-inhomogeneous.

Statement 1 is what makes the two-filter formula of Chapter 5 work: it licenses splitting \(p(s_t, y_{1:T})\) into a past factor and a future factor. Statement 4 is why backward sampling in Chapter 9 can draw an entire state trajectory in one sweep.

What is not true is that the observations are marginally independent. Marginalising out the states induces dependence at every lag, and that dependence is the whole signal an HMM exploits. Nor are the smoothed states independent across time, which is why the pairwise quantities \(\xi_t(i,j)\) have to be computed separately for EM.

Check your understanding. Someone tells you that in an HMM, \(Y_1\) and \(Y_2\) are independent, because emissions are conditionally independent given the states. Where is the error?

Show answer

The independence in (A2) is given the states. Statement 3 of Proposition 3.1 says \(Y_1 \perp\!\!\!\perp Y_2 \mid S_{1:T}\), not \(Y_1 \perp\!\!\!\perp Y_2\). Marginally, \[ \mathbb{P}(Y_1 \in A, Y_2 \in B) = \sum_{i,j} \delta_i \gamma_{ij} F_i(A) F_j(B), \] which factorises into \(\mathbb{P}(Y_1 \in A)\,\mathbb{P}(Y_2 \in B)\) only if \(\delta_i\gamma_{ij}\) factorises — that is, only if \(S_2\) is independent of \(S_1\). An HMM with a non-trivial transition matrix produces dependent observations by construction. Exercise 15.8 computes the autocorrelation explicitly.

3.4 The complete-data likelihood

If the states were observed, estimation would be counting. It is worth writing that case down first, because EM in Chapter 6 works by repeatedly pretending it holds.

By Theorem 3.1 the complete-data likelihood is \[ L_c(\theta; s_{1:T}, y_{1:T}) \;=\; \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}), \] with logarithm \[ \ell_c(\theta) \;=\; \log \delta_{s_1} \;+\; \sum_{t=2}^{T} \log \gamma_{s_{t-1} s_t} \;+\; \sum_{t=1}^{T} \log f_{s_t}(y_t; \phi_{s_t}) . \tag{3.1}\]

The three terms separate cleanly, and each depends on the state sequence only through counts. Writing \[ u_i = \mathbf{1}\{s_1 = i\}, \qquad n_{ij} = \sum_{t=2}^{T} \mathbf{1}\{s_{t-1} = i, \, s_t = j\}, \] 3.1 becomes \[ \ell_c(\theta) = \sum_i u_i \log \delta_i + \sum_{i,j} n_{ij} \log \gamma_{ij} + \sum_i \sum_{t: s_t = i} \log f_i(y_t;\phi_i) . \]

The initial-state indicator, the transition count matrix, and the observations grouped by state are complete-data sufficient statistics. Maximising over \(\boldsymbol{\Gamma}\) subject to the row constraints gives \(\hat\gamma_{ij} = n_{ij} / \sum_l n_{il}\), the empirical transition frequency, and the emission parameters split into \(K\) independent maximum likelihood problems. That is the whole content of the M-step; the work in Chapter 6 lies in replacing unobservable counts by their conditional expectations.

3.5 The observed-data likelihood

Definition 3.3 (Observed-data likelihood) For \(y_{1:T} \in \mathcal{Y}^T\), \[ L_T(\theta) \;=\; p_\theta(y_{1:T}) \;=\; \sum_{s_{1:T} \in E^T} \delta_{s_1} f_{s_1}(y_1) \prod_{t=2}^{T} \gamma_{s_{t-1}s_t} f_{s_t}(y_t) , \] and \(\ell_T(\theta) = \log L_T(\theta)\).

3.5.1 Why the naive sum is hopeless

The sum runs over \(K^T\) sequences. With \(K = 5\) states and \(T = 100\) — modest by any standard — that is \(5^{100} \approx 8 \times 10^{69}\) terms. No hardware helps; the growth is exponential in the length of the series, so evaluating the likelihood even once is impossible.

What redeems it is the structure of the summand. Each term is a product in which \(s_t\) appears in only three factors, so the sum distributes: summing over \(s_1\) first leaves an expression in which \(s_1\) no longer appears, and the process repeats. That is the distributive law applied systematically, and it is the whole idea behind the forward algorithm.

3.5.2 The matrix-product form

Definition 3.4 (Emission matrix) For \(y \in \mathcal{Y}\) set \(\mathbf{Q}(y) = \operatorname{diag}\bigl( f_1(y), \dots, f_K(y) \bigr)\), a \(K \times K\) diagonal matrix.

Theorem 3.2 (Matrix-product representation) For every \(T \ge 1\), \[ L_T(\theta) \;=\; \boldsymbol{\delta}^\top \mathbf{Q}(y_1) \, \boldsymbol{\Gamma} \mathbf{Q}(y_2) \, \boldsymbol{\Gamma} \mathbf{Q}(y_3) \cdots \boldsymbol{\Gamma} \mathbf{Q}(y_T) \, \mathbf{1} , \] where \(\mathbf{1} \in \mathbb{R}^K\) is the vector of ones.

Proof. Expand the product entrywise. Writing \(\mathbf{A} = \boldsymbol{\delta}^\top \mathbf{Q}(y_1) \prod_{t=2}^T \boldsymbol{\Gamma}\mathbf{Q}(y_t)\), the \(s_T\)-th entry of the row vector \(\mathbf{A}\) is obtained by summing over all intermediate indices: \[ A_{s_T} = \sum_{s_1,\dots,s_{T-1}} \delta_{s_1} f_{s_1}(y_1) \, \gamma_{s_1 s_2} f_{s_2}(y_2) \cdots \gamma_{s_{T-1}s_T} f_{s_T}(y_T), \] because each factor \(\boldsymbol{\Gamma}\mathbf{Q}(y_t)\) contributes \(\gamma_{s_{t-1}s_t} f_{s_t}(y_t)\) and matrix multiplication sums over the shared index. Right-multiplying by \(\mathbf{1}\) sums over \(s_T\) as well, giving exactly Definition 3.3.

This is Zucchini et al. (2016, eq. 2.12), and it repays study. The likelihood is a product of \(T\) matrices of size \(K \times K\), with the observation at time \(t\) entering only through the diagonal reweighting \(\mathbf{Q}(y_t)\).

Matrix multiplication is associative, so the product can be evaluated in any order — but the cost differs enormously.

Right to left, matrix by matrix. Each \(\boldsymbol{\Gamma}\mathbf{Q}(y_t)\) product is \(K \times K\) times \(K \times K\), costing \(\mathcal{O}(K^3)\) per step and \(\mathcal{O}(K^3 T)\) overall, with \(\mathcal{O}(K^2)\) memory.

Left to right, vector by matrix. Start with the row vector \(\boldsymbol{\delta}^\top \mathbf{Q}(y_1)\) and multiply by one matrix at a time. Each step is a row vector times a \(K\times K\) matrix, costing \(\mathcal{O}(K^2)\), for \(\mathcal{O}(K^2T)\) overall and only \(\mathcal{O}(K)\) memory.

The second is the forward algorithm. It is not a separate invention — it is the sensible evaluation order for this product, and the intermediate row vectors it produces are exactly the forward variables of Chapter 5.

Read the product from the left as a belief being pushed through time.

\(\boldsymbol{\delta}^\top\) is what you believe about \(S_1\) before seeing anything.
\(\mathbf{Q}(y_1)\) reweights each state by how well it explains the observation just seen.
\(\boldsymbol{\Gamma}\) pushes that belief one step forward through the dynamics.
The final \(\mathbf{1}\) sums over the last state, which you do not care about.

Alternating predict and update is exactly the structure of the Kalman filter in Chapter 12, and of Bayesian filtering in general. The finite-state case is special only in that a belief is a vector of \(K\) numbers, so the operations are matrix algebra.

Take \(K = 2\), \(\boldsymbol{\delta} = (0.5,\ 0.5)^\top\), and

\[ \boldsymbol{\Gamma} = \begin{pmatrix} 0.9 & 0.1 \\ 0.2 & 0.8 \end{pmatrix}, \qquad \mathbf{Q}(y_1) = \begin{pmatrix} 0.4 & 0 \\ 0 & 0.1 \end{pmatrix}, \qquad \mathbf{Q}(y_2) = \begin{pmatrix} 0.3 & 0 \\ 0 & 0.6 \end{pmatrix}. \]

Then \(\boldsymbol{\delta}^\top\mathbf{Q}(y_1) = (0.20,\ 0.05)\), and

\[ (0.20,\ 0.05)\,\boldsymbol{\Gamma} = (0.19,\ 0.06), \qquad (0.19,\ 0.06)\,\mathbf{Q}(y_2) = (0.057,\ 0.036), \]

so \(L_2(\theta) = 0.057 + 0.036 = 0.093\). Summing over all \(2^2 = 4\) state sequences directly gives the same number, and doing that by hand once is the fastest way to believe Theorem 3.2.

Two consequences follow immediately.

Corollary 3.1 (Smoothness) If each \(\phi \mapsto f_i(y;\phi)\) is continuous, then \(\theta \mapsto L_T(\theta)\) is continuous; if each is \(k\) times continuously differentiable, so is \(L_T\).

The likelihood is a finite sum of products of the parameters, so smoothness is inherited term by term. This licenses the Taylor expansions in Chapter 8 and the gradient-based optimisation of Chapter 7.

Corollary 3.2 (Multilinearity) \(L_T(\theta)\) is a multilinear function of \(\boldsymbol{\delta}\), of the rows of \(\boldsymbol{\Gamma}\), and of the vectors \((f_i(y_t))_i\) — linear in each when the others are held fixed.

Multilinearity is what the tensor arguments behind Theorem 4.1 exploit. It is also a warning: linearity in each block separately says nothing about joint concavity.

3.6 Symmetry and the shape of the likelihood

3.6.1 Label switching

Proposition 3.2 (Permutation invariance) Let \(\sigma\) be a permutation of \(E\) and define \(\theta^\sigma\) by \(\delta^\sigma_i = \delta_{\sigma^{-1}(i)}\), \(\gamma^\sigma_{ij} = \gamma_{\sigma^{-1}(i)\sigma^{-1}(j)}\) and \(\phi^\sigma_i = \phi_{\sigma^{-1}(i)}\). Then \(\theta\) and \(\theta^\sigma\) induce the same law for \(Y_{1:T}\), for every \(T\).

Proof. In the sum of Definition 3.3, reindex by \(s'_t = \sigma(s_t)\). Each factor maps to the corresponding factor of \(\theta^\sigma\), and the sum runs over the same set \(E^T\).

State labels carry no meaning, so nothing scientific is lost. What follows is that the best any procedure can do is identify \(\theta\) as an unordered collection of states. The practical consequences recur in Chapter 6 and, much more sharply, in Chapter 9.

3.6.2 Non-concavity

Proposition 3.3 (The log-likelihood is not concave) For \(K \ge 2\), \(\ell_T\) is in general a non-concave function of \(\theta\) with multiple local maxima.

Proof. By Proposition 3.2, \(\ell_T(\theta) = \ell_T(\theta^\sigma)\) for each of the \(K!\) permutations \(\sigma\). If \(\theta\) is a strict local maximum with distinct states, so is every \(\theta^\sigma\), and these are distinct points of \(\Theta\). A concave function on a convex set cannot have two strict local maxima, since its maximiser set is convex and any point on the segment between them would have to attain at least the common value.

The permutation symmetry alone guarantees at least \(K!\) global maxima, which is harmless — they are relabellings of each other. What is not harmless is that HMM likelihoods routinely have local maxima that are not relabellings, corresponding to genuinely worse fits from which local optimisation cannot escape. Any procedure reporting a single answer from a single starting point is reporting a local property of the surface.

Check your understanding. A colleague fits a three-state HMM twice from different random starts and gets log-likelihoods of \(-412.7\) both times, but completely different parameter values. Is something wrong?

Show answer

Almost certainly not. With \(K = 3\) there are \(3! = 6\) relabellings of any fit, all with identical log-likelihood. Check whether the two parameter sets are permutations of one another — sort the states by their emission means and compare. If they match, both runs found the same optimum and there is nothing to fix.

Had the log-likelihoods differed, say \(-412.7\) against \(-418.3\), that would be the more interesting case: two genuinely different local maxima, and a sign that more starting values are needed.

3.7 The parameter space

The parameter decomposes as \[ \Theta \;=\; \Delta^{K-1} \times \mathcal{G}_K \times \Phi_1 \times \cdots \times \Phi_K , \] where \(\mathcal{G}_K = (\Delta^{K-1})^K\) is the set of row-stochastic \(K \times K\) matrices.

Counting free parameters: \(\boldsymbol{\delta}\) contributes \(K-1\), or \(0\) under stationarity; \(\boldsymbol{\Gamma}\) contributes \(K(K-1)\); the emissions contribute \(K \dim(\Phi)\) for a common family.

Free parameters in a stationary univariate Gaussian HMM
\(K\)	Transition parameters	Emission parameters	Total
2	2	4	6
3	6	6	12
5	20	10	30
10	90	20	110

The quadratic growth in \(K\) is the practical reason state numbers stay small, and the reason model selection in Chapter 8 and Chapter 11 is not an afterthought.

The constraints are simplex constraints, and optimisers dislike them. Two standard devices:

Reparameterise. Write \(\gamma_{ij} = \exp(\tau_{ij}) / \sum_l \exp(\tau_{il})\) with \(\tau_{ii} = 0\) fixed for identification, and optimise the unconstrained \(\tau\). This is what Zucchini et al. (2016) do to hand the problem to a general-purpose optimiser, and it extends naturally to covariate-dependent transitions.
Solve the constrained problem in closed form. The M-step of EM maximises \(\sum_j n_{ij}\log\gamma_{ij}\) subject to \(\sum_j \gamma_{ij} = 1\), and a Lagrange multiplier gives \(\gamma_{ij} \propto n_{ij}\) directly. This is why EM never needs a constrained optimiser.

3.8 Summary

An HMM is defined by two conditional independence assumptions (Definition 3.1): a Markov latent chain, and observations conditionally independent given the states.
Those assumptions force the joint factorisation of Theorem 3.1, from which every algorithm in these notes follows.
With complete data, estimation is counting: the complete-data log-likelihood depends on the states only through initial, transition and emission sufficient statistics.
The observed-data likelihood is a sum over \(K^T\) sequences and equals a product of \(T\) matrices of size \(K\times K\) (Theorem 3.2). Evaluating that product left to right costs \(\mathcal{O}(K^2T)\) and is the forward algorithm.
The likelihood is smooth and multilinear in its blocks, but not concave: label switching alone forces \(K!\) global maxima (Proposition 3.2, Proposition 3.3).

3.9 Further reading

Zucchini et al. (2016, ch. 2) is the closest match to this chapter and derives the matrix representation in the same notation. Cappé et al. (2005, chs. 1–2) give the general state-space construction, where \(\boldsymbol{\Gamma}\) becomes a Markov kernel and the matrix product a composition of operators; the algebra is identical. Rabiner (1989) states the three canonical problems — evaluation, decoding, learning — that organise Chapter 5 and Chapter 6. Murphy (2012, ch. 17) and Bishop (2006, ch. 13) place HMMs in the wider graphical model framework, which is the cleanest way to see which independences are structural and which are incidental.