Алгоритм ожидания – максимизации

В статистике , ожидание максимизации ( ЭМ ) Алгоритм представляет собой итерационный метод , чтобы найти (локальное) максимальное правдоподобие или максимальный апостериорный (MAP) оценки параметров в статистических моделях , где модель зависит от ненаблюдаемого латентных переменных . Итерация EM чередуется между выполнением этапа ожидания (E), который создает функцию для ожидания логарифмической вероятности, оцениваемой с использованием текущей оценки параметров, и этапом максимизации (M), на котором вычисляются параметры, максимизирующие ожидаемое значение log. вероятность найти наШаг E. Эти оценки параметров затем используются для определения распределения скрытых переменных на следующем этапе E.

EM-кластеризация данных об извержениях Old Faithful . Случайная исходная модель (которая из-за разных масштабов осей кажется двумя очень плоскими и широкими сферами) соответствует наблюдаемым данным. На первых итерациях модель существенно меняется, но затем сходится к двум режимам гейзера . Визуализируется с помощью ELKI .

История [ править ]

Алгоритм EM был объяснен и получил свое название в классической статье 1977 года Артура Демпстера , Нэн Лэрд и Дональда Рубина . ^[1] Они указали, что этот метод «много раз предлагался при особых обстоятельствах» более ранними авторами. Одним из первых является метод подсчета генов для оценки частот аллелей Седрика Смита . ^[2] Очень подробное рассмотрение метода ЭМ для экспоненциальных семейств было опубликовано Рольфом Сундбергом в его диссертации и нескольких статьях ^{[3] [4] [5]} после его сотрудничества с Пером Мартином-Лёфом и Андерсом Мартином-Лёфом . ^{[6] [7][8] [9] [10] [11] [12]} Статья Демпстера – Лэрда – Рубина 1977 г. обобщила метод и наметила анализ сходимости для более широкого класса проблем. В статье Демпстера – Лэрда – Рубина ЭМ-метод определен как важный инструмент статистического анализа.

Анализ сходимости алгоритма Демпстера – Лэрда – Рубина был ошибочным, и правильный анализ сходимости был опубликован CF Jeff Wu в 1983 году. ^[13] Доказательство Ву установило сходимость метода EM вне экспоненциального семейства , как утверждал Демпстер – Лейрд– Вбивать в голову. ^[13]

Введение [ править ]

Алгоритм EM используется для поиска (локальных) параметров максимального правдоподобия статистической модели в случаях, когда уравнения не могут быть решены напрямую. Обычно эти модели включают скрытые переменные в дополнение к неизвестным параметрам и известным данным наблюдений. То есть либо среди данных существуют отсутствующие значения , либо модель может быть сформулирована более просто, допуская существование дополнительных ненаблюдаемых точек данных. Например, модель смеси можно описать более просто, если предположить, что каждая наблюдаемая точка данных имеет соответствующую ненаблюдаемую точку данных или скрытую переменную, определяющую компонент смеси, которому принадлежит каждая точка данных.

Нахождение решения с максимальным правдоподобием обычно требует взятия производных от функции правдоподобия по всем неизвестным значениям, параметрам и скрытым переменным и одновременного решения полученных уравнений. В статистических моделях со скрытыми переменными это обычно невозможно. Вместо этого результатом обычно является набор взаимосвязанных уравнений, в которых решение параметров требует значений скрытых переменных и наоборот, но замена одного набора уравнений в другой дает неразрешимое уравнение.

Алгоритм EM исходит из наблюдения, что существует способ решить эти две системы уравнений численно. Можно просто выбрать произвольные значения для одного из двух наборов неизвестных, использовать их для оценки второго набора, затем использовать эти новые значения, чтобы найти лучшую оценку первого набора, а затем продолжать чередовать два, пока оба результата не получат сходятся к неподвижным точкам. Не очевидно, что это сработает, но можно доказать, что в данном контексте это так, и что производная вероятности (произвольно близка) к нулю в этой точке, что, в свою очередь, означает, что точка является либо максимумом, либо седловой точки . ^[13] Как правило, могут возникать множественные максимумы без гарантии, что глобальный максимум будет найден. Некоторые вероятности также имеютособенности в них, т. е. бессмысленные максимумы. Например, одно из решений, которое может найти EM в смешанной модели, включает в себя установку одного из компонентов на нулевую дисперсию и параметр среднего значения для того же компонента, равный одной из точек данных.

