Optimization1977intermediate12 min read

Maximum Likelihood from Incomplete Data via the EM Algorithm

تقدير الأرجحية القصوى من بيانات ناقصة عبر خوارزمية EM

Dempster, A. P. · Laird, N. M. · Rubin, D. B. — Journal of the Royal Statistical Society, Series B

The problem

In many real-world statistical problems, some data are missing or hidden. A hospital records patient outcomes but some patients drop out before the study ends. A sensor reports values but some fall below the detection threshold. A dataset contains observations from multiple subgroups but nobody labeled which observation belongs to which group. In all these cases, directly computing the maximum estimate of the model parameters is analytically intractable — the likelihood function involves integrals or sums over the missing data that cannot be solved in closed form.

The contribution

The EM algorithm: a simple two-step iterative method that turns the intractable incomplete-data problem into a sequence of tractable complete-data problems. In the , you compute the expected value of the complete-data given the observed data and current estimates — essentially "filling in" the missing data probabilistically. In the , you maximize that expected log-likelihood to get updated parameters. Each iteration is guaranteed to increase (or maintain) the observed-data likelihood, and the algorithm converges to a stationary point. The paper unified dozens of previously ad-hoc methods under one principled framework.

The impact

The second-most cited paper in the history of statistics, with over 80,000 citations. EM became the default algorithm for fitting mixture models, handling missing data, training Hidden Markov Models, medical image reconstruction (PET/CT), and dozens of other problems. It directly inspired in Bayesian and the training procedures of Variational Autoencoders. Every time a modern ML system deals with latent variables — from topic models to to — the EM principle is at work.

You're a teacher grading essays that arrived in a shuffled pile, and several pages are missing. You can't compute exact grades. So you guess what the missing pages probably said based on the writing style you see — that's the E-step. Then you re-grade every essay using your best guess of the full text — that's the M-step. You repeat: better grades help you guess the missing pages better, and better page guesses help you grade more accurately, until grades and guesses stabilize. That's EM.

The problem: maximizing likelihood when data is missing

In statistics, given observed data $X$ and a model with parameters $\theta$ , we want the parameter values that make the data most probable — the maximum likelihood estimate (MLE). When data is complete, this often has a clean .

But what if some data is missing? The "missing data" concept is broader than it sounds. It includes:

Literally missing values — patients who dropped out of a clinical trial, sensors below a detection limit, survey non-responses.
Latent variables — cluster labels in a mixture model, hidden states in a sequence model, unobserved causes in a causal graph.
Censored observations — you know a component survived past 1,000 hours, but not exactly when it failed.

In all these cases, the observed-data likelihood $L(\theta | X)$ involves summing or integrating over all possible values of the missing data — and that sum or integral is usually analytically intractable.

Open in Lab

Toggle coins between "observed" and "hidden" to see how missing data makes the likelihood surface harder to optimize.

The demo wakes as you arrive…

The idea: guess, then improve

The insight behind EM is disarmingly simple: if having the would make the problem easy, then pretend you have it. Use your current best guess of $\theta$ to estimate what the missing data "probably" is — not by filling in a single value, but by computing the full probability of the missing data. Then, use that soft estimate to update $\theta$ , as if the data were complete.

Concretely, define:

$X$ = the observed (incomplete) data
$Z$ = the missing (unobserved) data
$(X, Z)$ = the complete data
$\theta$ = the parameters you want to estimate

The complete-data log-likelihood $\log p(X, Z | \theta)$ is usually easy to work with. The problem is that $Z$ is unknown. EM handles this by working with the expected complete-data log-likelihood — the expectation taken over $Z$ given the data you do observe and your current parameter guess.

The two steps: E and M

The algorithm alternates between two steps, starting from an initial guess $\theta^{(0)}$ :

E-step (Expectation): Compute the expected complete-data log-likelihood, treating the missing data as random variables whose distribution is determined by the current parameters $\theta^{(t)}$ and the observed data $X$ . Define the :

Q(\theta | \theta^{(t)}) = E_{Z|X,\theta^{(t)}} \left[ \log p(X, Z | \theta) \right]

The Q-function — the expected complete-data log-likelihood — Take the log-likelihood you *would* compute if you knew Z, and average it over all possible Z's, weighted by how probable each Z is given what you currently believe

Think of the E-step as asking: "Given what I believe about the parameters right now, what is the best picture I can construct of the missing data?" It is like a detective reconstructing a crime scene from partial evidence — not picking one scenario, but assigning probabilities to every possible scenario.

M-step (Maximization): Find the parameters that maximize the Q-function:

\theta^{(t+1)} = \arg\max_{\theta} \; Q(\theta | \theta^{(t)})

The M-step — maximize Q to get better parameters — Choose the

\theta

