6  Maximum likelihood and EM

Chapter 5 assumed the parameters were known. This chapter drops that assumption. The likelihood is available in \(\mathcal{O}(K^2T)\) from the forward recursion, so in principle any optimiser will do; in practice the dedicated algorithm — EM, known in this context as Baum–Welch — has properties worth understanding, and its failure modes are as instructive as its guarantees.

Two things are worth stating at the outset. EM never decreases the likelihood, and that is a theorem with a two-line proof. EM does not find the global maximum, and no amount of care changes that. Both facts follow from the same decomposition.

By the end of this chapter you will be able to

• derive the EM algorithm as coordinate ascent on a lower bound, and prove it cannot decrease the likelihood;
• write down the Baum–Welch updates and explain why two of them are closed-form;
• diagnose the standard failure modes: local maxima, symmetric initialisation, collapsing states, slow convergence;
• choose between EM and direct numerical maximisation for a given problem, and obtain standard errors either way.

6.1 Maximum likelihood

Definition 6.1 (Maximum likelihood estimator) Given \(y_{1:T}\), a maximum likelihood estimator is any \[ \hat\theta_T \in \arg\max_{\theta \in \Theta} \ell_T(\theta), \qquad \ell_T(\theta) = \log p_\theta(y_{1:T}) , \] with \(\Theta\) the parameter space of Chapter 3.

The indefinite article is deliberate. By Proposition 3.2 the maximiser is never unique: each of the \(K!\) relabellings attains the same value. And by Proposition 3.3 the surface has local maxima that are not relabellings of one another.

6.1.1 Direct numerical maximisation

The forward recursion with scaling returns \(\ell_T(\theta)\) for any \(\theta\), so a general-purpose quasi-Newton optimiser can be pointed at it. The constraints are handled by reparameterisation, as in Chapter 3: transition probabilities through a multinomial logit, variances through logarithms, and so on. Zucchini et al. (2016, ch. 3) treat this as the primary method and give it careful practical treatment, which is a useful corrective to literature that presents EM as the only option.

Direct maximisation has real advantages. It converges superlinearly near the optimum where EM crawls, it produces a numerical Hessian and hence standard errors as a by-product, and it accommodates constraints and covariate effects that would complicate an M-step. Its disadvantages are that it can step outside the parameter space if the reparameterisation is careless, and that it offers no monotonicity guarantee, so a bad step can move the iterate to a worse point.

6.1.2 The gradient

Even when using EM, it is worth knowing that the score has a closed form in terms of quantities already computed.

Proposition 6.1 (Fisher’s identity) Under regularity conditions permitting differentiation under the integral sign, \[ \nabla_\theta \, \ell_T(\theta) \;=\; \mathbb{E}_\theta \bigl[ \nabla_\theta \log p_\theta(S_{1:T}, Y_{1:T}) \bigm| Y_{1:T} = y_{1:T} \bigr] . \]

Proof. Write \(p_\theta(y) = \int p_\theta(s,y) \, ds\) and differentiate: \[ \nabla_\theta \log p_\theta(y) = \frac{\nabla_\theta p_\theta(y)}{p_\theta(y)} = \int \frac{\nabla_\theta p_\theta(s,y)}{p_\theta(y)} \, ds = \int \frac{\nabla_\theta p_\theta(s,y)}{p_\theta(s,y)} \cdot \frac{p_\theta(s,y)}{p_\theta(y)} \, ds , \] which is the conditional expectation of \(\nabla_\theta \log p_\theta(s,y)\) under \(p_\theta(s \mid y)\).

The complete-data log-likelihood 3.1 is a sum of terms depending on one or two states at a time, so its conditional expectation involves only \(\gamma_t(i)\) and \(\xi_t(i,j)\) — the smoothed quantities from Chapter 5. A single forward–backward pass therefore yields both the likelihood and its gradient, which is what makes gradient-based optimisation of HMMs cheap.

6.2 The EM algorithm

6.2.1 General form

Definition 6.2 (EM algorithm) With observed data \(Y\), latent data \(S\) and current iterate \(\theta^{(k)}\):

E-step. Compute \(\ Q(\theta \mid \theta^{(k)}) = \mathbb{E}_{\theta^{(k)}}\bigl[ \log p_\theta(S, Y) \bigm| Y \bigr]\).

M-step. Set \(\ \theta^{(k+1)} \in \arg\max_{\theta \in \Theta} Q(\theta \mid \theta^{(k)})\).

The algorithm is due to Dempster et al. (1977) in general and, for HMMs specifically, predates it: Baum et al. (1970) derived the recursions and proved monotonicity for probabilistic functions of Markov chains.