Описание [ править ]

Учитывая статистическую модель , которая генерирует набор наблюдаемых данных, набор ненаблюдаемых латентных данных или пропущенные значения , и вектор неизвестных параметров , а также функция правдоподобия , то оценка максимального правдоподобия (MLE) неизвестных параметров определяются путем максимизации предельная вероятность наблюдаемых данных${\displaystyle \mathbf {X} }$ ${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$ ${\displaystyle L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )=p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})}$

${\displaystyle L({\boldsymbol {\theta }};\mathbf {X} )=p(\mathbf {X} \mid {\boldsymbol {\theta }})=\int p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})\,d\mathbf {Z} }$

Однако с этой величиной часто трудно справиться (например, если это последовательность событий, так что количество значений растет экспоненциально с увеличением длины последовательности, точное вычисление суммы будет чрезвычайно трудным).${\displaystyle \mathbf {Z} }$

Алгоритм EM пытается найти MLE предельного правдоподобия, итеративно применяя эти два шага:

Шаг Ожидание (Е шаг) : Определить , как ожидаемое значение из журнала функции правдоподобия из , по отношению к текущему условного распределения из дается и текущих оценок параметров :${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {X} }$${\displaystyle {\boldsymbol {\theta }}^{(t)}}$

${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)}}\left[\log L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )\right]\,}$

Шаг максимизации (шаг M) : найдите параметры, которые максимизируют это количество:

${\displaystyle {\boldsymbol {\theta }}^{(t+1)}={\underset {\boldsymbol {\theta }}{\operatorname {arg\,max} }}\ Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\,}$

Типичные модели, к которым применяется EM, используют как скрытую переменную, указывающую на принадлежность к одной из множества групп:${\displaystyle \mathbf {Z} }$

1. Наблюдаемые точки данных могут быть дискретными (принимающими значения в конечном или счетно бесконечном множестве) или непрерывными (принимающими значения в несчетно бесконечном множестве). С каждой точкой данных может быть связан вектор наблюдений.${\displaystyle \mathbf {X} }$
2. В пропущенные значения ( так называемые латентные переменные ) являются дискретными , взяты из фиксированного числа значений, и с одной скрытой переменной в наблюдаемой единицы.${\displaystyle \mathbf {Z} }$
3. Параметры являются непрерывными и бывают двух видов: параметры, связанные со всеми точками данных, и параметры, связанные с конкретным значением скрытой переменной (т. Е. Связанные со всеми точками данных, у которых соответствующая скрытая переменная имеет это значение).

Однако ЭМ можно применять к другим типам моделей.

Мотив таков. Если значение параметров известно, обычно значение скрытых переменных можно найти путем максимизации логарифма правдоподобия по всем возможным значениям , либо просто путем повторения, либо с помощью такого алгоритма, как алгоритм Баума – Велча для скрытого Маркова. модели . И наоборот, если мы знаем значение скрытых переменных , мы можем довольно легко найти оценку параметров , обычно просто группируя наблюдаемые точки данных в соответствии со значением связанной скрытой переменной и усредняя значения или некоторую функцию от значения баллов в каждой группе. Это предполагает итерационный алгоритм в случае, когда и и${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle {\boldsymbol {\theta }}}$${\displaystyle \mathbf {Z} }$ неизвестны:

