Incomplete data are a map, not merely empty cells

A table with blank entries is the obvious incomplete-data case, but it is not the paper’s governing abstraction. Dempster, Laird and Rubin begin with a complete-data space X, an observed-data space Y and a many-to-one map y=y(x). The record y tells us that the unseen complete record x belongs to the fibre X(y), the collection of complete configurations compatible with what was observed. Maximising the observed likelihood g(y|θ) may be awkward because it sums or integrates the complete density f(x|θ) over that fibre.

This definition reaches beyond accidental nonresponse. A mixture observation has a measured value but an unknown component label. A censored lifetime is known only to cross a threshold. A grouped count reports an interval rather than the underlying value. A random effect is not a missing questionnaire response, yet treating it as an unobserved part of a larger record can simplify likelihood calculation. “Incomplete” names an inferential relation between two spaces, not a defect visible in a spreadsheet.

The mapping view also separates three quantities that introductory accounts often blur. The observed values y are fixed for the likelihood calculation. The hidden complete record X is a random variable under the conditional distribution induced by y and the current parameter. The parameter θ indexes both the complete density and its marginal observed density. EM does not alternately estimate “the data” and “the parameter” as if they had equal status. It alternates an expectation over possible complete records with an optimisation over the parameter. That distinction matters when the hidden object is a whole path, a factor score or a count of unobserved truncated cases rather than one scalar blank.

It also prevents a common error about evidence. Introducing latent variables can make an optimisation elegant, but it does not create new observations. The conditional weights in the E-step are derived from the same y and the current model. When the observed record poorly distinguishes compatible complete states, the resulting uncertainty remains large. EM can propagate that uncertainty into sufficient statistics, yet no number of iterations recovers information that the observation map discarded. Computational completion is not empirical completion.

The complete-data density must induce the same observed-data model when its hidden part is integrated away. This constraint is firm; the augmentation itself may not be unique. The analyst is therefore designing a computational representation. A useful complete record makes its maximum-likelihood calculation simple and its conditional distribution given y tractable. A poor choice can preserve the observed likelihood yet produce difficult expectations or painfully small updates.

The augmentation is part of the algorithm

Given g(y|θ), many complete-data specifications can integrate to it. That freedom explains why it is misleading to call EM a single fixed recipe. An EM algorithm is defined by the observed model, the chosen augmentation and the resulting update map. Two augmentations can have the same target likelihood but different Q functions, computational costs and rates of convergence.

The paper develops the method at three levels. In a regular exponential family, the E-step needs only the conditional expectations of complete-data sufficient statistics. For a curved exponential family, those expectations still drive a constrained maximisation. At the broadest level, no exponential-family form is required: one computes the expected complete-data log likelihood Q(θ′|θ) and maximises it in θ′.

At the exponential-family level the design principle is especially transparent. Complete-data fitting depends on t(X), a finite vector of sufficient statistics. The E-step replaces t(X) by its conditional expectation tᵖ=E[t(X)|y,θᵖ]. The M-step solves the same complete-data likelihood equations, but with tᵖ in place of an observed statistic. This is why old routines for counts, regressions or covariance matrices can often be reused. The general Q formulation is broader, but the computational appeal is strongest when the hidden record compresses to a small expected statistic.

This is not ordinary single imputation. Replacing every missing value by one conditional mean can lose uncertainty and distort nonlinear sufficient statistics. If a normal observation is missing, estimating a mean may use E[X|y,θ], but estimating a variance also needs E[X²|y,θ]. Squaring the conditional mean is generally not the conditional second moment. EM expects the complete-data log likelihood, or the sufficient statistics it requires, rather than pretending one completed dataset was observed.

