Expectation–maximization algorithm


In statistics, an expectation–maximization algorithm is an iterative method to find maximum likelihood or maximum a posteriori estimates of parameters in statistical models, where the model depends on unobserved latent variables. The EM iteration alternates between performing an expectation step, which creates a function for the expectation of the log-likelihood evaluated using the current estimate for the parameters, and a maximization step, which computes parameters maximizing the expected log-likelihood found on the E step. These parameter-estimates are then used to determine the distribution of the latent variables in the next E step.
eruption data. The random initial model is fit to the observed data. In the first iterations, the model changes substantially, but then converges to the two modes of the geyser. Visualized using ELKI.

History

The EM algorithm was explained and given its name in a classic 1977 paper by Arthur Dempster, Nan Laird, and Donald Rubin. They pointed out that the method had been "proposed many times in special circumstances" by earlier authors. One of the earliest is the gene-counting method for estimating allele frequencies by Cedric Smith. A very detailed treatment of the EM method for exponential families was published by Rolf Sundberg in his thesis and several papers following his collaboration with Per Martin-Löf and Anders Martin-Löf. The Dempster–Laird–Rubin paper in 1977 generalized the method and sketched a convergence analysis for a wider class of problems. Regardless of earlier inventions, the innovative Dempster–Laird–Rubin paper in the Journal of the Royal Statistical Society received an enthusiastic discussion at the Royal Statistical Society meeting with Sundberg calling the paper "brilliant". The Dempster–Laird–Rubin paper established the EM method as an important tool of statistical analysis.
The convergence analysis of the Dempster–Laird–Rubin algorithm was flawed and a correct convergence analysis was published by C. F. Jeff Wu in 1983.
Wu's proof established the EM method's convergence outside of the exponential family, as claimed by Dempster–Laird–Rubin.

Introduction

The EM algorithm is used to find maximum likelihood parameters of a statistical model in cases where the equations cannot be solved directly. Typically these models involve latent variables in addition to unknown parameters and known data observations. That is, either missing values exist among the data, or the model can be formulated more simply by assuming the existence of further unobserved data points. For example, a mixture model can be described more simply by assuming that each observed data point has a corresponding unobserved data point, or latent variable, specifying the mixture component to which each data point belongs.
Finding a maximum likelihood solution typically requires taking the derivatives of the likelihood function with respect to all the unknown values, the parameters and the latent variables, and simultaneously solving the resulting equations. In statistical models with latent variables, this is usually impossible. Instead, the result is typically a set of interlocking equations in which the solution to the parameters requires the values of the latent variables and vice versa, but substituting one set of equations into the other produces an unsolvable equation.
The EM algorithm proceeds from the observation that there is a way to solve these two sets of equations numerically. One can simply pick arbitrary values for one of the two sets of unknowns, use them to estimate the second set, then use these new values to find a better estimate of the first set, and then keep alternating between the two until the resulting values both converge to fixed points. It's not obvious that this will work, but it can be proven that in this context it does, and that the derivative of the likelihood is zero at that point, which in turn means that the point is either a maximum or a saddle point. In general, multiple maxima may occur, with no guarantee that the global maximum will be found. Some likelihoods also have singularities in them, i.e., nonsensical maxima. For example, one of the solutions that may be found by EM in a mixture model involves setting one of the components to have zero variance and the mean parameter for the same component to be equal to one of the data points.

Description

Given the statistical model which generates a set of observed data, a set of unobserved latent data or missing values, and a vector of unknown parameters, along with a likelihood function, the maximum likelihood estimate of the unknown parameters is determined by maximizing the marginal likelihood of the observed data
However, this quantity is often intractable.
The EM algorithm seeks to find the MLE of the marginal likelihood by iteratively applying these two steps:
The typical models to which EM is applied use as a latent variable indicating membership in one of a set of groups:
  1. The observed data points may be discrete or continuous. Associated with each data point may be a vector of observations.
  2. The missing values are discrete, drawn from a fixed number of values, and with one latent variable per observed unit.
  3. The parameters are continuous, and are of two kinds: Parameters that are associated with all data points, and those associated with a specific value of a latent variable.
However, it is possible to apply EM to other sorts of models.
The motive is as follows. If the value of the parameters is known, usually the value of the latent variables can be found by maximizing the log-likelihood over all possible values of, either simply by iterating over or through an algorithm such as the Baum–Welch algorithm for hidden Markov models. Conversely, if we know the value of the latent variables, we can find an estimate of the parameters fairly easily, typically by simply grouping the observed data points according to the value of the associated latent variable and averaging the values, or some function of the values, of the points in each group. This suggests an iterative algorithm, in the case where both and are unknown:
  1. First, initialize the parameters to some random values.
  2. Compute the probability of each possible value of , given.
  3. Then, use the just-computed values of to compute a better estimate for the parameters.
  4. Iterate steps 2 and 3 until convergence.