6.2.2 Why it works

The mechanism is a single identity. For any distribution \(q\) over the latent states, \[ \log p_\theta(y) \;=\; \underbrace{\mathbb{E}_q\bigl[ \log p_\theta(S, y) \bigr] + H(q)}_{\textstyle \mathcal{F}(q,\theta)} \;+\; \mathrm{KL}\bigl( q \, \big\Vert \, p_\theta(\cdot \mid y) \bigr), \tag{6.1}\] where \(H(q) = -\mathbb{E}_q[\log q(S)]\). To verify it, substitute \(p_\theta(s,y) = p_\theta(s \mid y) p_\theta(y)\) into \(\mathcal{F}\) and observe that the \(\log p_\theta(y)\) term comes out of the expectation as a constant.

Since \(\mathrm{KL} \ge 0\) by Proposition 1.3, \(\mathcal{F}(q,\theta)\) is a lower bound on the log-likelihood — the evidence lower bound — and it is tight exactly when \(q = p_\theta(\cdot \mid y)\).

Theorem 6.1 (EM increases the likelihood) The EM iterates satisfy \(\ \ell(\theta^{(k+1)}) \ge \ell(\theta^{(k)})\) for every \(k\), with equality only if \(\theta^{(k+1)}\) leaves \(Q(\cdot \mid \theta^{(k)})\) at the same value.

Proof. Set \(q^{(k)} = p_{\theta^{(k)}}(\cdot \mid y)\). The E-step makes the bound tight: \(\mathcal{F}(q^{(k)}, \theta^{(k)}) = \ell(\theta^{(k)})\), because the KL term vanishes. The M-step maximises \(\mathcal{F}(q^{(k)}, \cdot)\), since \(\mathcal{F}(q^{(k)},\theta) = Q(\theta \mid \theta^{(k)}) + H(q^{(k)})\) and the entropy does not involve \(\theta\). Therefore \[ \ell(\theta^{(k+1)}) \;\ge\; \mathcal{F}(q^{(k)}, \theta^{(k+1)}) \;\ge\; \mathcal{F}(q^{(k)}, \theta^{(k)}) \;=\; \ell(\theta^{(k)}) , \] the first inequality because \(\mathcal{F}\) bounds \(\ell\) below and the second because \(\theta^{(k+1)}\) maximises it.

Read as coordinate ascent on \(\mathcal{F}(q,\theta)\) — E-step over \(q\), M-step over \(\theta\) — the algorithm loses its air of mystery, and the generalisations follow immediately. Restricting \(q\) to a tractable family gives variational EM; improving rather than maximising in either step gives generalised EM, and Theorem 6.1 survives unchanged since only the inequality \(\mathcal{F}(q^{(k)},\theta^{(k+1)}) \ge \mathcal{F}(q^{(k)},\theta^{(k)})\) was used. This is the view of Neal and Hinton (1998), and it is the bridge to the approximate methods of Chapter 7.

6.2.3 The E-step for HMMs

Taking the expectation of 3.1 term by term under \(p_{\theta^{(k)}}(\cdot \mid y_{1:T})\), \[ Q(\theta \mid \theta^{(k)}) \;=\; \underbrace{\sum_{i=1}^K \gamma_1^{(k)}(i) \log \delta_i}_{\text{initial}} \;+\; \underbrace{\sum_{t=2}^{T} \sum_{i,j} \xi_{t-1}^{(k)}(i,j) \log \gamma_{ij}}_{\text{transitions}} \;+\; \underbrace{\sum_{t=1}^{T} \sum_{i} \gamma_t^{(k)}(i) \log f_i(y_t;\phi_i)}_{\text{emissions}} , \] where \(\gamma_t^{(k)}\) and \(\xi_t^{(k)}\) are the smoothed and pairwise probabilities of Chapter 5 evaluated at \(\theta^{(k)}\). So the E-step is exactly one forward–backward pass.

The structure to notice is that the three groups of parameters appear in three separate terms. The M-step therefore decouples into three independent maximisations.

6.2.4 The M-step

Proposition 6.2 (Baum–Welch updates) Write \(\hat n_{ij} = \sum_{t=2}^{T} \xi^{(k)}_{t-1}(i,j)\) for the expected transition counts and \(w_t(i) = \gamma^{(k)}_t(i)\) for the smoothed weights. The M-step gives \[ \delta_i^{(k+1)} = \gamma_1^{(k)}(i), \qquad \gamma_{ij}^{(k+1)} = \frac{\hat n_{ij}}{\sum_{l=1}^K \hat n_{il}}, \] and \(\phi_i^{(k+1)}\) maximises the weighted log-likelihood \(\sum_{t=1}^T w_t(i) \log f_i(y_t;\phi_i)\) separately for each \(i\).