1. Сначала инициализируйте параметры некоторыми случайными значениями.${\displaystyle {\boldsymbol {\theta }}}$
2. Вычислите вероятность каждого возможного значения , данного .${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$
3. Затем используйте только что вычисленные значения, чтобы вычислить более точную оценку параметров .${\displaystyle \mathbf {Z} }$${\displaystyle {\boldsymbol {\theta }}}$
4. Повторите шаги 2 и 3 до схождения.

Алгоритм, как только что описанный, монотонно приближается к локальному минимуму функции стоимости.

Свойства [ править ]

Говорить о шаге ожидания (E) немного неправильно . То , что вычисляются на первом этапе являются фиксированными, данные , зависящие от параметров функции Q . Как только параметры Q известны, они полностью определяются и максимизируются на втором (M) шаге EM-алгоритма.

Хотя итерация EM увеличивает наблюдаемые данные (т. Е. Предельную) функцию правдоподобия, нет гарантии, что последовательность сходится к оценщику максимального правдоподобия . Для мультимодальных распределений это означает, что алгоритм EM может сходиться к локальному максимуму наблюдаемой функции правдоподобия данных в зависимости от начальных значений. Существует множество эвристических или метаэвристических подходов, позволяющих избежать локального максимума, например, восхождение на холм со случайным перезапуском (начиная с нескольких различных случайных начальных оценок θ ^{( t )} ) или применение методов моделирования отжига .

EM особенно полезен, когда вероятность является экспоненциальным семейством : шаг E становится суммой ожиданий достаточной статистики , а шаг M включает максимизацию линейной функции. В таком случае обычно можно получить обновления выражений в закрытой форме для каждого шага, используя формулу Сундберга (опубликованную Рольфом Сундбергом с использованием неопубликованных результатов Пера Мартина-Лёфа и Андерса Мартина-Лёфа ). ^{[4] [5] [8] [9] [10] [11] [12]}

Метод EM был изменен для вычисления максимальных апостериорных оценок (MAP) для байесовского вывода в исходной статье Демпстера, Лэрда и Рубина.

Существуют и другие методы для поиска оценок максимального правдоподобия, такие как градиентный спуск , сопряженный градиент или варианты алгоритма Гаусса – Ньютона . В отличие от EM, такие методы обычно требуют оценки первой и / или второй производной функции правдоподобия.

Доказательство правильности [ править ]

Максимизация ожиданий работает скорее на улучшение , чем на прямое улучшение . Здесь показано, что улучшения первого подразумевают улучшения второго. ^[14]${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$${\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}$

Для любого с ненулевой вероятностью мы можем написать${\displaystyle \mathbf {Z} }$${\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})}$

${\displaystyle \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 ${\displaystyle \mathbf {Z} }$ under the current parameter estimate ${\displaystyle \theta ^{(t)}}$ by multiplying both sides by ${\displaystyle p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})}$ and summing (or integrating) over ${\displaystyle \mathbf {Z} }$. The left-hand side is the expectation of a constant, so we get:

${\displaystyle {\begin{aligned}\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{aligned}}}$

where ${\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$ is defined by the negated sum it is replacing. This last equation holds for every value of ${\displaystyle {\boldsymbol {\theta }}}$ including ${\displaystyle {\boldsymbol {\theta }}={\boldsymbol {\theta }}^{(t)}}$,

${\displaystyle \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

${\displaystyle \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 ${\displaystyle H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\geq H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})}$, so we can conclude that

${\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})-\log p(\mathbf {X} \mid {\boldsymbol {\theta }}^{(t)})\geq Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})-Q({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)}).}$

In words, choosing ${\displaystyle {\boldsymbol {\theta }}}$ to improve ${\displaystyle Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})}$ causes ${\displaystyle \log p(\mathbf {X} \mid {\boldsymbol {\theta }})}$ to improve at least as much.

As a maximization–maximization procedure[edit]

The EM algorithm can be viewed as two alternating maximization steps, that is, as an example of coordinate descent.^[15][16] Consider the function:

${\displaystyle F(q,\theta ):=\operatorname {E} _{q}[\log L(\theta ;x,Z)]+H(q),}$

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

