Expectation–maximization algorithm

In statistics, an expectation–maximization (EM) algorithm is an iterative method to find (local) maximum likelihood or maximum a posteriori (MAP) estimates of parameters in statistical models, where the model depends on unobserved latent variables. The EM iteration alternates between performing an expectation (E) step, which creates a function for the expectation of the log-likelihood evaluated using the current estimate for the parameters, and a maximization (M) step, which computes parameters maximizing the expected log-likelihood found on the step. These parameter-estimates are then used to determine the distribution of the latent variables in the next E step. It can be used, for example, to estimate a mixture of gaussians, or to solve the multiple linear regression problem.

right|frame|EM clustering of [[Old Faithful eruption data. The random initial model (which, due to the different scales of the axes, appears to be two very flat and wide ellipses) 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. Another was proposed by H.O. Hartley in 1958, and Hartley and Hocking in 1977, from which many of the ideas in the Dempster–Laird–Rubin paper originated. Another one by S.K Ng, Thriyambakam Krishnan and G.J McLachlan in 1977. Hartley's ideas can be broadened to any grouped discrete distribution. 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. The Dempster–Laird–Rubin paper established the EM method as an important tool of statistical analysis. See also Meng and van Dyk (1997).

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 also outside of the exponential family, as claimed by Dempster–Laird–Rubin. 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 (proved and published by Rolf Sundberg, based on unpublished results of Per Martin-Löf and Anders Martin-Löf).

For any <math>\mathbf{Z}</math> with non-zero probability , we can write

\log p(\mathbf{X}\mid\boldsymbol\theta) = \log p(\mathbf{X},\mathbf{Z}\mid\boldsymbol\theta) - \log p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta).

</math>

We take the expectation over possible values of the unknown data <math>\mathbf{Z}</math> under the current parameter estimate <math>\theta^{(t)}</math> by multiplying both sides by <math>p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)})</math> and summing (or integrating) over <math>\mathbf{Z}</math>. The left-hand side is the expectation of a constant, so we get:

\begin{align}

\log p(\mathbf{X}\mid\boldsymbol\theta) &

= \sum_{\mathbf{Z p{\left(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}\right)} \log p(\mathbf{X},\mathbf{Z}\mid\boldsymbol\theta)

- \sum_{\mathbf{Z p{\left(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}\right)} \log p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta) \\

& = Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}),

\end{align}

</math>

where <math>H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)})</math> is defined by the negated sum it is replacing.

This last equation holds for every value of <math>\boldsymbol\theta</math> including <math>\boldsymbol\theta = \boldsymbol\theta^{(t)}</math>,

\log p(\mathbf{X}\mid\boldsymbol\theta^{(t)})

= Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}),

</math>

and subtracting this last equation from the previous equation gives

\log p(\mathbf{X}\mid\boldsymbol\theta) - \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)})

= Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)})

+ H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}).

</math>

However, Gibbs' inequality tells us that <math>H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \ge H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)})</math>, so we can conclude that

\log p(\mathbf{X}\mid\boldsymbol\theta) - \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)})

\ge Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}).

</math>

In words, choosing <math>\boldsymbol\theta</math> to improve <math>Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)})</math> causes <math>\log p(\mathbf{X}\mid\boldsymbol\theta)</math> 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:

<math display="block"> F(q,\theta) := \operatorname{E}_q [ \log L (\theta ; x,Z) ] + H(q), </math>

where is an arbitrary probability distribution over the unobserved data and is the entropy of the distribution . This function can be written as