A valid augmentation must pass two tests. First, marginalising f(x|θ) over X(y) must reproduce g(y|θ) exactly. Otherwise the algorithm climbs the wrong likelihood. Second, both the conditional expectation and the complete-data optimisation must be feasible enough to justify the detour. An augmentation can pass the first test yet fail the second. For truncated samples, for example, one complete-data choice can leave an iterative maximum-likelihood problem inside every M-step. Extending the hidden record to include the truncated sample points may restore a simpler calculation, but it also changes the missing-information fraction and computational burden.

Expectation and maximisation are two different commitments

Let θᵖ be the current parameter. The E-step forms Q(θ′|θᵖ)=E[log f(X|θ′) | Y=y, θᵖ]. The first argument is the candidate to be optimised; the second fixes the conditional distribution used for the expectation. Confusing those roles breaks the method. During the E-step the current parameter is held fixed while hidden configurations receive posterior weights. During the M-step those weights are held fixed while the candidate parameter changes.

The heuristic “fill in, then fit” is safe only when read as expected complete-data calculation. In a finite mixture, responsibilities are fractional membership indicators. They are not claims that 0.73 of a person belongs to a component. They are the current conditional probabilities used to average the log likelihood over unknown labels. The M-step then recomputes mixture proportions and component parameters from weighted sufficient statistics.

Bayes’ rule supplies those weights. For component k and observation yᵢ, responsibility rᵢk is proportional to the current mixing weight times the component density at yᵢ, then normalised across components. The M-step reads ∑ᵢrᵢk as an expected component count. A weighted mean divides ∑ᵢrᵢkyᵢ by that count. Thus the two steps are coupled: the current fit determines which observations each component provisionally owns, and that provisional ownership determines the next fit. The observed values never move; only their conditional allocation does.

The notation Q(θ′|θᵖ) also blocks a tempting but invalid shortcut. Responsibilities computed at θᵖ belong to that conditional distribution. If an optimiser silently recomputes them while changing θ′ inside the same M-step, it is no longer maximising the Q function that the ascent decomposition uses. Some algorithms deliberately interleave or approximate the steps, but then their guarantee must be established for that algorithm. Reusing the name EM does not inherit the theorem automatically.

Dempster, Laird and Rubin also define generalised EM, or GEM. A GEM step need not find the global maximum of Q. It need only raise Q relative to the current parameter. This matters when the exact M-step is itself difficult: a partial conditional optimisation can preserve likelihood ascent. But “partial” does not mean arbitrary. If the candidate lowers Q, the ascent proof no longer applies.

Why every exact step climbs the observed likelihood

The central result is an inequality, not a metaphor. Write L(θ)=log g(y|θ). The complete conditional density of X within the observed fibre is k(x|y,θ). The expected complete log likelihood can be decomposed as Q(θ′|θ)=L(θ′)+H(θ′|θ), where H is the conditional expectation of log k. Rearrangement expresses the observed-likelihood gain as a Q gain plus an H difference.

The M-step makes the first bracket nonnegative. Jensen’s inequality, equivalently the nonnegativity of a conditional Kullback-Leibler divergence under regular conditions, makes the second bracket nonnegative. Therefore any GEM move that raises Q cannot lower L. The proof is broad because it depends on the relationship between complete and observed densities, not on one particular mixture or missing-value table.

The second term deserves attention because it is what turns a surrogate improvement into an observed-likelihood improvement. At the current parameter, the conditional density k(x|y,θ) is the reference distribution. Comparing it with k(x|y,θ′) produces an expected log ratio. Jensen’s inequality says this reference-to-candidate divergence cannot be negative. Equality requires the two conditional distributions to agree almost everywhere on the observed fibre. The M-step therefore gains observed likelihood through any Q improvement and may gain more because the hidden conditional distribution changes.

This argument is sometimes restated as a lower-bound construction. Holding the current conditional distribution fixed produces a tangent lower bound to the observed log likelihood, up to a constant. Maximising Q raises that bound, and the bound touches the likelihood at θᵖ. The language anticipates variational optimisation, but exact EM uses the actual conditional distribution under the current model, not a freely chosen approximate family. If the E-step substitutes an approximation, the bound and monotonicity properties must be re-derived.