Proof. The three terms of \(Q\) involve disjoint parameter blocks, so each is maximised separately. For the initial distribution, maximise \(\sum_i \gamma_1^{(k)}(i)\log\delta_i\) subject to \(\sum_i \delta_i = 1\). The Lagrangian gives \(\gamma_1^{(k)}(i)/\delta_i = \lambda\), so \(\delta_i \propto \gamma_1^{(k)}(i)\), and since the \(\gamma_1^{(k)}(i)\) already sum to one, \(\delta_i = \gamma_1^{(k)}(i)\). The same argument applied row by row to \(\sum_j \hat n_{ij} \log \gamma_{ij}\) subject to \(\sum_j \gamma_{ij} = 1\) gives the transition update. The emission term is already in the stated form.

Both closed-form updates are instances of the same fact: maximising \(\sum_j c_j \log p_j\) over the simplex returns the normalised weights \(p_j = c_j / \sum_l c_l\). The complete-data problem of Chapter 3 was counting; the incomplete-data problem is the same counting with expected counts in place of observed ones.

For Gaussian emissions the weighted problem has the expected solution, \[ \mu_i^{(k+1)} = \frac{\sum_t w_t(i) y_t}{\sum_t w_t(i)}, \qquad (\sigma_i^2)^{(k+1)} = \frac{\sum_t w_t(i)\bigl(y_t - \mu_i^{(k+1)}\bigr)^2}{\sum_t w_t(i)} , \] and for Poisson emissions \(\lambda_i^{(k+1)} = \sum_t w_t(i) y_t / \sum_t w_t(i)\). More generally, for an exponential family in natural form the M-step solves \(\nabla A(\eta_i) = \sum_t w_t(i) T(y_t) / \sum_t w_t(i)\), matching model moments to weighted empirical ones. Zucchini et al. (2016, ch. 4) tabulate the updates for the families in common use.

WarningStationarity breaks the closed form

If \(\boldsymbol{\delta}\) is constrained to equal the stationary distribution of \(\boldsymbol{\Gamma}\), the initial and transition terms of \(Q\) are no longer separate, since \(\boldsymbol{\delta}\) becomes a function of \(\boldsymbol{\Gamma}\). The M-step then requires numerical optimisation over \(\boldsymbol{\Gamma}\). Some implementations simply ignore the constraint during the M-step and re-impose it afterwards; this breaks the monotonicity guarantee of Theorem 6.1. Direct numerical maximisation of \(\ell_T\) is the cleaner option when stationarity matters.

6.2.5 Convergence

Theorem 6.2 (Limit points are stationary points) Suppose \(\ell\) is continuous on \(\Theta\), \(Q(\cdot \mid \cdot)\) is continuous in both arguments, and the level set \(\{\theta : \ell(\theta) \ge \ell(\theta^{(0)})\}\) is compact. Then the sequence \(\ell(\theta^{(k)})\) converges monotonically to some \(\ell^*\), and every limit point of \(\{\theta^{(k)}\}\) is a stationary point of \(\ell\).

This is Wu (1983), correcting the convergence claims of Dempster et al. (1977). Three qualifications matter more than the theorem.

The likelihood values converge; the iterates need not. Convergence of the sequence \(\theta^{(k)}\) itself requires extra conditions, and on a ridge — which occurs whenever two states are nearly indistinguishable — the iterates can drift for a long time while \(\ell\) barely moves.

A stationary point is not a maximum. Saddle points are stationary, and EM initialised at a symmetric configuration will sit on one indefinitely: if two states start with identical parameters, their E-step weights are identical and the M-step keeps them identical, forever. Random perturbation of the starting values is not a nicety.

The rate is linear, with a factor governed by the ratio of missing to complete information. When states are well separated the smoothed probabilities \(\gamma_t(i)\) are near \(0\) or \(1\), little information is missing, and EM converges quickly. When states overlap, convergence is slow — and this is exactly the regime where the answer is most uncertain. Switching to a quasi-Newton method after a few dozen EM iterations is a standard remedy.

6.2.6 Initialisation

Since the algorithm finds a local maximum near where it starts, the starting point is part of the method rather than an implementation detail. What works in practice:

