View source for Expectation–maximization algorithm

{{Short description|Iterative method for finding maximum likelihood estimates in statistical models}}
In [[Statistics|statistics]], an '''expectation–maximization''' ('''EM''') '''algorithm''' is an [[Iterative method|iterative method]] to find (local) maximum likelihood or maximum a posteriori (MAP) estimates of [[Parameter|parameter]]s in [[Statistical model|statistical model]]s, where the model depends on unobserved [[Social:Latent variable|latent variable]]s.<ref>{{cite journal |last1=Meng |first1=X.-L. |last2=van Dyk |first2=D. |title=The EM algorithm – an old folk-song sung to a fast new tune |journal=J. Royal Statist. Soc. B |date=1997 |volume=59 |issue=3 |pages=511–567|doi=10.1111/1467-9868.00082 |s2cid=17461647 |doi-access=free }}</ref> The EM iteration alternates between performing an expectation (E) step, which creates a function for the expectation of the [[Likelihood function#Log-likelihood|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 ''E'' 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 [[Normal distribution|gaussians]], or to solve the multiple linear regression problem.<ref>Jeongyeol Kwon, Constantine Caramanis

''Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics'', PMLR 108:1727-1736, 2020.</ref>

[[File:EM Clustering of Old Faithful data.gif|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 elipses) 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 [[Biography:Arthur P. Dempster|Arthur Dempster]], [[Biography:Nan Laird|Nan Laird]], and [[Biography:Donald Rubin|Donald Rubin]].<ref name="Dempster1977">
{{cite journal
 |title=Maximum Likelihood from Incomplete Data via the EM Algorithm
 |journal=Journal of the Royal Statistical Society, Series B
 |year=1977 |volume=39 |issue=1 |pages=1–38
 |jstor=2984875 |mr= 0501537
}}
</ref> 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 [[Biography:Cedric Smith (statistician)|Cedric Smith]].<ref>{{cite journal |last1=Ceppelini |first1=R.M. |title=The estimation of gene frequencies in a random-mating population |journal=Ann. Hum. Genet. |volume=20 |issue=2 |pages=97–115 |doi=10.1111/j.1469-1809.1955.tb01360.x|pmid=13268982 |year=1955 |s2cid=38625779 }}</ref>  Another was proposed by [[Biography:Herman Otto Hartley|H.O. Hartley]] in 1958, and Hartley and Hocking in 1977, from which many of the ideas in the Dempster–Laird–Rubin paper originated.<ref>{{Cite journal |last=Hartley |first=Herman Otto |date=1958 |title=Maximum Likelihood estimation from incomplete data |journal=Biometrics |volume=14 |issue=2 |pages=174–194|doi=10.2307/2527783 |jstor=2527783 }}</ref> Another one by S.K Ng, Thriyambakam Krishnan and G.J McLachlan in 1977.<ref>{{Citation |last1=Ng |first1=Shu Kay |title=The EM Algorithm |date=2011-12-21 |url=http://dx.doi.org/10.1007/978-3-642-21551-3_6 |work=Handbook of Computational Statistics |pages=139–172 |place=Berlin, Heidelberg |publisher=Springer Berlin Heidelberg |isbn=978-3-642-21550-6 |access-date=2022-10-15 |last2=Krishnan |first2=Thriyambakam |last3=McLachlan |first3=Geoffrey J.|doi=10.1007/978-3-642-21551-3_6 |s2cid=59942212 }}</ref> 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,<ref name="Sundberg1974">{{cite journal
 |last=Sundberg |first=Rolf
 |title=Maximum likelihood theory for incomplete data from an exponential family
 |journal=Scandinavian Journal of Statistics
 |volume=1 |year=1974 |issue=2 |pages=49–58
 |jstor=4615553 |mr= 381110
}}</ref><ref name="Sundberg1971">
Rolf 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.</ref><ref name="Sundberg1976">{{cite journal
 |last=Sundberg |first=Rolf
 |year=1976
 |title=An iterative method for solution of the likelihood equations for incomplete data from exponential families
 |journal=Communications in Statistics – Simulation and Computation
 |volume=5 |issue=1 |pages=55–64
 |doi=10.1080/03610917608812007 |mr=443190
}}</ref> following his collaboration with [[Biography:Per Martin-Löf|Per Martin-Löf]] and [[Biography:Anders Martin-Löf|Anders Martin-Löf]].<ref>See the acknowledgement by Dempster, Laird and Rubin on pages 3, 5 and 11.</ref><ref name="Martin-Löf1966">[[Biography:Per Martin-Löf|Per 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).</ref><ref name="Martin-Löf1970">[[Biography:Per Martin-Löf|Per 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. </ref><!-- * Martin-Löf, P. "Exact tests, confidence regions and estimates", with a discussion by [[Biography:A. W. F. Edwards|A. W. F. Edwards]], G. A. Barnard, D. A. Sprott, O. Barndorff-Nielsen, D. Basu and [[Rasch model|G. Rasch]]. ''Proceedings of Conference on Foundational Questions in Statistical Inference'' (Aarhus, 1973), pp. 121–138. Memoirs, No. 1, Dept. Theoret. Statist., Inst. Math., Univ. Aarhus, Aarhus, 1974. --><ref name="Martin-Löf1974a">Martin-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, [[Biography:Arthur P. Dempster|A. P. Dempster]], D. Basu, D. R. Cox, [[Biography:A. W. F. Edwards|A. W. F. Edwards]], D. A. Sprott, G. A. Barnard, O. Barndorff-Nielsen, J. D. Kalbfleisch and [[Rasch model|G. 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.</ref><ref name="Martin-Löf1974b">{{cite journal |last=Martin-Löf |first=Per |title=The notion of redundancy and its use as a quantitative measure of the discrepancy between a statistical hypothesis and a set of observational data |journal=Scand. J. Statist. |volume=1 |year=1974 |issue=1 |pages=3–18 }}</ref> 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 [[Biography:C. F. Jeff Wu|C. F. Jeff Wu]] in 1983.<ref name="Wu">
{{cite journal
|first=C. F. Jeff
|last=Wu
|title=On the Convergence Properties of the EM Algorithm
|journal=[[Annals of Statistics]]
|volume=11
|issue=1
|date=Mar 1983
|pages=95–103
|jstor=2240463
|doi=10.1214/aos/1176346060
|mr= 684867
|doi-access=free
}}</ref>
Wu's proof established the EM method's convergence also outside of the [[Exponential family|exponential family]], as claimed by Dempster–Laird–Rubin.<ref name="Wu" />

== Introduction ==
The EM algorithm is used to find (local) maximum likelihood parameters of a [[Statistical model|statistical model]] in cases where the equations cannot be solved directly.  Typically these models involve [[Social:Latent variable|latent variable]]s 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|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 [[Derivative|derivative]]s of the [[Likelihood function|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 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 a [[Saddle point|saddle point]].<ref name="Wu" /> 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 ==

=== The symbols ===
Given the [[Statistical model|statistical model]] which generates a set <math>\mathbf{X}</math> of observed data, a set of unobserved latent data or missing values <math>\mathbf{Z}</math>, and a vector of unknown parameters <math>\boldsymbol\theta</math>, along with a [[Likelihood function|likelihood function]] <math>L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}\mid\boldsymbol\theta)</math>, the maximum likelihood estimate (MLE) of the unknown parameters is determined by maximizing the [[Marginal likelihood|marginal likelihood]] of the observed data

:<math>L(\boldsymbol\theta; \mathbf{X}) = p(\mathbf{X}\mid\boldsymbol\theta) = \int  p(\mathbf{X},\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} = \int  p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol\theta) p(\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} </math>

However, this quantity is often intractable since <math>\mathbf{Z}</math> is unobserved and the distribution of <math>\mathbf{Z}</math> is unknown before attaining <math>\boldsymbol\theta</math>.

=== The EM algorithm ===
The EM algorithm seeks to find the MLE of the marginal likelihood by iteratively applying these two steps:
:''Expectation step (E step)'': Define <math>Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)})</math> as the [[Expected value|expected value]] of the log [[Likelihood function|likelihood function]] of <math>\boldsymbol\theta</math>, with respect to the current [[Conditional probability distribution|conditional distribution]] of <math>\mathbf{Z}</math> given <math>\mathbf{X}</math> and the current estimates of the parameters <math>\boldsymbol\theta^{(t)}</math>:
::<math>Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) = \operatorname{E}_{\mathbf{Z} \sim p(\cdot | \mathbf{X},\boldsymbol\theta^{(t)})}\left[ \log p (\mathbf{X},\mathbf{Z} | \boldsymbol\theta)  \right] \,</math>

:''Maximization step (M step)'': Find the parameters that maximize this quantity:
::<math>\boldsymbol\theta^{(t+1)} = \underset{\boldsymbol\theta}{\operatorname{arg\,max}} \ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \, </math>

More succinctly, we can write it as one equation:<math display="block">\boldsymbol\theta^{(t+1)} = \underset{\boldsymbol\theta}{\operatorname{arg\,max}} 
\operatorname{E}_{\mathbf{Z} \sim p(\cdot | \mathbf{X},\boldsymbol\theta^{(t)})}\left[ \log p (\mathbf{X},\mathbf{Z} | \boldsymbol\theta)  \right] \, </math>

=== Interpretation of the variables ===
The typical models to which EM is applied use <math>\mathbf{Z}</math> as a latent variable indicating membership in one of a set of groups:
#The observed data points <math>\mathbf{X}</math> may be discrete (taking values in a finite or countably infinite set) or continuous (taking values in an uncountably infinite set). Associated with each data point may be a vector of observations.
#The missing values (aka latent variables) <math>\mathbf{Z}</math> are discrete, drawn from a fixed number of values, and with one latent variable per observed unit.
#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 <math>\boldsymbol\theta</math> is known, usually the value of the latent variables <math>\mathbf{Z}</math> can be found by maximizing the log-likelihood over all possible values of <math>\mathbf{Z}</math>, either simply by iterating over <math>\mathbf{Z}</math> or through an algorithm such as the [[Viterbi algorithm]] for [[Hidden Markov model|hidden Markov model]]s.  Conversely, if we know the value of the latent variables <math>\mathbf{Z}</math>, we can find an estimate of the parameters <math>\boldsymbol\theta</math> 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 <math>\boldsymbol\theta</math> and <math>\mathbf{Z}</math> are unknown:
# First, initialize the parameters <math>\boldsymbol\theta</math> to some random values.
# Compute the probability of each possible value of <math>\mathbf{Z}</math> , given <math>\boldsymbol\theta</math>.
# Then, use the just-computed values of <math>\mathbf{Z}</math> to compute a better estimate for the parameters <math>\boldsymbol\theta</math>.
# Iterate steps 2 and 3 until convergence.
The algorithm as just described monotonically approaches a local minimum of the cost function.

== Properties ==
Although an EM iteration does increase the observed data (i.e., marginal) 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|metaheuristic]] approaches exist to escape a local maximum, such as random-restart [[Hill climbing|hill climbing]] (starting with several different random initial estimates <math>\boldsymbol\theta^{(t)}</math>), or applying [[Simulated annealing|simulated annealing]] methods.

EM is especially useful when the likelihood is an [[Exponential family|exponential family]], see Sundberg (2019, Ch. 8) for a comprehensive treatment:<ref>{{cite book |last1=Sundberg |first1=Rolf |title=Statistical Modelling by Exponential Families |date=2019 |publisher=Cambridge University Press |isbn=9781108701112}}</ref> the E step becomes the sum of expectations of [[Sufficient statistic|sufficient statistic]]s, and the M step involves maximizing a linear function. In such a case, it is usually possible to derive [[Closed-form expression|closed-form expression]] updates for each step, using the Sundberg formula<ref>{{cite book |last1=Laird |first1=Nan |title=Encyclopedia of Statistical Sciences |chapter=Sundberg formulas |chapter-url=https://doi.org/10.1002/0471667196.ess2643.pub2 |publisher=Wiley |date=2006|doi=10.1002/0471667196.ess2643.pub2 |isbn=0471667196 }}</ref> (proved and published by Rolf Sundberg, based on unpublished results of [[Biography:Per Martin-Löf|Per Martin-Löf]] and [[Biography:Anders Martin-Löf|Anders Martin-Löf]]).<ref name="Sundberg1971"/><ref name="Sundberg1976"/><ref name="Martin-Löf1966"/><ref name="Martin-Löf1970"/><ref name="Martin-Löf1974a"/><ref name="Martin-Löf1974b"/>

The EM method was modified to compute maximum a posteriori (MAP) 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|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 <math>Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)})</math> rather than directly improving <math>\log p(\mathbf{X}\mid\boldsymbol\theta)</math>.  Here it is shown that improvements to the former imply improvements to the latter.<ref name="Little1987">{{cite book |last1=Little |first1= Roderick J.A. |last2= Rubin |first2= Donald B. |title= Statistical Analysis with Missing Data |url=https://archive.org/details/statisticalanaly00litt |url-access=limited | series = Wiley Series in Probability and Mathematical Statistics |year= 1987 |publisher= John Wiley & Sons |location= New York |isbn= 978-0-471-80254-9 |pages= [https://archive.org/details/statisticalanaly00litt/page/n145 134]–136}}</ref>

For any <math>\mathbf{Z}</math> with non-zero probability <math>p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta)</math>, we can write
: <math>
\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:
: <math>
\begin{align}
\log p(\mathbf{X}\mid\boldsymbol\theta) &
= \sum_{\mathbf{Z}} p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) \log p(\mathbf{X},\mathbf{Z}\mid\boldsymbol\theta)
- \sum_{\mathbf{Z}} p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) \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>,
: <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
: <math>
\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
: <math>
\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|coordinate descent]].<ref name="neal1999">{{cite book|last1=Neal |first1=Radford |last2=Hinton |first2=Geoffrey |title=A view of the EM algorithm that justifies incremental, sparse, and other variants |journal=Learning in Graphical Models |editor=Michael I. Jordan |pages= 355–368 |publisher= MIT Press |location=Cambridge, MA |year=1999 |isbn=978-0-262-60032-3 |url=http://ftp.cs.toronto.edu/pub/radford/emk.pdf |access-date=2009-03-22}}</ref><ref name="hastie2001">{{cite book|last1=Hastie |first1= Trevor|author-link1=Trevor Hastie|last2=Tibshirani|first2=Robert|author-link2=Robert Tibshirani|last3=Friedman|first3=Jerome |year=2001 |title=The Elements of Statistical Learning |url=https://archive.org/details/elementsstatisti00thas_842 |url-access=limited |isbn=978-0-387-95284-0 |publisher=Springer |location=New York |chapter=8.5 The EM algorithm |pages=[https://archive.org/details/elementsstatisti00thas_842/page/n237 236]–243}}</ref> Consider the function:
:<math> F(q,\theta) := \operatorname{E}_q [ \log L (\theta ; x,Z) ] + H(q), </math> 
where ''q'' is an arbitrary probability distribution over the unobserved data ''z'' and ''H(q)'' is the [[Entropy (information theory)|entropy]] of the distribution ''q''. This function can be written as
:<math> 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>F</math>:
::<math> q^{(t)} = \operatorname{arg\,max}_q \ F(q,\theta^{(t)}) </math>
:''Maximization step'': Choose <math>\theta</math> to maximize <math>F</math>:
::<math> \theta^{(t+1)} = \operatorname{arg\,max}_\theta \ F(q^{(t)},\theta) </math>

== Applications ==
*EM is frequently used for parameter estimation of [[Mixed model|mixed model]]s,<ref>{{cite journal |doi=10.1080/01621459.1988.10478693 |title=Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data |journal=Journal of the American Statistical Association |volume=83 |issue=404 |pages=1014 |year=1988 |last1=Lindstrom |first1=Mary J |last2=Bates |first2=Douglas M }}</ref><ref>{{cite journal |doi=10.2307/1390614 |jstor=1390614 |title=Fitting Mixed-Effects Models Using Efficient EM-Type Algorithms |journal=Journal of Computational and Graphical Statistics |volume=9 |issue=1 |pages=78–98 |year=2000 |last1=Van Dyk |first1=David A }}</ref> notably in [[Biology:Quantitative genetics|quantitative genetics]].<ref>{{cite journal |doi=10.1111/anzs.12208 |title=A new REML (parameter expanded) EM algorithm for linear mixed models |journal=Australian & New Zealand Journal of Statistics |volume=59 |issue=4 |pages=433 |year=2017 |last1=Diffey |first1=S. M |last2=Smith |first2=A. B |last3=Welsh |first3=A. H |last4=Cullis |first4=B. R |doi-access=free |hdl=1885/211365 |hdl-access=free }}</ref>
*In [[Psychometrics|psychometrics]], EM is an important tool for estimating item parameters and latent abilities of [[Item response theory|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.{{Citation needed|date=November 2017}}
*The EM algorithm (and its faster variant [[Ordered subset expectation maximization|ordered subset expectation maximization]]) is also widely used in [[Medicine:Medical imaging|medical image]] reconstruction, especially in [[Physics:Positron emission tomography|positron emission tomography]], [[Biology:Single-photon emission computed tomography|single-photon emission computed tomography]], and x-ray computed tomography. See below for other faster variants of EM.
*In [[Engineering:Structural engineering|structural engineering]], the Structural Identification using Expectation Maximization (STRIDE)<ref>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</ref> algorithm is an output-only method for identifying natural vibration properties of a structural system using sensor data (see [[Physics:Operational Modal Analysis|Operational Modal Analysis]]).
*EM is also used for data clustering. In [[Natural language processing|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 grammar|probabilistic context-free grammar]]s.
*In the analysis of intertrade waiting times i.e. the time between subsequent trades in [[Finance:Share|shares of stock]] at a [[Finance:Stock exchange|stock exchange]] the EM algorithm has proved to be very useful.<ref>{{Cite journal|last1=Kreer|first1=Markus|last2=Kizilersu|first2=Ayse|last3=Thomas|first3=Anthony W.|date=2022|title=Censored expectation maximization algorithm for mixtures: Application to intertrade waiting times|journal=Physica A: Statistical Mechanics and Its Applications|volume=587|issue=1|pages=126456|doi=10.1016/j.physa.2021.126456|bibcode=2022PhyA..58726456K |s2cid=244198364 |issn=0378-4371|url=https://www.sciencedirect.com/science/article/pii/S0378437121007299}}</ref>

== 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>\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>\widehat{\sigma}^2_w =   \frac{1}{N} \sum_{k=1}^N {(\widehat{x}_{k+1}-\widehat{F}\widehat{{x}}_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>\widehat{F} = \frac{\sum_{k=1}^N {(\widehat{x}_{k+1}-\widehat{F} \widehat{x}_k)}^2}{\sum_{k=1}^N \widehat{x}_k^2}.</math>

The convergence of parameter estimates such as those above are well studied.<ref>{{Cite journal
|last1 = Einicke
|first1 =  G. A.
|last2 = Malos
|first2 = J. T.
|last3 = Reid
|first3 =  D. C.
|last4 = Hainsworth
|first4 = D. W.
|title = Riccati Equation and EM Algorithm Convergence for Inertial Navigation Alignment
|journal = IEEE Trans. Signal Process.
|volume = 57
|issue = 1
|pages = 370–375
|date=January 2009
|doi = 10.1109/TSP.2008.2007090
|bibcode = 2009ITSP...57..370E
|s2cid = 1930004
}}</ref><ref>{{Cite journal
|last1 = Einicke
|first1 =  G. A.
|last2 = Falco
|first2 =  G.
|last3 = Malos
|first3 = J. T.
|title = EM Algorithm State Matrix Estimation for Navigation
|journal = IEEE Signal Processing Letters
|volume = 17
|issue = 5
|pages = 437–440
|date=May 2010
|doi = 10.1109/LSP.2010.2043151
|bibcode = 2010ISPL...17..437E |s2cid = 14114266
}}</ref><ref>{{Cite journal
|last1 = Einicke
|first1 =  G. A.
|last2 = Falco
|first2 =  G.
|last3 = Dunn
|first3 = M. T.
|last4 = Reid
|first4 = D. C.
|title = Iterative Smoother-Based Variance Estimation
|journal = IEEE Signal Processing Letters
|volume = 19
|issue = 5
|pages = 275–278
|date=May 2012
|bibcode = 2012ISPL...19..275E
|doi = 10.1109/LSP.2012.2190278
|s2cid = 17476971
}}</ref><ref>{{Cite journal
|last = Einicke
|first =  G. A.
|title = Iterative Filtering and Smoothing of Measurements Possessing Poisson Noise
|journal = IEEE Transactions on Aerospace and Electronic Systems
|volume = 51
|issue = 3
|pages = 2205–2011
|date=Sep 2015
|doi = 10.1109/TAES.2015.140843
|bibcode = 2015ITAES..51.2205E
|s2cid = 32667132
}}</ref>

== 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 method]]s (Newton–Raphson).<ref>{{cite journal |first1= Mortaza |last1=Jamshidian |first2=Robert I. |last2=Jennrich|title=Acceleration of the EM Algorithm by using Quasi-Newton Methods |year=1997 |journal=Journal of the Royal Statistical Society, Series B |volume=59 |issue=2 |pages=569–587 |doi=10.1111/1467-9868.00083 |mr=1452026 |s2cid=121966443 }}</ref> 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".<ref>{{cite journal |doi=10.1093/biomet/85.4.755 |title=Parameter expansion to accelerate EM: The PX-EM algorithm |journal=Biometrika |volume=85 |issue=4 |pages=755–770 |year=1998 |last1=Liu |first1=C |citeseerx=10.1.1.134.9617 }}</ref>

''Expectation conditional maximization (ECM)'' replaces each M step with a sequence of conditional maximization (CM) steps in which each parameter ''θ''<sub>''i''</sub> is maximized individually, conditionally on the other parameters remaining fixed.<ref>{{cite journal|last1=Meng  |first1= Xiao-Li |last2=Rubin |first2=Donald B. |s2cid= 40571416 |title=Maximum likelihood estimation via the ECM algorithm: A general framework |year=1993 |journal=[[Biometrika]] |volume=80 |issue=2 |pages=267–278 |doi=10.1093/biomet/80.2.267 |mr=1243503}}</ref> Itself can be extended into the ''Expectation conditional maximization either (ECME)'' algorithm.<ref>{{cite journal |doi= 10.1093/biomet/81.4.633|jstor=2337067 |title=The ECME Algorithm: A Simple Extension of EM and ECM with Faster Monotone Convergence |journal=Biometrika |volume=81 |issue=4 |pages=633 |year=1994 |last1=Liu |first1=Chuanhai |last2=Rubin |first2=Donald B }}</ref>

This idea is further extended in ''generalized expectation maximization (GEM)'' 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.<ref name="neal1999"/> GEM is further developed in a distributed environment and shows promising results.<ref>
{{cite journal
 |author1=Jiangtao Yin |author2=Yanfeng Zhang |author3=Lixin Gao |title=Accelerating Expectation–Maximization Algorithms with Frequent Updates
 |journal=Proceedings of the IEEE International Conference on Cluster Computing
 |year=2012
 |url=http://rio.ecs.umass.edu/mnilpub/papers/cluster2012-yin.pdf
}}
</ref>

It is also possible to consider the EM algorithm as a subclass of the '''[[MM algorithm|MM]]''' (Majorize/Minimize or Minorize/Maximize, depending on context) algorithm,<ref>Hunter DR and Lange K (2004), [http://www.stat.psu.edu/~dhunter/papers/mmtutorial.pdf A Tutorial on MM Algorithms], The American Statistician, 58: 30–37</ref> 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<ref>
{{cite journal
|last=Matsuyama |first=Yasuo
|title=The α-EM algorithm: Surrogate likelihood maximization using α-logarithmic information measures
|journal=IEEE Transactions on Information Theory
|volume=49 | year=2003 |pages=692–706 |issue=3
|doi=10.1109/TIT.2002.808105
}}
</ref>
which contains the log-EM algorithm as its subclass. Thus, the α-EM algorithm by [[Biography:Yasuo Matsuyama|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.
<ref>
{{cite journal
|last=Matsuyama |first=Yasuo
|title=Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs
|journal=International Joint Conference on Neural Networks
| year=2011 |pages=808–816
}}
</ref>

== Relation to variational Bayes methods ==
EM is a partially non-Bayesian, maximum likelihood method.  Its final result gives a [[Probability distribution|probability distribution]] over the latent variables (in the Bayesian style) together with a point estimate for ''θ'' (either a [[Maximum likelihood estimation|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, ''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 connection|affine connection]]s, called the e-connection and the m-connection; the [[Kullback–Leibler divergence]] can also be understood in these terms.

== Examples ==
=== Gaussian mixture === <!--This section is linked from [[Matrix calculus]] -->
[[File:ClusterAnalysis_Mouse.svg|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 diagram|Voronoi]]-cells. The cluster center is indicated by the lighter, bigger symbol.]]

[[File:Em old faithful.gif|thumb|240px|An animation demonstrating the EM algorithm fitting a two component Gaussian [[Mixture model|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 model|mixture]] of two [[Multivariate normal distribution|multivariate normal distribution]]s 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.<ref name="hastie2001"/>
: <math>X_i \mid(Z_i = 1) \sim \mathcal{N}_d(\boldsymbol{\mu}_1,\Sigma_1)</math> and <math>X_i \mid(Z_i = 2) \sim \mathcal{N}_d(\boldsymbol{\mu}_2,\Sigma_2),</math>
where
: <math>\operatorname{P} (Z_i = 1 ) = \tau_1 \, </math> and <math>\operatorname{P} (Z_i=2) = \tau_2 = 1-\tau_1.</math>

The aim is to estimate the unknown parameters representing the ''mixing'' value between the Gaussians and the means and covariances of each:
: <math>\theta = \big( \boldsymbol{\tau},\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma_1,\Sigma_2 \big),</math>
where the incomplete-data likelihood function is
: <math>L(\theta;\mathbf{x}) =  \prod_{i=1}^n \sum_{j=1}^2 \tau_j \ f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j),</math>

and the complete-data likelihood function is
: <math>L(\theta;\mathbf{x},\mathbf{z}) = p(\mathbf{x},\mathbf{z} \mid \theta) = \prod_{i=1}^n  \prod_{j=1}^2  \ [f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j) \tau_j] ^{\mathbb{I}(z_i=j)},</math>

or

: <math>L(\theta;\mathbf{x},\mathbf{z}) = \exp \left\{ \sum_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big] \right\},</math>

where <math>\mathbb{I}</math> is an [[Indicator function|indicator function]] and <math>f</math> is the [[Physics:Probability density function|probability density function]] of a multivariate normal.

In the last equality, for each {{math|''i''}}, one indicator <math>\mathbb{I}(z_i=j)</math> 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 ''θ''<sup>(''t'')</sup>, the conditional distribution of the ''Z''<sub>''i''</sub> is determined by Bayes theorem to be the proportional height of the normal [[Physics:Probability density function|density]] weighted by ''τ'':
: <math>T_{j,i}^{(t)} := \operatorname{P}(Z_i=j \mid X_i=\mathbf{x}_i ;\theta^{(t)}) = \frac{\tau_j^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)})}{\tau_1^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_1^{(t)},\Sigma_1^{(t)}) + \tau_2^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_2^{(t)},\Sigma_2^{(t)})}.</math>

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:
: <math>\begin{align}Q(\theta\mid\theta^{(t)})
&= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} [\log L(\theta;\mathbf{x},\mathbf{Z}) ] \\
&= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} [\log \prod_{i=1}^{n}L(\theta;\mathbf{x}_i,Z_i) ] \\
&= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} [\sum_{i=1}^n \log L(\theta;\mathbf{x}_i,Z_i) ] \\
&= \sum_{i=1}^n\operatorname{E}_{Z_i\mid X_i=x_i;\mathbf{\theta}^{(t)}} [\log L(\theta;\mathbf{x}_i,Z_i) ] \\
&= \sum_{i=1}^n \sum_{j=1}^2 P(Z_i =j \mid X_i = \mathbf{x}_i; \theta^{(t)}) \log L(\theta_j;\mathbf{x}_i,j) \\
&= \sum_{i=1}^n \sum_{j=1}^2 T_{j,i}^{(t)} \big[ \log \tau_j  -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big].
\end{align}</math>
The expectation of <math>\log L(\theta;\mathbf{x}_i,Z_i)</math> inside the sum is taken with respect to the probability density function <math>P(Z_i \mid X_i = \mathbf{x}_i; \theta^{(t)})</math>, which might be different for each  <math>\mathbf{x}_i</math> of the training set. Everything in the E step is known before the step is taken except <math>T_{j,i}</math>, 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''(''θ''&nbsp;|&nbsp;''θ''<sup>(''t'')</sup>) being quadratic in form means that determining the maximizing values of ''θ'' is relatively straightforward. Also, ''τ'', ('''μ'''<sub>1</sub>,''Σ''<sub>1</sub>) and ('''μ'''<sub>2</sub>,''Σ''<sub>2</sub>) may all be maximized independently since they all appear in separate linear terms.