Monotonicity is computationally valuable. A falling observed log likelihood in an exact, correctly implemented EM routine is a diagnostic: the E-step, M-step, numerical integration or likelihood evaluation is inconsistent. Approximate variants require more care. Monte Carlo noise, truncated inner optimisation or numerical underflow can produce apparent falls unless additional controls are used.

Even exact arithmetic requires the observed likelihood to be evaluated stably. Mixture densities can be tiny, so implementations should use log-sum-exp calculations rather than taking logarithms after probabilities underflow to zero. A reported decline can be a likelihood-evaluation defect rather than an EM defect. Conversely, a monotone trace cannot validate the responsibilities: two consistently wrong steps may optimise a different surrogate. Unit tests need algebraic invariants as well as the trace.

A staircase can stop below the highest summit

A nondecreasing likelihood sequence is not the same as a parameter sequence converging to the global maximum. The likelihood may be bounded while parameters move along a ridge. It may have several local maxima. A fixed point can be stationary without being globally best. The 1977 paper supplies sufficient conditions for convergence and relates a converged regular solution to likelihood stationarity, but its discussion itself notes ridges, saddle behaviour and local rather than global maxima.

C. F. Jeff Wu’s 1983 analysis is essential reading discipline here. Wu identified a flaw in the implication used in the original proof of Theorem 2 and rebuilt convergence results under clearer topological and differentiability conditions. For curved exponential families with compact parameter space, limit points are stationary under his conditions; with a unimodal likelihood and suitable differentiability, the sequence converges to the unique maximum. The later theorem is stronger precisely because the conditions are visible.

The practical rule is to run multiple starts, evaluate the observed likelihood, inspect parameters for degeneracy and distinguish a small update from a solved optimisation problem. In mixture models, swapping component labels gives an equivalent maximum; collapsing a variance towards zero can create pathological likelihood behaviour in more flexible specifications. No stopping tolerance can decide model identifiability or global optimality by itself.

There are also two kinds of convergence to distinguish. The scalar sequence L(θᵖ) may converge because it is nondecreasing and bounded, while θᵖ itself fails to approach a single point. On a likelihood ridge, separated parameter values have the same observed fit. A sequence can move within that set while the likelihood appears perfectly stable. Conditions that control the update distance or make stationary points isolated are what convert scalar convergence into pointwise parameter convergence.

A fixed point must be interpreted with the parameter space in view. Interior differentiable maxima have zero score. Boundary solutions need constrained optimality conditions instead; a zero mixture weight or variance component may be meaningful, degenerate or a sign of misspecification. The original discussion raised these cases for contingency tables and factor analysis. Declaring “gradient near zero” without checking constraints can therefore miss both a valid boundary maximum and an improper solution.

The missing-information fraction sets the pace

Near a regular fixed point, EM is generally linearly convergent. Dempster, Laird and Rubin connect its local rate to observed information and information in the unobserved part of the complete record. If little information is lost by observing y instead of x, the conditional hidden distribution is sharp and updates can be large. If the hidden part remains highly uncertain, successive E-steps change little and the largest local rate eigenvalue can sit close to one.

This explains both stability and slowness. EM uses the geometry of an easier complete-data problem and avoids aggressive curvature inversion. It often respects parameter constraints naturally: mixture weights remain probabilities and covariance updates can remain positive semidefinite. The price is that much of each complete-data step can be spent re-estimating information that remains missing.

The rate matrix can be read direction by direction. One combination of parameters may be well determined by y and update rapidly, while another is closely tied to an uncertain hidden allocation and crawls. A single norm of the parameter change hides this anisotropy. Inspecting component-wise traces, the observed score or an estimate of the local rate can reveal that one apparently settled parameter is masking another still moving slowly.

