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. 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.

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. 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. Given the statistical model which generates a set \mathbf{X} of observed data, a set of unobserved latent data or missing values \mathbf{Z} , and a vector of unknown parameters \boldsymbol\theta , along with a likelihood function L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}\mid\boldsymbol\theta) , the maximum likelihood estimate (MLE) of the unknown parameters is determined by maximizing the marginal likelihood of the observed data \begin{align} 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} \end{align} However, this quantity is often intractable since \mathbf{Z} is unobserved and the distribution of \mathbf{Z} is unknown before attaining \boldsymbol\theta . {{block indent | em = 1.5 | text = Expectation step (E step): Define Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) as the expected value of the log likelihood function of , with respect to the current conditional distribution of \mathbf{Z} given \mathbf{X} and the current estimates of the parameters : 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] = \int \log p (\mathbf X, \mathbf Z | \boldsymbol \theta) \, p (\mathbf Z | \mathbf X , \boldsymbol \theta^{(t)}) \, d \mathbf Z \, Maximization step (M step): Find the parameters that maximize this quantity: \boldsymbol\theta^{(t+1)} = \mathop{\arg\max}_{\boldsymbol\theta} Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \, }} More succinctly, we can write it as one equation: \boldsymbol\theta^{(t+1)} = \mathop{\arg\max}_{\boldsymbol\theta} \operatorname{E}_{\mathbf{Z} \sim p(\cdot | \mathbf{X},\boldsymbol\theta^{(t)})} \left[ \log p (\mathbf{X},\mathbf{Z} | \boldsymbol\theta) \right] \, \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). We take the expectation over possible values of the unknown data \mathbf{Z} under the current parameter estimate \theta^{(t)} by multiplying both sides by p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) and summing (or integrating) over \mathbf{Z} . The left-hand side is the expectation of a constant, so we get: \begin{align} \log p(\mathbf{X}\mid\boldsymbol\theta) & This last equation holds for every value of \boldsymbol\theta including \boldsymbol\theta = \boldsymbol\theta^{(t)} , \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)}) = Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}), 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)}). However, Gibbs' inequality tells us that H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \ge H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) , 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)}). In words, choosing \boldsymbol\theta to improve Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) causes \log p(\mathbf{X}\mid\boldsymbol\theta) to improve at least as much. F(q,\theta) := \operatorname{E}_q [ \log L (\theta ; x,Z) ] + H(q), where is an arbitrary probability distribution over the unobserved data and is the entropy of the distribution . This function can be written as F(q,\theta) = -D_{\mathrm{KL}}\big(q \parallel p_{Z\mid X}(\cdot\mid x;\theta ) \big) + \log L(\theta;x), where p_{Z\mid X}(\cdot\mid x;\theta ) is the conditional distribution of the unobserved data given the observed data x and D_{KL} is the Kullback–Leibler divergence. Then the steps in the EM algorithm may be viewed as: Expectation step: Choose q to maximize : q^{(t)} = \mathop{\arg\max}_q F{\left(q,\theta^{(t)}\right)} Maximization step: Choose \theta to maximize : \theta^{(t+1)} = \mathop{\arg\max}_\theta F{\left(q^{(t)},\theta\right)} \widehat{\sigma}^2_v = \frac{1}{N} \sum_{k=1}^N {(z_k-\widehat{x}_k)}^2, where \widehat{x}_k are scalar output estimates calculated by a filter or a smoother from N scalar measurements z_k . 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 \widehat{\sigma}^2_w = \frac{1}{N} \sum_{k=1}^N {(\widehat{x}_{k+1}-\widehat{F}\widehat _k)}^2, where \widehat{x}_k and \widehat{x}_{k+1} are scalar state estimates calculated by a filter or a smoother. The updated model coefficient estimate is obtained via \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}. The convergence of parameter estimates such as those above are well studied. Let \mathbf{x} = (\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n) be a sample of n independent observations from a mixture of two multivariate normal distributions of dimension d , and let \mathbf{z} = (z_1,z_2,\ldots,z_n) be the latent variables that determine the component from which the observation originates. \begin{align} X_i \mid(Z_i = 1) &\sim \mathcal{N}_d(\boldsymbol{\mu}_1,\Sigma_1), \\ X_i \mid(Z_i = 2) &\sim \mathcal{N}_d(\boldsymbol{\mu}_2,\Sigma_2), \end{align} where \operatorname{P} (Z_i = 1) = \tau_1 \, \quad \text{and} \quad \operatorname{P} (Z_i=2) = \tau_2 = 1-\tau_1. The aim is to estimate the unknown parameters representing the mixing value between the Gaussians and the means and covariances of each: \theta = \big( \boldsymbol{\tau},\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma_1,\Sigma_2 \big), where the incomplete-data likelihood function is L(\theta;\mathbf{x}) = \prod_{i=1}^n \sum_{j=1}^2 \tau_j \ f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j), and the complete-data likelihood function is \begin{align} L(\theta;\mathbf{x},\mathbf{z}) &= p(\mathbf{x},\mathbf{z} \mid \theta) \\ &= \prod_{i=1}^n \prod_{j=1}^2 \left[f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j) \tau_j\right] ^{\mathbb{I}(z_i=j)}, \end{align} or \log L(\theta;\mathbf{x},\mathbf{z}) = \sum_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \left[ \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) \right] , where \mathbb{I} is an indicator function and f is the probability density function of a multivariate normal. In the last equality, for each , one indicator \mathbb{I}(z_i=j) is equal to zero, and one indicator is equal to one. The inner sum thus reduces to one term. \begin{align} T_{j,i}^{(t)} :={}& \operatorname{P}(Z_i=j \mid X_i=\mathbf{x}_i ;\theta^{(t)}) \\ ={}& \frac{\tau_j^{(t)} \, f{\left(\mathbf{x}_i;\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)}\right)}}{\tau_1^{(t)} \, f{\left(\mathbf{x}_i;\boldsymbol{\mu}_1^{(t)},\Sigma_1^{(t)}\right)} + \tau_2^{(t)} \, f{\left(\mathbf{x}_i;\boldsymbol{\mu}_2^{(t)},\Sigma_2^{(t)}\right)}}. \end{align} These are called the "membership probabilities", which are normally considered the output of the E step (although this is not the function of below). This E step corresponds with setting up this function for : \begin{align} Q(\theta\mid\theta^{(t)}) &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} \left[\log L(\theta;\mathbf{x},\mathbf{Z}) \right] \\ &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} \left[\log \prod_{i=1}^{n}L(\theta;\mathbf{x}_i,Z_i) \right] \\ &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X}=\mathbf{x};\mathbf{\theta}^{(t)}} \left[\sum_{i=1}^n \log L(\theta;\mathbf{x}_i,Z_i) \right] \\ &= \sum_{i=1}^n\operatorname{E}_{Z_i\mid X_i=x_i;\mathbf{\theta}^{(t)}} \left[\log L(\theta;\mathbf{x}_i,Z_i) \right] \\ &= \sum_{i=1}^n \sum_{j=1}^2 P(Z_i =j \mid X_i = \mathbf{x}_i; \theta^{(t)}) \log L(\theta;\mathbf{x}_i,j) \\ &= \sum_{i=1}^n \sum_{j=1}^2 T_{j,i}^{(t)} \left[ \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) \right]. \end{align} The expectation of \log L(\theta;\mathbf{x}_i,Z_i) inside the sum is taken with respect to the probability density function P(Z_i \mid X_i = \mathbf{x}_i; \theta^{(t)}) , which might be different for each \mathbf{x}_i of the training set. Everything in the E step is known before the step is taken except T_{j,i} , 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. To begin, consider \tau , which has the constraint \tau_1+\tau_2=1 : \begin{align} \boldsymbol{\tau}^{(t+1)} &= \mathop{\arg\max}_\boldsymbol{\tau} Q{\left(\theta \mid \theta^{(t)} \right)} \\ &= \mathop{\arg\max}_\boldsymbol{\tau} \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} This has the same form as the maximum likelihood estimate for the binomial distribution, so \tau^{(t+1)}_j = \frac{\sum\limits_{i=1}^n T_{j,i}^{(t)}}{\sum\limits_{i=1}^n \left(T_{1,i}^{(t)} + T_{2,i}^{(t)} \right) } = \frac{1}{n} \sum_{i=1}^n T_{j,i}^{(t)}. For the next estimates of (\boldsymbol{\mu}_1,\Sigma_1) : \begin{align} \left(\boldsymbol{\mu}_1^{(t+1)},\Sigma_1^{(t+1)}\right) &= \mathop{\arg\max}_{\boldsymbol{\mu}_1,\Sigma_1} Q{\left(\theta \mid \theta^{(t)} \right)} \\ &= \mathop{\arg\max}_{\boldsymbol{\mu}_1,\Sigma_1} \sum_{i=1}^n T_{1,i}^{(t)} \left\{ -\tfrac{1}{2} \log \left|\Sigma_1\right| - \tfrac{1}{2} \left(\mathbf{x}_i-\boldsymbol{\mu}_1\right)^\top \Sigma_1^{-1} \left(\mathbf{x}_i-\boldsymbol{\mu}_1\right) \right\}. \end{align} This has the same form as a weighted maximum likelihood estimate for a normal distribution, so \begin{align} \boldsymbol{\mu}_1^{(t+1)} &= \frac{\sum\limits_{i=1}^n T_{1,i}^{(t)} \mathbf{x}_i}{\sum\limits_{i=1}^n T_{1,i}^{(t)}}, \\[1ex] \Sigma_1^{(t+1)} &= \frac{\sum\limits_{i=1}^n T_{1,i}^{(t)} \left(\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}\right) \left(\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}\right)^\top }{\sum\limits_{i=1}^n T_{1,i}^{(t)}} \end{align} and, by symmetry, \begin{align} \boldsymbol{\mu}_2^{(t+1)} &= \frac{\sum\limits_{i=1}^n T_{2,i}^{(t)} \mathbf{x}_i}{\sum\limits_{i=1}^n T_{2,i}^{(t)}}, \\[1ex] \Sigma_2^{(t+1)} &= \frac{\sum\limits_{i=1}^n T_{2,i}^{(t)} \left(\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)}\right) \left(\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)}\right)^\top }{\sum\limits_{i=1}^n T_{2,i}^{(t)}}. \end{align}

Text derived from https://en.wikipedia.org/w/index.php?oldid=1369278721 under CC BY-SA 4.0; modifications have been made.