${\displaystyle F(q,\theta )=-D_{\mathrm {KL} }{\big (}q\parallel p_{Z\mid X}(\cdot \mid x;\theta ){\big )}+\log L(\theta ;x),}$

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

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

Expectation step: Choose ${\displaystyle q}$ to maximize ${\displaystyle F}$:

${\displaystyle q^{(t)}=\operatorname {arg\,max} _{q}\ F(q,\theta ^{(t)})}$

Maximization step: Choose ${\displaystyle \theta }$ to maximize ${\displaystyle F}$:

${\displaystyle \theta ^{(t+1)}=\operatorname {arg\,max} _{\theta }\ F(q^{(t)},\theta )}$

Applications[edit]

EM is frequently used for parameter estimation of mixed models,^[17][18] notably in quantitative genetics.^[19]

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.^{[citation needed]}

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)^[20] 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 frequently used for data clustering, computer vision and in machine learning. 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.

Filtering and smoothing EM algorithms[edit]

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

${\displaystyle {\widehat {\sigma }}_{v}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{(z_{k}-{\widehat {x}}_{k})}^{2},}$

where ${\displaystyle {\widehat {x}}_{k}}$ are scalar output estimates calculated by a filter or a smoother from N scalar measurements ${\displaystyle 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

${\displaystyle {\widehat {\sigma }}_{w}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}^{2},}$

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

${\displaystyle {\widehat {F}}={\frac {\sum _{k=1}^{N}({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}{\sum _{k=1}^{N}{\widehat {x}}_{k}^{2}}}.}$

The convergence of parameter estimates such as those above are well studied.^{[21][22][23][24]}

Variants[edit]

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).^[25] 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".^[26]

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.^[27] Itself can be extended into the Expectation conditional maximization either (ECME) algorithm.^[28]

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.^[15] GEM is further developed in a distributed environment and shows promising results.^[29]

It is also possible to consider the EM algorithm as a subclass of the MM (Majorize/Minimize or Minorize/Maximize, depending on context) algorithm,^[30] and therefore use any machinery developed in the more general case.

α-EM algorithm[edit]

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^[31]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.^[32]

Relation to variational Bayes methods[edit]

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, 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[edit]

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[edit]

Gaussian mixture[edit]

Let ${\displaystyle \mathbf {x} =(\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{n})}$ be a sample of ${\displaystyle n}$ independent observations from a mixture of two multivariate normal distributions of dimension ${\displaystyle d}$, and let ${\displaystyle \mathbf {z} =(z_{1},z_{2},\ldots ,z_{n})}$ be the latent variables that determine the component from which the observation originates.^[16]

${\displaystyle X_{i}\mid (Z_{i}=1)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{1},\Sigma _{1})}$ and ${\displaystyle X_{i}\mid (Z_{i}=2)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{2},\Sigma _{2}),}$

where

${\displaystyle \operatorname {P} (Z_{i}=1)=\tau _{1}\,}$ and ${\displaystyle \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:

${\displaystyle \theta ={\big (}{\boldsymbol {\tau }},{\boldsymbol {\mu }}_{1},{\boldsymbol {\mu }}_{2},\Sigma _{1},\Sigma _{2}{\big )},}$

where the incomplete-data likelihood function is

${\displaystyle 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

${\displaystyle 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)},}$

or

${\displaystyle 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\},}$

where ${\displaystyle \mathbb {I} }$ is an indicator function and ${\displaystyle f}$ is the probability density function of a multivariate normal.

In the last equality, for each i, one indicator ${\displaystyle \mathbb {I} (z_{i}=j)}$ is equal to zero, and one indicator is equal to one. The inner sum thus reduces to one term.

E step[edit]

Given our current estimate of the parameters θ^(t), the conditional distribution of the Z_i is determined by Bayes theorem to be the proportional height of the normal density weighted by τ:

${\displaystyle 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)})}}.}$

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:

${\displaystyle {\begin{aligned}Q(\theta \mid \theta ^{(t)})&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,\mathbf {\theta } ^{(t)}}[\log L(\theta ;\mathbf {x} ,\mathbf {Z} )]\\&=\operatorname {E} _{\mathbf {Z} \mid \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 {\theta } ^{(t)}}[\sum _{i=1}^{n}\log L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\sum _{i=1}^{n}\operatorname {E} _{Z_{i}\mid \mathbf {X} ;\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{aligned}}}$

The expectation of ${\displaystyle \log L(\theta ;\mathbf {x} _{i},Z_{i})}$ inside the sum is taken with respect to the probability density function ${\displaystyle P(Z_{i}\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})}$, which might be different for each ${\displaystyle \mathbf {x} _{i}}$ of the training set. Everything in the E step is known before the step is taken except ${\displaystyle 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.

M step[edit]

Q(θ | θ^(t)) being quadratic in form means that determining the maximizing values of θ is relatively straightforward. Also, τ, (μ₁,Σ₁) and (μ₂,Σ₂) may all be maximized independently since they all appear in separate linear terms.

To begin, consider τ, which has the constraint τ₁ + τ₂=1:

${\displaystyle {\begin{aligned}{\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{aligned}}}$

This has the same form as the MLE for the binomial distribution, so

${\displaystyle \tau _{j}^{(t+1)}={\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)}.}$

For the next estimates of (μ₁,Σ₁):

${\displaystyle {\begin{aligned}({\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{aligned}}.}$

This has the same form as a weighted MLE for a normal distribution, so

${\displaystyle {\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)}}}}$ and ${\displaystyle \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)}}}}$

and, by symmetry,

${\displaystyle {\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)}}}}$ and ${\displaystyle \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)}}}.}$

Termination[edit]

Conclude the iterative process if ${\displaystyle 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 }$ for ${\displaystyle \varepsilon }$ below some preset threshold.

Generalization[edit]

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

Truncated and censored regression[edit]

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

Alternatives[edit]

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^[34] or the so-called spectral techniques^{[35][36][citation needed]}. Moment-based approaches to learning the parameters of a probabilistic model are of increasing interest recently since they enjoy guarantees such as global convergence under certain conditions unlike EM which is often plagued by the issue of getting stuck in local optima. Algorithms with guarantees for learning can be derived for a number of important models such as mixture models, HMMs etc. For these spectral methods, no spurious local optima occur, and the true parameters can be consistently estimated under some regularity conditions^{[citation needed]}.

See also[edit]

• mixture distribution
• compound distribution
• density estimation
• total absorption spectroscopy
• The EM algorithm can be viewed as a special case of the majorize-minimization (MM) algorithm.^[37]

References[edit]

Further reading[edit]

• Hogg, Robert; McKean, Joseph; Craig, Allen (2005). Introduction to Mathematical Statistics. Upper Saddle River, NJ: Pearson Prentice Hall. pp. 359–364.
• Dellaert, Frank (2002). "The Expectation Maximization Algorithm". CiteSeerX 10.1.1.9.9735. Cite journal requires |journal= (help) gives an easier explanation of EM algorithm as to lowerbound maximization.
• Bishop, Christopher M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2.CS1 maint: ref duplicates default (link)
• Gupta, M. R.; Chen, Y. (2010). "Theory and Use of the EM Algorithm". Foundations and Trends in Signal Processing. 4 (3): 223–296. CiteSeerX 10.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.
• Bilmes, Jeff (1998). "A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models". CiteSeerX 10.1.1.28.613. Cite journal requires |journal= (help) includes a simplified derivation of the EM equations for Gaussian Mixtures and Gaussian Mixture Hidden Markov Models.
• McLachlan, Geoffrey J.; Krishnan, Thriyambakam (2008). The EM Algorithm and Extensions (2nd ed.). Hoboken: Wiley. ISBN 978-0-471-20170-0.

External links[edit]

• 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.
• k-MLE: A fast algorithm for learning statistical mixture models
• 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 k-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.