• Run from many dispersed starts and keep the best result. This is the only reliable protection against local maxima, and the number of starts should scale with \(K\).
• Initialise emission parameters by clustering the observations, ignoring time order, and set \(\boldsymbol{\Gamma}\) to a diagonally dominant matrix reflecting expected persistence.
• Fit \(K\) states by starting from a good \(K-1\) state fit with one state split, rather than from scratch.
• Watch for degeneracy. A state whose expected occupancy \(\sum_t \gamma_t(i)\) falls towards zero, or whose variance collapses, signals either too many states or the unboundedness discussed in Chapter 4.

Check your understanding. You run EM and the log-likelihood goes \(-812.4\), \(-806.1\), \(-806.9\), \(-804.2\). What do you conclude?

TipShow answer

There is a bug. Theorem 6.1 says the log-likelihood cannot decrease at any iteration, so the drop from \(-806.1\) to \(-806.9\) is impossible for a correctly implemented EM.

The usual culprits, in order of likelihood: a stationarity constraint being imposed after the M-step rather than inside it, which breaks the monotonicity guarantee (see the warning above); an M-step that maximises the wrong objective, often by using filtered \(\alpha_t(i)\) where smoothed \(\gamma_t(i)\) is required; or a mismatch between the scaling used in the forward and backward passes, so that \(\xi_t(i,j)\) is wrong by a \(t\)-dependent factor. Check the marginal consistency conditions of Proposition 5.2 first — they localise the last of these immediately.

6.3 Standard errors

The observed information matrix \(-\nabla^2 \ell_T(\hat\theta_T)\) inverts to an approximate covariance, justified by the asymptotic theory of Chapter 8. EM does not produce it, which is a genuine drawback. Three ways round it:

• Numerical differentiation of \(\ell_T\), cheap to implement since the forward recursion evaluates \(\ell_T\) at any \(\theta\), but sensitive to step size in high dimensions.
• Louis’s identity, which expresses the observed information as the complete-data information minus the conditional variance of the complete-data score. Both pieces are computable from forward–backward quantities.
• The parametric bootstrap, simulating series from \(\hat\theta_T\) and refitting. Expensive, but it also captures the finite-sample behaviour of the estimator, which is what Zucchini et al. (2016, ch. 3) recommend when the asymptotics look doubtful.

On the boundary — a transition probability estimated at zero, say — none of these are valid, because the asymptotic normality argument assumes an interior maximum. Bootstrap intervals are the pragmatic choice there, with the caveat that they too degrade.

6.4 Summary

• The likelihood is computable in \(\mathcal{O}(K^2T)\), so both direct numerical maximisation and EM are viable; the choice is a trade-off between monotonicity and speed near the optimum.
• EM is coordinate ascent on the evidence lower bound 6.1, which is why it cannot decrease the likelihood (Theorem 6.1).
• The E-step is one forward–backward pass; the M-step decouples into three independent problems, two of which have closed-form solutions obtained by normalising expected counts (Proposition 6.2).
• Limit points are stationary points, not maxima (Theorem 6.2). Multiple dispersed starts are the only protection against local optima, and symmetric initialisation is a trap.
• Standard errors require extra work: numerical Hessians, Louis’s identity, or the bootstrap.

6.5 Further reading

Zucchini et al. (2016, chs. 3–4) cover both direct maximisation and EM with practical guidance and R code. Cappé et al. (2005, chs. 10–11) give the rigorous treatment, including Louis’s identity and recursive variants. Dempster et al. (1977) is the original EM paper and Wu (1983) the definitive convergence analysis; Neal and Hinton (1998) is the source of the lower-bound view used here and the shortest route to understanding the generalisations. Rydén (2008) compares EM against MCMC on the same problems and is a useful reality check before committing to either.

Baum, Leonard E., Ted Petrie, George Soules, and Norman Weiss. 1970. “A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains.” The Annals of Mathematical Statistics 41 (1): 164–71. https://doi.org/10.1214/aoms/1177697196.

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.

Dempster, Arthur P., Nan M. Laird, and Donald B. Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” Journal of the Royal Statistical Society, Series B 39 (1): 1–38. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x.

Neal, Radford M., and Geoffrey E. Hinton. 1998. “A View of the EM Algorithm That Justifies Incremental, Sparse, and Other Variants.” In Learning in Graphical Models, edited by Michael I. Jordan. Springer. https://doi.org/10.1007/978-94-011-5014-9_12.

Rydén, Tobias. 2008. “EM Versus Markov Chain Monte Carlo for Estimation of Hidden Markov Models: A Computational Perspective.” Bayesian Analysis 3 (4): 659–88. https://doi.org/10.1214/08-BA326.

Wu, C. F. Jeff. 1983. “On the Convergence Properties of the EM Algorithm.” The Annals of Statistics 11 (1): 95–103. https://doi.org/10.1214/aos/1176346060.

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.