that makes the "filled-in" data most probable. Since Q decomposes nicely for exponential families, this often has a closed-form solution

The M-step is analogous to a student who, having filled in the blanks on an incomplete answer sheet (E-step), now recalculates the best-fit model from the complete answer sheet. Since working with complete data is typically easy, this step is often straightforward.

The two steps feed each other: better parameters produce a more accurate picture of the missing data, and a more accurate picture of the missing data yields better parameters. This feedback loop continues until the parameters converge.

Open in Lab

Step through E and M iterations on a 2-component Gaussian mixture. Watch the parameters converge.

The demo wakes as you arrive…

Why it works: the likelihood never decreases

The central theoretical result of the paper is the monotonicity property: every EM iteration guarantees that the observed-data log-likelihood does not decrease:

\log p(X | \theta^{(t+1)}) \geq \log p(X | \theta^{(t)})

Monotone ascent — each iteration climbs or stays level — The observed-data log-likelihood never drops across iterations, guaranteeing convergence to a stationary point (though not necessarily the global maximum)

The proof rests on a decomposition of the log-likelihood into the Q-function and an term. The E-step fixes the entropy term at its maximum (by using the true of the missing data), and the M-step increases Q. Together, the observed-data log-likelihood must go up.

However, EM comes with an important caveat: is to a local maximum or , not necessarily the global maximum. The algorithm is also linearly convergent — its rate of convergence is proportional to the fraction of missing information. When most of the data is missing, convergence can be painfully slow. Practitioners typically run EM from multiple random initializations and take the solution with the highest likelihood.

Open in Lab

Watch the log-likelihood climb with each EM iteration. Try different starting points to see how initialization affects convergence.

The demo wakes as you arrive…

Concrete example: Gaussian mixture models

The most iconic application of EM is fitting a (GMM). Imagine you have data points scattered in space that seem to form distinct clusters, but you don't know which cluster each point belongs to. You model the data as coming from $K$ Gaussian distributions, each with its own mean $\mu_k$ , covariance $\Sigma_k$ , and mixing weight $\pi_k$ .

The for each data point $x_i$ is the cluster label $z_i \in \{1, \ldots, K\}$ — you don't observe which Gaussian generated each point. This makes direct MLE intractable. But EM makes it elegant:

E-step — For each point $x_i$ and each cluster $k$ , compute the $\gamma_{ik}$ : the posterior probability that cluster $k$ generated point $x_i$ , given current parameters. This is just :