To begin, consider ''τ'', which has the constraint ''τ''<sub>1</sub> + ''τ''<sub>2</sub>=1:
: <math>\begin{align}\boldsymbol{\tau}^{(t+1)}
&= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}}\  Q(\theta \mid \theta^{(t)} ) \\
&= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}} \ \left\{ \left[  \sum_{i=1}^n T_{1,i}^{(t)} \right] \log \tau_1 + \left[  \sum_{i=1}^n T_{2,i}^{(t)} \right] \log \tau_2  \right\}.
\end{align}</math>
This has the same form as the MLE for the [[Binomial distribution|binomial distribution]], so
: <math>\tau^{(t+1)}_j = \frac{\sum_{i=1}^n T_{j,i}^{(t)}}{\sum_{i=1}^n (T_{1,i}^{(t)} + T_{2,i}^{(t)} ) } = \frac{1}{n} \sum_{i=1}^n T_{j,i}^{(t)}.</math>

For the next estimates of ('''μ'''<sub>1</sub>,Σ<sub>1</sub>):
: <math>\begin{align}(\boldsymbol{\mu}_1^{(t+1)},\Sigma_1^{(t+1)})
&= \underset{\boldsymbol{\mu}_1,\Sigma_1} \operatorname{arg\,max}\ Q(\theta \mid \theta^{(t)} ) \\
&= \underset{\boldsymbol{\mu}_1,\Sigma_1} \operatorname{arg\,max}\ \sum_{i=1}^n T_{1,i}^{(t)} \left\{ -\tfrac{1}{2} \log |\Sigma_1| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_1)^\top\Sigma_1^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_1) \right\}
\end{align}.</math>
This has the same form as a weighted MLE for a normal distribution, so
: <math>\boldsymbol{\mu}_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{1,i}^{(t)}} </math> and <math>\Sigma_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)})^\top }{\sum_{i=1}^n T_{1,i}^{(t)}} </math>
and, by symmetry,
: <math>\boldsymbol{\mu}_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{2,i}^{(t)}} </math> and <math>\Sigma_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)})^\top }{\sum_{i=1}^n T_{2,i}^{(t)}}.</math>