The algorithm as just described monotonically approaches a local minimum of the cost function.

Properties

Speaking of an expectation step is a bit of a misnomer. What are calculated in the first step are the fixed, data-dependent parameters of the function Q. Once the parameters of Q are known, it is fully determined and is maximized in the second step of an EM algorithm.
Although an EM iteration does increase the observed data likelihood function, no guarantee exists that the sequence converges to a maximum likelihood estimator. For multimodal distributions, this means that an EM algorithm may converge to a local maximum of the observed data likelihood function, depending on starting values. A variety of heuristic or metaheuristic approaches exist to escape a local maximum, such as random-restart hill climbing, or applying simulated annealing methods.
EM is especially useful when the likelihood is an exponential family: the E step becomes the sum of expectations of sufficient statistics, and the M step involves maximizing a linear function. In such a case, it is usually possible to derive closed-form expression updates for each step, using the Sundberg formula.
The EM method was modified to compute maximum a posteriori estimates for Bayesian inference in the original paper by Dempster, Laird, and Rubin.
Other methods exist to find maximum likelihood estimates, such as gradient descent, conjugate gradient, or variants of the Gauss–Newton algorithm. Unlike EM, such methods typically require the evaluation of first and/or second derivatives of the likelihood function.

Proof of correctness

Expectation-maximization works to improve rather than directly improving. Here is shown that improvements to the former imply improvements to the latter.
For any with non-zero probability, we can write
We take the expectation over possible values of the unknown data under the current parameter estimate by multiplying both sides by and summing over. The left-hand side is the expectation of a constant, so we get:
where is defined by the negated sum it is replacing.
This last equation holds for every value of including,
and subtracting this last equation from the previous equation gives
However, Jensen's inequality tells us that, so we can conclude that
In words, choosing to improve causes to improve at least as much.

As a maximization–maximization procedure

The EM algorithm can be viewed as two alternating maximization steps, that is, as an example of coordinate descent. Consider the function:
where q is an arbitrary probability distribution over the unobserved data z and H is the entropy of the distribution q. This function can be written as
where is the conditional distribution of the unobserved data given the observed data and is the Kullback–Leibler divergence.
Then the steps in the EM algorithm may be viewed as:

Applications

EM is frequently used for data clustering in machine learning and computer vision. In natural language processing, two prominent instances of the algorithm are the Baum–Welch algorithm for hidden Markov models, and the inside-outside algorithm for unsupervised induction of probabilistic context-free grammars.
EM is frequently used for parameter estimation of mixed models, notably in quantitative genetics.
In psychometrics, EM is almost indispensable for estimating item parameters and latent abilities of item response theory models.
With the ability to deal with missing data and observe unidentified variables, EM is becoming a useful tool to price and manage risk of a portfolio.
The EM algorithm is also widely used in medical image reconstruction, especially in positron emission tomography, single photon emission computed tomography, and x-ray computed tomography. See below for other faster variants of EM.
In structural engineering, the Structural Identification using Expectation Maximization algorithm is an output-only method for identifying natural vibration properties of a structural system using sensor data.

Filtering and smoothing EM algorithms

A Kalman filter is typically used for on-line state estimation and a minimum-variance smoother may be employed for off-line or batch state estimation. However, these minimum-variance solutions require estimates of the state-space model parameters. EM algorithms can be used for solving joint state and parameter estimation problems.
Filtering and smoothing EM algorithms arise by repeating this two-step procedure:
;E-step
;M-step
Suppose that a Kalman filter or minimum-variance smoother operates on measurements of a single-input-single-output system that possess additive white noise. An updated measurement noise variance estimate can be obtained from the maximum likelihood calculation
where are scalar output estimates calculated by a filter or a smoother from N scalar measurements. The above update can also be applied to updating a Poisson measurement noise intensity. Similarly, for a first-order auto-regressive process, an updated process noise variance estimate can be calculated by
where and are scalar state estimates calculated by a filter or a smoother. The updated model coefficient estimate is obtained via
The convergence of parameter estimates such as those above are well studied.

Variants

A number of methods have been proposed to accelerate the sometimes slow convergence of the EM algorithm, such as those using conjugate gradient and modified Newton's methods. Also, EM can be used with constrained estimation methods.
Parameter-expanded expectation maximization algorithm often provides speed up by "us a `covariance adjustment' to correct the analysis of the M step, capitalising on extra information captured in the imputed complete data".
Expectation conditional maximization replaces each M step with a sequence of conditional maximization steps in which each parameter θi is maximized individually, conditionally on the other parameters remaining fixed. Itself can be extended into the Expectation conditional maximization either algorithm.
This idea is further extended in generalized expectation maximization algorithm, in which is sought only an increase in the objective function F for both the E step and M step as described in the [|As a maximization-maximization procedure] section. GEM is further developed in a distributed environment and shows promising results.
It is also possible to consider the EM algorithm as a subclass of the MM algorithm, and therefore use any machinery developed in the more general case.