Uncertainty estimation is a separate task. The M-step’s complete-data curvature is not the observed-data curvature because hidden information has been averaged away. Treating expected records as truly observed will make standard errors too small. The missing-information principle relates complete and observed information, and later work gives practical observed-information formulae. Bootstrap or direct likelihood curvature can also be used, but the chosen method must reflect the observed model, not merely the convenient complete-data routine.

The discussants at the 1976 Royal Statistical Society meeting made the trade-off concrete. Several praised simplicity and numerical stability while warning that Newton-Raphson or scoring could be faster and supply variance estimates more directly. Sundberg reported cases with local convergence factors above 0.98. The authors’ reply accepted that EM is not an all-purpose optimal numerical method and stressed robustness against model failure and inspection of multiple maxima.

One pattern spans very different statistical problems

The paper’s examples are not decorative. They argue that a general level of description exposes a shared mechanism across techniques that had appeared separate. In missing multinomial data, expected cell indicators replace unknown classifications. In a normal linear model, expected missing values can support least squares, but multivariate covariance estimation also needs expected squares and cross-products. In grouped, censored and truncated samples, the unseen exact values or truncated cases become complete data.

For finite mixtures, component labels are hidden. The E-step computes posterior component probabilities and the M-step performs weighted complete-data fitting. The paper connects this to earlier genetic counting and Baum’s Markov-process method. In variance-component and empirical-Bayes settings, latent random quantities make the complete-data problem simpler. Iteratively reweighted least squares appears when unobserved scale variables turn a heavy-tailed marginal model into conditionally normal observations. Factor analysis treats factor scores as hidden variables and alternates familiar regressions.

Generality does not erase application assumptions. Accidental missingness requires a defensible missing-data mechanism; the paper invokes Rubin’s missing-at-random condition for the examples in Section 4.1. A mixture requires a component model and chosen number of components. Censoring and truncation require correct observation rules. Factor analysis can have rotational non-identifiability and boundary solutions. EM computes within a model; it does not validate the model.

Missing at random is not synonymous with “missing values look random”. It is a conditional statement about how the probability of a missingness pattern may depend on observed information but not on the unseen value after conditioning. The distinctness of missingness parameters also matters for ignoring that mechanism in likelihood inference. If a measurement disappears because of its unseen magnitude even after recorded predictors are held fixed, the ordinary observed-data likelihood omitting that mechanism can be wrong. EM can maximise that wrong likelihood flawlessly.

Censoring and truncation further illustrate why the observation map must be explicit. A censored case remains in the sample with only a bound or interval observed. A truncated case is absent because its value fell outside the sampling region. Their fibres and likelihood normalisations differ. Treating truncation as censoring invents observed cases; treating censoring as ordinary missingness discards the information carried by the threshold. The EM architecture can handle either, but only after the sampling story is correctly encoded.

An executable mixture laboratory

The laboratory fits a two-component normal mixture to twelve fixed one-dimensional observations. Component variance is known and shared, while the mixing proportion and two means are estimated. The hidden complete data are twelve component indicators. The E-step computes their responsibilities; the M-step turns those fractional indicators into a mixing weight and two weighted means.

Choose a clearly separated or overlapping dataset and alter both starting means. The chart shows the final responsibility of every observation and the observed log-likelihood trace. Multiple starts matter most when groups overlap or the starting means are symmetric. Swapping the two labels changes parameter names but not the fitted density or final likelihood.

Assumptions encoded. Observations are independent, there are exactly two normal components, variance is fixed at 1.44, every value is recorded, and the likelihood is evaluated without rounding. These favourable choices isolate EM’s mechanics. A real analysis must justify component number, distribution, independence and stopping criteria.

How to inspect it. First run the separated data from -1 and 1. Then set both means to 0 and observe the collapsed symmetric solution. Finally switch to overlap and start at -3 and 3. Compare final likelihoods, not just mean movement. The three self-tests check responsibility normalisation, monotone likelihood and label-swap invariance. Passing validates this finite implementation, not every EM program or mixture model.