<math display="block"> F(q,\theta) = -D_{\mathrm{KL\big(q \parallel p_{Z\mid X}(\cdot\mid x;\theta ) \big) + \log L(\theta;x), </math>

where <math>p_{Z\mid X}(\cdot\mid x;\theta )</math> is the conditional distribution of the unobserved data given the observed data <math>x</math> and <math>D_{KL}</math> is the Kullback–Leibler divergence.

Then the steps in the EM algorithm may be viewed as:

Expectation step: Choose <math>q</math> to maximize : <math display="block"> q^{(t)} = \mathop{\arg\max}_q F{\left(q,\theta^{(t)}\right)} </math>
Maximization step: Choose <math>\theta</math> to maximize : <math display="block"> \theta^{(t+1)} = \mathop{\arg\max}_\theta F{\left(q^{(t)},\theta\right)} </math>

Applications

EM is frequently used for parameter estimation of mixed models, notably in quantitative genetics.
In psychometrics, EM is an important tool 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 (and its faster variant ordered subset expectation maximization) 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 (STRIDE) algorithm is an output-only method for identifying natural vibration properties of a structural system using sensor data (see Operational Modal Analysis).
EM is also used for data clustering. 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.
In the analysis of intertrade waiting times the EM algorithm has proved to be very useful.

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

: Operate a Kalman filter or a minimum-variance smoother designed with current parameter estimates to obtain updated state estimates.

;M-step

: Use the filtered or smoothed state estimates within maximum-likelihood calculations to obtain updated parameter estimates.

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

<math display="block">\widehat{\sigma}^2_v = \frac{1}{N} \sum_{k=1}^N {(z_k-\widehat{x}_k)}^2,</math>

where <math>\widehat{x}_k</math> are scalar output estimates calculated by a filter or a smoother from N scalar measurements <math>z_k</math>. 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

<math display="block">\widehat{\sigma}^2_w = \frac{1}{N} \sum_{k=1}^N {(\widehat{x}_{k+1}-\widehat{F}\widehat_k)}^2,</math>

where <math>\widehat{x}_k</math> and <math>\widehat{x}_{k+1}</math> are scalar state estimates calculated by a filter or a smoother. The updated model coefficient estimate is obtained via

<math display="block">\widehat{F} = \frac{\sum\limits_{k=1}^N {\left(\widehat{x}_{k+1} - \widehat{F} \widehat{x}_k\right)}^2}{\sum\limits_{k=1}^N \widehat{x}_k^2}.</math>

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 (Newton–Raphson). Also, EM can be used with constrained estimation methods.

Parameter-expanded expectation maximization (PX-EM) algorithm often provides speed up by "us[ing] a `covariance adjustment' to correct the analysis of the M step, capitalising on extra information captured in the imputed complete data".

Expectation conditional maximization (ECM) replaces each M step with a sequence of conditional maximization (CM) steps in which each parameter is maximized individually, conditionally on the other parameters remaining fixed. Itself can be extended into the Expectation conditional maximization either (ECME) algorithm.

This idea is further extended in generalized expectation maximization (GEM) algorithm, in which is sought only an increase in the objective function for both the E step and M step as described in the As a maximization–maximization procedure section.

It is also possible to consider the EM algorithm as a subclass of the MM (Majorize/Minimize or Minorize/Maximize, depending on context) 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 (in the Bayesian style) together with a point estimate for (either a maximum likelihood estimate or a posterior mode). 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 (variational Bayes), solving can iterate over each latent variable (now including ) and optimize them one at a time. Now, steps per iteration are needed, where is the number of latent variables. For graphical models this is easy to do as each variable's new 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 ===

thumb|400px|Comparison of [[K-means clustering|k-means and EM on artificial data visualized with ELKI. Using the variances, the EM algorithm can describe the normal distributions exactly, while k-means splits the data in Voronoi-cells. The cluster center is indicated by the lighter, bigger symbol.]]

thumb|240px|An animation demonstrating the EM algorithm fitting a two component Gaussian [[mixture model to the Old Faithful dataset. The algorithm steps through from a random initialization to convergence. ]]

Let <math>\mathbf{x} = (\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n)</math> be a sample of <math>n</math> independent observations from a mixture of two multivariate normal distributions of dimension <math>d</math>, and let <math>\mathbf{z} = (z_1,z_2,\ldots,z_n)</math> be the latent variables that determine the component from which the observation originates. Special cases of this model include censored or truncated observations from one normal distribution. or the so-called spectral techniques. Moment-based approaches to learning the parameters of a probabilistic model 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.

References

External links

Various 1D, 2D and 3D demonstrations of EM together with Mixture Modeling are provided as part of the paired SOCR activities and applets. These applets and activities show empirically the properties of the EM algorithm for parameter estimation in diverse settings.
Class hierarchy in C++ (GPL) including Gaussian Mixtures
The on-line textbook: Information Theory, Inference, and Learning Algorithms, by David J.C. MacKay includes simple examples of the EM algorithm such as clustering using the soft -means algorithm, and emphasizes the variational view of the EM algorithm, as described in Chapter 33.7 of version 7.2 (fourth edition).
Variational Algorithms for Approximate Bayesian Inference, by M. J. Beal includes comparisons of EM to Variational Bayesian EM and derivations of several models including Variational Bayesian HMMs (chapters).
The Expectation Maximization Algorithm: A short tutorial, A self-contained derivation of the EM Algorithm by Sean Borman.
The EM Algorithm, by Xiaojin Zhu.
EM algorithm and variants: an informal tutorial by Alexis Roche. A concise and very clear description of EM and many interesting variants.

History

As a maximization–maximization procedure

Applications

Filtering and smoothing EM algorithms

Variants

α-EM algorithm

Relation to variational Bayes methods

Geometric interpretation

Examples

See also

References

Further reading

External links