α-EM algorithm

The Q-function used in the EM algorithm is based on the log likelihood. Therefore, it is regarded as the log-EM algorithm. The use of the log likelihood can be generalized to that of the α-log likelihood ratio. Then, the α-log likelihood ratio of the observed data can be exactly expressed as equality by using the Q-function of the α-log likelihood ratio and the α-divergence. Obtaining this Q-function is a generalized E step. Its maximization is a generalized M step. This pair is called the α-EM algorithm
which contains the log-EM algorithm as its subclass. Thus, the α-EM algorithm by Yasuo Matsuyama is an exact generalization of the log-EM algorithm. No computation of gradient or Hessian matrix is needed. The α-EM shows faster convergence than the log-EM algorithm by choosing an appropriate α. The α-EM algorithm leads to a faster version of the Hidden Markov model estimation algorithm α-HMM.

Relation to variational Bayes methods

EM is a partially non-Bayesian, maximum likelihood method. Its final result gives a probability distribution over the latent variables together with a point estimate for θ. A fully Bayesian version of this may be wanted, giving a probability distribution over θ and the latent variables. The Bayesian approach to inference is simply to treat θ as another latent variable. In this paradigm, the distinction between the E and M steps disappears. If using the factorized Q approximation as described above, solving can iterate over each latent variable and optimize them one at a time. Now, k steps per iteration are needed, where k is the number of latent variables. For graphical models this is easy to do as each variable's new Q depends only on its Markov blanket, so local message passing can be used for efficient inference.

Geometric interpretation

In information geometry, the E step and the M step are interpreted as projections under dual affine connections, called the e-connection and the m-connection; the Kullback–Leibler divergence can also be understood in these terms.

Examples

Gaussian mixture

Let be a sample of independent observations from a mixture of two multivariate normal distributions of dimension, and let be the latent variables that determine the component from which the observation originates.
where
The aim is to estimate the unknown parameters representing the mixing value between the Gaussians and the means and covariances of each:
where the incomplete-data likelihood function is
and the complete-data likelihood function is
or
where is an indicator function and is the probability density function of a multivariate normal.
In the last equality, for each, one indicator is equal to zero, and one indicator is equal to one. The inner sum thus reduces to one term.

E step

Given our current estimate of the parameters θ, the conditional distribution of the Zi is determined by Bayes theorem to be the proportional height of the normal density weighted by τ:
These are called the "membership probabilities", which are normally considered the output of the E step.
This E step corresponds with setting up this function for Q:
The expectation of inside the sum is taken with respect to the probability density function, which might be different for each of the training set. Everything in the E step is known before the step is taken except, which is computed according to the equation at the beginning of the E step section.
This full conditional expectation does not need to be calculated in one step, because τ and μ/Σ appear in separate linear terms and can thus be maximized independently.

M step

Q being quadratic in form means that determining the maximizing values of θ is relatively straightforward. Also, τ, and may all be maximized independently since they all appear in separate linear terms.
To begin, consider τ, which has the constraint τ1 + τ2=1:
This has the same form as the MLE for the binomial distribution, so
For the next estimates of :
This has the same form as a weighted MLE for a normal distribution, so
and, by symmetry,

Termination

Conclude the iterative process if for below some preset threshold.

Generalization

The algorithm illustrated above can be generalized for mixtures of more than two multivariate normal distributions.

Truncated and censored regression

The EM algorithm has been implemented in the case where an underlying linear regression model exists explaining the variation of some quantity, but where the values actually observed are censored or truncated versions of those represented in the model. Special cases of this model include censored or truncated observations from one normal distribution.

Alternatives

EM typically converges to a local optimum, not necessarily the global optimum, with no bound on the convergence rate in general. It is possible that it can be arbitrarily poor in high dimensions and there can be an exponential number of local optima. Hence, a need exists for alternative methods for guaranteed learning, especially in the high-dimensional setting. Alternatives to EM exist with better guarantees for consistency, which are termed moment-based approaches or the so-called spectral techniques. Moment-based approaches to learning the parameters of a probabilistic model are of increasing interest recently since they enjoy guarantees such as global convergence under certain conditions unlike EM which is often plagued by the issue of getting stuck in local optima. Algorithms with guarantees for learning can be derived for a number of important models such as mixture models, HMMs etc. For these spectral methods, no spurious local optima occur, and the true parameters can be consistently estimated under some regularity conditions.