Instatistics, anexpectation–maximization (EM)algorithm is aniterative method to find (local)maximum likelihood ormaximum a posteriori (MAP) estimates ofparameters instatistical models, where the model depends on unobservedlatent variables.[1] The EM iteration alternates between performing an expectation (E) step, which creates a function for the expectation of thelog-likelihood evaluated using the current estimate for the parameters, and a maximization (M) step, which computes parameters maximizing the expected log-likelihood found on theE 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 ofgaussians, or to solve the multiple linear regression problem.[2]
EM clustering ofOld 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 thegeyser. Visualized usingELKI.
The EM algorithm was explained and given its name in a classic 1977 paper byArthur Dempster,Nan Laird, andDonald Rubin.[3] 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 byCedric Smith.[4] Another was proposed byH.O. Hartley in 1958, and Hartley and Hocking in 1977, from which many of the ideas in the Dempster–Laird–Rubin paper originated.[5] Another one by S.K Ng, Thriyambakam Krishnan and G.J McLachlan in 1977.[6] 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,[7][8][9] following his collaboration withPer Martin-Löf andAnders Martin-Löf.[10][11][12][13][14] 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 byC. F. Jeff Wu in 1983.[15]Wu's proof established the EM method's convergence also outside of theexponential family, as claimed by Dempster–Laird–Rubin.[15]
The EM algorithm is used to find (local)maximum likelihood parameters of astatistical model in cases where the equations cannot be solved directly. Typically these models involvelatent variables in addition to unknownparameters and known data observations. That is, eithermissing values exist among the data, or the model can be formulated more simply by assuming the existence of further unobserved data points. For example, amixture 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 thederivatives of thelikelihood 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 in this context. Additionally, it can be proven that the derivative of the likelihood is (arbitrarily close to) zero at that point, which in turn means that the point is either a local maximum or asaddle point.[15] In general, multiple maxima may occur, with no guarantee that the global maximum will be found. Some likelihoods also havesingularities in them, i.e., nonsensical maxima. For example, one of thesolutions 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.
The typical models to which EM is applied use as a latent variable indicating membership in one of a set of groups:
The observed data points may bediscrete (taking values in a finite or countably infinite set) orcontinuous (taking values in an uncountably infinite set). Associated with each data point may be a vector of observations.
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 (i.e., associated with all data points whose corresponding latent variable has that value).
However, it is possible to apply EM to other sorts of models.
The motivation 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 theViterbi algorithm forhidden 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:
First, initialize the parameters to some random values.
Compute the probability of each possible value of , given.
Then, use the just-computed values of to compute a better estimate for the parameters.
Iterate steps 2 and 3 until convergence.
The algorithm as just described monotonically approaches a local minimum of the cost function.
Although an EM iteration does increase the observed data (i.e., marginal) likelihood function, no guarantee exists that the sequence converges to amaximum likelihood estimator. Formultimodal distributions, this means that an EM algorithm may converge to alocal maximum of the observed data likelihood function, depending on starting values. A variety of heuristic ormetaheuristic approaches exist to escape a local maximum, such as random-restarthill climbing (starting with several different random initial estimates), or applyingsimulated annealing methods.
EM is especially useful when the likelihood is anexponential family, see Sundberg (2019, Ch. 8) for a comprehensive treatment:[16] the E step becomes the sum of expectations ofsufficient statistics, and the M step involves maximizing a linear function. In such a case, it is usually possible to deriveclosed-form expression updates for each step, using the Sundberg formula[17] (proved and published by Rolf Sundberg, based on unpublished results ofPer Martin-Löf andAnders Martin-Löf).[8][9][11][12][13][14]
Other methods exist to find maximum likelihood estimates, such asgradient descent,conjugate gradient, or variants of theGauss–Newton algorithm. Unlike EM, such methods typically require the evaluation of first and/or second derivatives of the likelihood function.
Expectation-Maximization works to improve rather than directly improving. Here it is shown that improvements to the former imply improvements to the latter.[18]
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 (or integrating) 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
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.[citation needed]
Instructural engineering, the Structural Identification using Expectation Maximization (STRIDE)[24] algorithm is an output-only method for identifying natural vibration properties of a structural system using sensor data (seeOperational Modal Analysis).
AKalman 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 aKalman 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 themaximum 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.[26][27][28][29]
A number of methods have been proposed to accelerate the sometimes slow convergence of the EM algorithm, such as those usingconjugate gradient and modifiedNewton's methods (Newton–Raphson).[30] 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".[31]
Expectation conditional maximization (ECM) replaces each M step with a sequence of conditional maximization (CM) steps in which each parameterθi is maximized individually, conditionally on the other parameters remaining fixed.[32] Itself can be extended into theExpectation conditional maximization either (ECME) algorithm.[33]
This idea is further extended ingeneralized expectation maximization (GEM) algorithm, in which is sought only an increase in the objective functionF for both the E step and M step as described in theAs a maximization–maximization procedure section.[19] GEM is further developed in a distributed environment and shows promising results.[34]
It is also possible to consider the EM algorithm as a subclass of theMM (Majorize/Minimize or Minorize/Maximize, depending on context) algorithm,[35] and therefore use any machinery developed in the more general case.
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[36]which contains the log-EM algorithm as its subclass. Thus, the α-EM algorithm byYasuo 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.[37]
EM is a partially non-Bayesian, maximum likelihood method. Its final result gives aprobability distribution over the latent variables (in the Bayesian style) together with a point estimate forθ (either amaximum 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,k steps per iteration are needed, wherek is the number of latent variables. Forgraphical models this is easy to do as each variable's newQ depends only on itsMarkov blanket, so localmessage passing can be used for efficient inference.
Comparison ofk-means and EM on artificial data visualized withELKI. Using the variances, the EM algorithm can describe the normal distributions exactly, while k-means splits the data inVoronoi-cells. The cluster center is indicated by the lighter, bigger symbol.An animation demonstrating the EM algorithm fitting a two component Gaussianmixture model to theOld Faithful dataset. The algorithm steps through from a random initialization to convergence.
Let be a sample of independent observations from amixture of twomultivariate normal distributions of dimension, and let be the latent variables that determine the component from which the observation originates.[20]
and
where
and
The aim is to estimate the unknown parameters representing themixing value between the Gaussians and the means and covariances of each:
Given our current estimate of the parametersθ(t), the conditional distribution of theZi is determined byBayes' theorem to be the proportional height of the normaldensity weighted byτ:
These are called the "membership probabilities", which are normally considered the output of the E step (although this is not the Q function of below).
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.
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:
This has the same form as the maximum likelihood estimate for thebinomial distribution, so
For the next estimates of:
This has the same form as a weighted maximum likelihood estimate for a normal distribution, so
The EM algorithm has been implemented in the case where an underlyinglinear 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.[38] Special cases of this model include censored or truncated observations from onenormal distribution.[38]
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 termedmoment-based approaches[39] or the so-calledspectral techniques.[40][41] 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.[citation needed]
^Jeongyeol Kwon, Constantine CaramanisProceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, PMLR 108:1727-1736, 2020.
^Sundberg, Rolf (1974). "Maximum likelihood theory for incomplete data from an exponential family".Scandinavian Journal of Statistics.1 (2):49–58.JSTOR4615553.MR0381110.
^abRolf Sundberg. 1971.Maximum likelihood theory and applications for distributions generated when observing a function of an exponential family variable. Dissertation, Institute for Mathematical Statistics, Stockholm University.
^abSundberg, Rolf (1976). "An iterative method for solution of the likelihood equations for incomplete data from exponential families".Communications in Statistics – Simulation and Computation.5 (1):55–64.doi:10.1080/03610917608812007.MR0443190.
^See the acknowledgement by Dempster, Laird and Rubin on pages 3, 5 and 11.
^abPer Martin-Löf. 1966.Statistics from the point of view of statistical mechanics. Lecture notes, Mathematical Institute, Aarhus University. ("Sundberg formula", credited to Anders Martin-Löf).
^abPer Martin-Löf. 1970.Statistiska Modeller (Statistical Models): Anteckningar från seminarier läsåret 1969–1970 (Lecture notes 1969-1970), with the assistance of Rolf Sundberg. Stockholm University.
^abMartin-Löf, P. The notion of redundancy and its use as a quantitative measure of the deviation between a statistical hypothesis and a set of observational data. With a discussion by F. Abildgård,A. P. Dempster,D. Basu,D. R. Cox,A. W. F. Edwards, D. A. Sprott,G. A. Barnard, O. Barndorff-Nielsen, J. D. Kalbfleisch andG. Rasch and a reply by the author.Proceedings of Conference on Foundational Questions in Statistical Inference (Aarhus, 1973), pp. 1–42. Memoirs, No. 1, Dept. Theoret. Statist., Inst. Math., Univ. Aarhus, Aarhus, 1974.
^abMartin-Löf, Per (1974). "The notion of redundancy and its use as a quantitative measure of the discrepancy between a statistical hypothesis and a set of observational data".Scand. J. Statist.1 (1):3–18.
^Lindstrom, Mary J; Bates, Douglas M (1988). "Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data".Journal of the American Statistical Association.83 (404): 1014.doi:10.1080/01621459.1988.10478693.
^Van Dyk, David A (2000). "Fitting Mixed-Effects Models Using Efficient EM-Type Algorithms".Journal of Computational and Graphical Statistics.9 (1):78–98.doi:10.2307/1390614.JSTOR1390614.
^Matarazzo, T. J., and Pakzad, S. N. (2016). “STRIDE for Structural Identification using Expectation Maximization: Iterative Output-Only Method for Modal Identification.” Journal of Engineering Mechanics.http://ascelibrary.org/doi/abs/10.1061/(ASCE)EM.1943-7889.0000951
^Einicke, G. A.; Malos, J. T.; Reid, D. C.; Hainsworth, D. W. (January 2009). "Riccati Equation and EM Algorithm Convergence for Inertial Navigation Alignment".IEEE Trans. Signal Process.57 (1):370–375.Bibcode:2009ITSP...57..370E.doi:10.1109/TSP.2008.2007090.S2CID1930004.
^Liu, Chuanhai; Rubin, Donald B (1994). "The ECME Algorithm: A Simple Extension of EM and ECM with Faster Monotone Convergence".Biometrika.81 (4): 633.doi:10.1093/biomet/81.4.633.JSTOR2337067.
^Matsuyama, Yasuo (2003). "The α-EM algorithm: Surrogate likelihood maximization using α-logarithmic information measures".IEEE Transactions on Information Theory.49 (3):692–706.doi:10.1109/TIT.2002.808105.
^Matsuyama, Yasuo (2011). "Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs".International Joint Conference on Neural Networks:808–816.
^Balle, Borja Quattoni, Ariadna Carreras, Xavier (2012-06-27).Local Loss Optimization in Operator Models: A New Insight into Spectral Learning.OCLC815865081.{{cite book}}: CS1 maint: multiple names: authors list (link)
Hogg, Robert; McKean, Joseph;Craig, Allen (2005).Introduction to Mathematical Statistics. Upper Saddle River, NJ: Pearson Prentice Hall. pp. 359–364.
Dellaert, Frank (February 2002).The Expectation Maximization Algorithm(PDF) (Technical Report number GIT-GVU-02-20). Georgia Tech College of Computing. gives an easier explanation of EM algorithm as to lowerbound maximization.
Gupta, M. R.; Chen, Y. (2010). "Theory and Use of the EM Algorithm".Foundations and Trends in Signal Processing.4 (3):223–296.CiteSeerX10.1.1.219.6830.doi:10.1561/2000000034. A well-written short book on EM, including detailed derivation of EM for GMMs, HMMs, and Dirichlet.
Various 1D, 2D and 3Ddemonstrations of EM together with Mixture Modeling are provided as part of the pairedSOCR activities and applets. These applets and activities show empirically the properties of the EM algorithm for parameter estimation in diverse settings.