\gamma_{ik} = \frac{\pi_k \, \mathcal{N}(x_i | \mu_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(x_i | \mu_j, \Sigma_j)}

Responsibility — how much cluster k "claims" point i — Each point gets a soft assignment: a vector of K probabilities summing to 1. Unlike K-means, no point is forced to belong to exactly one cluster

M-step — Update each cluster's parameters using the responsibilities as soft weights: the new mean of cluster $k$ is the responsibility-weighted average of all data points, and the new covariance is the responsibility-weighted covariance. The mixing weights update to the average responsibility for each cluster.

These two steps repeat. The responsibilities start vague and become sharper as the clusters separate. Watch this happen live in the interactive below.

Open in Lab

Click "Run EM" to watch clusters emerge from unlabeled data. The ellipses show each Gaussian's covariance.

The demo wakes as you arrive…

Under the hood: why maximizing Q maximizes L

The mathematical engine driving EM's convergence is a decomposition of the observed-data log-likelihood. For any distribution $q(Z)$ over the missing data, we can write:

\log p(X | \theta) = Q(\theta | \theta^{(t)}) + H(\theta | \theta^{(t)})

Log-likelihood decomposition — Q plus entropy —

H

is the negative KL divergence term. The E-step fixes

H

at its maximum; the M-step pushes

Q

up. Together,

\log p(X|\theta)

must rise

This decomposition reveals something profound: EM doesn't directly optimize the observed likelihood. Instead, it constructs and maximizes a that touches the likelihood at the current parameters. The E-step tightens the bound (makes it touch), and the M-step pushes the bound (and hence the likelihood) upward. This is the same principle that later inspired variational inference and the Evidence Lower Bound () in modern Bayesian machine learning.

The same idea in code

EM for Gaussian Mixture Model — completepython

Simplified to show the idea — not the real implementation.

import numpy as np
from scipy.stats import multivariate_normal

def em_gmm(X, K, max_iter=100, tol=1e-6):
    """Fit a K-component Gaussian mixture model using EM.
    X: (N, D) data matrix — N points in D dimensions.
    Returns: means, covariances, mixing weights."""
    N, D = X.shape

    # --- Initialize: random means, identity covariances, equal weights ---
    means = X[np.random.choice(N, K, replace=False)]
    covs = [np.eye(D) for _ in range(K)]
    weights = np.ones(K) / K

    log_likelihood_old = -np.inf

    for iteration in range(max_iter):
        # === E-STEP: compute responsibilities ===
        # For each point, how probable is each cluster?
        resp = np.zeros((N, K))
        for k in range(K):
            resp[:, k] = weights[k] * multivariate_normal.pdf(X, means[k], covs[k])
        resp /= resp.sum(axis=1, keepdims=True)  # normalize rows to sum to 1

        # === M-STEP: update parameters using soft assignments ===
        Nk = resp.sum(axis=0)  # effective number of points per cluster
        for k in range(K):
            means[k] = (resp[:, k] @ X) / Nk[k]
            diff = X - means[k]
            covs[k] = (resp[:, k:k+1] * diff).T @ diff / Nk[k]
        weights = Nk / N

        # Check convergence
        log_likelihood = np.sum(np.log(
            sum(weights[k] * multivariate_normal.pdf(X, means[k], covs[k])
                for k in range(K))
        ))
        if abs(log_likelihood - log_likelihood_old) < tol:
            break
        log_likelihood_old = log_likelihood

    return means, covs, weights

# That's it. The E-step is one matrix of Bayes' theorem applications.
# The M-step is weighted averages. Everything else is bookkeeping.

Beyond mixtures: where EM appears in AI

The paper's genius was in recognizing that dozens of existing ad-hoc algorithms were all special cases of one principle. EM appears whenever a model has latent variables:

Hidden Markov Models — the Baum-Welch algorithm for speech recognition is EM applied to sequential hidden states.
Medical imaging — PET and CT reconstruction algorithms use EM to infer emission intensities from incomplete detector counts.
Topic models — Latent Dirichlet Allocation uses a variational EM to discover topics in documents.
— image , motion estimation, and stereo matching all use EM variants.
— Variational Autoencoders maximize a lower bound on the log-likelihood, directly descending from EM's Q-function framework.

The common thread: whenever you have observed data plus hidden structure, and the complete-data problem is easier than the incomplete one, EM is the natural tool.

Open in Lab

Click each application to see how EM's E and M steps specialize for different problems.

The demo wakes as you arrive…

Properties and limitations

Strengths of EM:

Simplicity. Each step typically has a closed-form solution for models.
Numerical stability. No matrix inversions of the full observed-data are needed (unlike Newton-Raphson).
. The likelihood never decreases — you always make progress or stand still.
Natural handling of constraints. Parameters like mixing weights (which must sum to 1) emerge naturally from the M-step.

Limitations:

Local optima. Multiple runs with different initializations are essential.
Linear convergence. Slower than Newton-type methods near the solution, especially when the fraction of missing information is high.
No standard errors. EM produces point estimates but not standard errors; the supplemented EM (SEM) algorithm or bootstrapping is needed for uncertainty.
Model specification. EM assumes you've chosen the right model (e.g., the right number of mixture components $K$ ). Model selection requires separate tools like BIC or .

Why it mattered

1886
Newcomb's mixture problem
Simon Newcomb attempted to fit a mixture of two normal distributions by hand — essentially an informal EM — foreshadowing the need for a systematic algorithm.
1958
Hartley's iterative MLE
H. O. Hartley proposed an iterative maximum likelihood method for mixtures, one of the earliest recognizable precursors of EM.
1970
Baum-Welch algorithm
Baum et al. developed an iterative re-estimation procedure for Hidden Markov Models — later recognized as a special case of EM.
1977
Dempster, Laird & Rubin
The EM paper unified all prior ad-hoc methods, proved monotone convergence, and named the algorithm. It became the second most cited paper in statistics.
1983
Wu's convergence theory
C. F. Jeff Wu established rigorous convergence conditions for EM, strengthening the original paper's guarantees and laying the mathematical foundation.
1993
Variational EM & mean-field
Neal & Hinton recast EM as coordinate ascent on a free-energy functional, opening the door to variational Bayes and approximate inference.
2014
Variational Autoencoders
Kingma & Welling used the ELBO — EM's lower-bound principle — with neural network encoders and decoders, bridging classical statistics and deep generative models.

Dempster, Laird, and Rubin didn't invent any of the individual applications — many had been discovered independently. Their contribution was recognizing that all these algorithms shared the same underlying principle, proving it works, and giving it a name. That unification made it possible to apply the same framework to any new problem with missing data, immediately.

CitationDempster, Laird, Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society, Series B, 1977.

Terms in this paper

E-stepخطوة التوقع
M-stepخطوة التعظيم
Q-functionالدالة Q
Incomplete Dataالبيانات الناقصة
Complete Dataالبيانات الكاملة
Responsibilityالمسؤولية الاحتمالية
Monotone Convergenceالتقارب الرتيب
Gaussian Mixture Modelنموذج المزيج الغاوسي