Read the upper and lower bar segments as one probability allocation. A point with a nearly full indigo upper bar is strongly assigned to component one under the final model; a nearly full coral lower segment favours component two. Near-equal segments mark ambiguity, not an observed hybrid label. When the datasets overlap, more points remain ambiguous and the hidden labels carry more missing information. This generally makes the component parameters more coupled and the iteration less decisive.

The log-likelihood trace is the primary implementation diagnostic, but compare its scale with parameter movement. Late gains can fall below display precision while means still change. The laboratory reports the smallest numerical gain and uses a tight tolerance only for its deterministic fixed-size calculation. Production software needs tolerances scaled to likelihood magnitude, parameter scale and downstream decisions, plus an iteration cap and explicit non-convergence status.

Later work repairs, accelerates and approximates the loop

Wu’s convergence analysis clarifies what conditions turn likelihood ascent into stationary-point or unique-maximum convergence. Meng and Rubin’s 1993 ECM algorithm replaces a complicated M-step with a sequence of conditional maximisations and retains GEM-style ascent properties. When the E-step integral is intractable, Monte Carlo EM replaces it with simulation, but sampling error breaks deterministic monotonicity unless effort and acceptance controls are managed.

The same pattern connects to the surrounding canon. Rubin’s potential-outcome paper made a missing quantity central to causal definition; EM makes hidden or coarsened data central to computation, but it does not identify causal effects. Dawid’s conditional independence paper, next in the canon, gives a formal language for when conditional factorisations are legitimate. Rosenbaum and Rubin then use conditional treatment probabilities to balance observed covariates. Geman and Geman replace exact conditional expectation with stochastic local sampling in a different computational architecture.

ECM is not merely “EM with more steps”. It preserves the same E-step, then replaces one joint maximisation with conditional maximisation steps over parameter blocks while holding specified functions of the other blocks fixed. Each block step raises Q, so the whole cycle is GEM. The attraction is practical: a conditional complete-data fit may have a closed form even when the joint fit does not. The possible cost is more iterations or a different local rate.

Monte Carlo EM changes a different part of the architecture. It approximates Q when conditional expectations are analytically unavailable, often by sampling hidden records. A fixed small Monte Carlo sample injects persistent noise, so ordinary deterministic convergence is not available. Increasing simulation effort as the iteration approaches a solution is sensible because the true EM steps become smaller and must be distinguished from sampling error. Stochastic approximation and variational methods make other trade-offs; their names should not be allowed to blur their distinct guarantees.

Memory studio

Next three papers to read

1. Canon 0019, Dawid (1979). Treat conditional independence as a formal relation with its own logic.
2. Canon 0020, Rosenbaum and Rubin (1983). See how a scalar treatment probability balances observed covariates under strong assumptions.
3. Canon 0021, Geman and Geman (1984). Replace exact hidden-data averaging with local stochastic transitions for image restoration.

The achievement is a reusable computational representation

Dempster, Laird and Rubin did not invent every special iteration later recognised as EM. Their paper names the two-step method, defines incomplete data at a level broad enough to include latent and coarsened structures, proves a general ascent property, studies local rate and gathers applications that had lived in separate literatures. The unification lets an analyst search for useful complete data instead of attacking a difficult marginal likelihood directly.

The disciplined reading is conditional. First specify f(x|θ), y(x) and the induced g(y|θ). Then calculate Q with the current conditional distribution, perform a valid M or GEM step, monitor the observed likelihood and inspect the solution landscape. If the likelihood falls, debug the implementation. If it stalls, examine missing information and augmentation. If starts disagree, report the landscape rather than hiding it.

For research on intelligence, agency and consciousness, EM is valuable wherever internal states, strategies or types are latent. It can fit a probabilistic model of hidden state sequences or mixture membership. The inferred responsibilities remain model-conditioned distributions over latent variables. They do not by themselves show that the proposed states are real psychological kinds, conscious experiences or moral patients.