==== Termination ====

Conclude the iterative process if <math> E_{Z\mid\theta^{(t)},\mathbf{x}}[\log L(\theta^{(t)};\mathbf{x},\mathbf{Z})] \leq E_{Z\mid\theta^{(t-1)},\mathbf{x}}[\log L(\theta^{(t-1)};\mathbf{x},\mathbf{Z})] + \varepsilon</math> for <math> \varepsilon </math> below some preset threshold.
==== Generalization ====
The algorithm illustrated above can be generalized for mixtures of more than two [[Multivariate normal distribution|multivariate normal distribution]]s.

=== Truncated and censored regression ===
The EM algorithm has been implemented in the case where an underlying [[Linear regression|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.<ref name=Wolynetz>{{cite journal |title= Maximum likelihood estimation in a linear model from confined and censored normal data |last1= Wolynetz |first1= M.S. |journal= Journal of the Royal Statistical Society, Series C |year= 1979 |volume= 28 |issue= 2 |pages= 195–206 |doi= 10.2307/2346749|jstor= 2346749 }}</ref>  Special cases of this model include censored or truncated observations from one [[Normal distribution|normal distribution]].<ref name=Wolynetz/>

== 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''<ref>{{Cite journal|last=Pearson|first=Karl|date=1894|title=Contributions to the Mathematical Theory of Evolution|journal=Philosophical Transactions of the Royal Society of London A|volume=185|pages=71–110|issn=0264-3820|jstor=90667|doi=10.1098/rsta.1894.0003|bibcode=1894RSPTA.185...71P|doi-access=free}}</ref> or the so-called ''spectral techniques''<ref>{{Cite journal|last1=Shaban|first1=Amirreza|last2=Mehrdad|first2=Farajtabar|last3=Bo|first3=Xie|last4=Le|first4=Song|last5=Byron|first5=Boots|date=2015|title=Learning Latent Variable Models by Improving Spectral Solutions with Exterior Point Method|url=https://www.cc.gatech.edu/~bboots3/files/SpectralExteriorPoint-NIPSWorkshop.pdf|journal=UAI|pages=792–801}}</ref><ref>{{Cite book|title=Local Loss Optimization in Operator Models: A New Insight into Spectral Learning|last=Balle, Borja Quattoni, Ariadna Carreras, Xavier|date=2012-06-27|oclc=815865081}}</ref>{{Citation needed|date=April 2019}}. Moment-based approaches to learning the parameters of a probabilistic model are of increasing interest recently{{when|date=January 2022}} 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{{Citation needed|date=April 2019}}.

== See also ==
* mixture distribution
* compound distribution
* [[Density estimation|density estimation]]
* [[Principal component analysis]]
* [[Physics:Total absorption spectroscopy|total absorption spectroscopy]]
* The EM algorithm can be viewed as a special case of the [[MM algorithm|majorize-minimization (MM) algorithm]].<ref>{{cite web|last=Lange|first=Kenneth|title=The MM Algorithm|url=http://www.stat.berkeley.edu/~aldous/Colloq/lange-talk.pdf}}</ref>

== References ==
{{Reflist}}

== Further reading ==
* {{cite book |first1=Robert |last1=Hogg |first2=Joseph |last2=McKean |first3=Allen |last3=Craig |title=Introduction to Mathematical Statistics |pages=359–364 |location=Upper Saddle River, NJ |publisher=Pearson Prentice Hall |year=2005 }}
* {{cite journal |citeseerx= 10.1.1.9.9735 |title= The Expectation Maximization Algorithm|first=Frank |last= Dellaert |year= 2002}} gives an easier explanation of EM algorithm as to lowerbound maximization.
* {{cite book
 |last1 = Bishop |first1 = Christopher M.
 |title = Pattern Recognition and Machine Learning
 |year = 2006
 |publisher = Springer
 |isbn = 978-0-387-31073-2
 }}
* {{cite journal |doi=10.1561/2000000034 |title=Theory and Use of the EM Algorithm |journal=Foundations and Trends in Signal Processing |volume=4 |issue=3 |pages=223–296 |first1=M. R. |last1=Gupta |first2=Y. |last2=Chen |year=2010 |citeseerx=10.1.1.219.6830 }} A well-written short book on EM, including detailed derivation of EM for GMMs, HMMs, and Dirichlet. 
* {{cite journal |citeseerx= 10.1.1.28.613 |title= A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models|first=Jeff |last= Bilmes |year= 1998}} includes a simplified derivation of the EM equations for Gaussian Mixtures and Gaussian Mixture Hidden Markov Models.
* {{cite book |first1=Geoffrey J. |last1=McLachlan |first2=Thriyambakam |last2=Krishnan |title=The EM Algorithm and Extensions |location=Hoboken |publisher=Wiley |edition=2nd |year=2008 |isbn=978-0-471-20170-0 }}

== External links ==
* Various 1D, 2D and 3D [http://wiki.stat.ucla.edu/socr/index.php/SOCR_EduMaterials_Activities_2D_PointSegmentation_EM_Mixture 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.
* [https://github.com/l-/CommonDataAnalysis Class hierarchy] in [[C++]] (GPL) including Gaussian Mixtures
* [http://www.inference.phy.cam.ac.uk/mackay/itila/ 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 ''k''-means algorithm, and emphasizes the variational view of the EM algorithm, as described in Chapter 33.7 of version 7.2 (fourth edition).
* [http://www.cse.buffalo.edu/faculty/mbeal/papers/beal03.pdf 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  ([http://www.cse.buffalo.edu/faculty/mbeal/thesis/index.html chapters]).
* [http://www.seanborman.com/publications/EM_algorithm.pdf The Expectation Maximization Algorithm: A short tutorial], A self-contained derivation of the EM Algorithm by Sean Borman.
* [http://pages.cs.wisc.edu/~jerryzhu/cs838/EM.pdf The EM Algorithm], by Xiaojin Zhu.
* [https://arxiv.org/abs/1105.1476 EM algorithm and variants: an informal tutorial] by Alexis Roche.  A concise and very clear description of EM and many interesting variants.

{{DEFAULTSORT:Expectation-maximization Algorithm}}
[[Category:Estimation methods]]
[[Category:Machine learning algorithms]]
[[Category:Missing data]]
[[Category:Statistical algorithms]]
[[Category:Optimization algorithms and methods]]
[[Category:Cluster analysis algorithms]]

{{Sourceattribution|Expectation–maximization algorithm}}