EM-алгоритм

EM-алгоритм (англ. expectation-maximization) — алгоритм, используемый в математической статистике для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных. Каждая итерация алгоритма состоит из двух шагов. На E-шаге (expectation) вычисляется ожидаемое значение функции правдоподобия, при этом скрытые переменные рассматриваются как наблюдаемые. На M-шаге (maximization) вычисляется оценка максимального правдоподобия, таким образом увеличивается ожидаемое правдоподобие, вычисляемое на E-шаге. Затем это значение используется для E-шага на следующей итерации. Алгоритм выполняется до сходимости.

Идея, суть, смысл

Положим, что модель распределения наблюдаемых объектов $X$ описывается не только с помощью некоторых параметров $\theta$, но и с использованием скрытых переменных $Z$. В таком случае логарифм правдоподобия наблюдаемой выборки записывается в виде $L(\theta) = \log\left(\mathbb{P}(X \mid \theta)\right) = \log\left(\sum_{Z} \mathbb{P}(X, Z \mid \theta)\right)$. В методе максимального правдоподобия этот логарифм максимизируется по параметру $\theta$. Однако максимизировать логарифм от суммы некруто: получаются аналитически нерешаемые уравнения, локальные оптимумы и куча прочих неприятностей... Идея!

Давайте придумаем способ для каждого фиксированного $\hat{\theta}$ находить оптимальную оценку снизу для $L(\theta)$, т. е. $g(\theta): g(\hat{\theta}) = L(\hat{\theta}), g(\theta) \le L(\theta), \forall \theta$, максимизировать ее по параметру, а затем проделывать все заново до сходимости!

Ок, пускай у нас есть некоторое распределение на скрытых переменных: $q(Z): q(Z) \ge 0, \sum_{Z} q(Z) = 1$. Тогда $L(\theta )=\log \left(\mathbb {P} (X\mid \theta )\right)=\log \left(\sum _{Z}\mathbb {P} (X,Z\mid \theta )\right)=\log \left(\sum _{Z}q(Z){\frac {\mathbb {P} (X,Z\mid \theta )}{q(Z)}}\right)\geq \sum _{Z}q(Z)\log \left({\frac {\mathbb {P} (X,Z\mid \theta )}{q(Z)}}\right)=g(\theta )$. Неравенство здесь откуда? А, это ж Йенсен для строго вогнутого логарифма от случайной величины $U$, которая принимает значения $\frac{\mathbb{P}(X, Z \mid \theta)}{q(Z)})$ с вероятностями $q(Z)$. Вот мы получили нижнюю оценку для $L(\theta)$.

Давайте теперь эту оценку сделаем оптимальной для некоторого фиксированного $\hat{\theta}$. Вспомним, что Йенсен обращается в равенство для строго вогнутых функций тогда и только тогда, когда случайная величина, относительно которой записывается неравенство, почти всюду совпадает со своим матожиданием. Ага, ну тогда: $\frac{\mathbb{P}(X, Z \mid \hat{\theta})}{q(Z)} = \mathbb{E}U = \sum_{Z} \mathbb{P}(X, Z \mid \hat{\theta}) = \mathbb{P}(X \mid \hat{\theta}), \forall Z$. Из чего немедленно следует, что $q(Z) = \frac{\mathbb{P}(X, Z \mid \hat{\theta})}{\mathbb{P}(X \mid \hat{\theta})} = \mathbb{P}(Z \mid X,\hat{\theta})$.

Отлично! Теперь эту оптимальную оценку мы, как и обещали, будем максимизировать по параметру: ${\hat {\theta }}_{new}={\underset {\theta }{\arg \max }}g(\theta )={\underset {\theta }{\arg \max }}\sum _{Z}q(Z)\log \left({\frac {\mathbb {P} (X,Z\mid \theta )}{q(Z)}}\right)={\underset {\theta }{\arg \max }}\sum _{Z}q(Z)\log \left(\mathbb {P} (X,Z\mid \theta )\right)={\underset {\theta }{\arg \max }}\mathbb {E} _{Z\sim q(Z)=\mathbb {P} (Z\mid X,{\hat {\theta }}_{old})}\log \left(\mathbb {P} (X,Z\mid \theta )\right)$. Мы там выкинули часть не зависящую от параметра, разбив логарифм частного на разность логарифмов.

Все, вот, собственно, и весь алгоритм. Часть с настройкой распределения на скрытых переменных при фиксированном параметре называется E-шагом, а часть с оптимизацией по параметру нижней оценки при фиксированном распределении на скрытых переменных — M-шагом. Сам алгоритм приведен тут.

Свойства, комментарии, замечания

Ну и что, норм алгоритм? А то! Оказывается, получаемые оценки параметра $\hat{\theta_{0}}, \hat{\theta_{1}}, \hat{\theta_{2}}, \dots$ дают неубывающую последовательность значений логарифма правдоподобия $L(\hat{\theta_{0}}) \le L(\hat{\theta_{1}}) \le L(\hat{\theta_{2}}) \le \dots$. Мда? Да! Ведь для любых соседних $\hat{\theta}_{old}$ и $\hat{\theta}_{new}$ выполняются следующие соотношения: $L(\hat{\theta}_{old}) = g_{old}(\hat{\theta}_{old}) \le g_{old}(\hat{\theta}_{new}) \le L(\hat{\theta}_{new})$. Это справедливо в силу того, что оценка $g_{old}(\theta)$ точна при $\theta = \hat{\theta}_{old}$ и является обычной оценкой снизу при прочих значениях параметра, а также в силу M-шага.

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

Что можно сказать по этому поводу?

• Имеем в виду, что распределение на скрытых переменных не зависит от параметра на M-шаге, так как зафиксировано на E-шаге.
• Итерировать можно до сходимости по аргументу, по правдоподобию, по числу итераций.
• Так как EM сходится к локальному оптимому, было бы неплохо запустить его из нескольких начальных приближений, а затем взять лучший результат (по величине логарифма правдоподобия) в качестве окончательной оценки параметра.
• Есть еще такая штука, как Generalized EM. Тут идея такова, что для сходимости на M-шаге не обязательно искать прям точь-в-точь аргмаксимум. По сути достаточно даже грубого одношагового приближения. Ну лан, будем знать.
• EM, как вы, мальчики и девочки, наверное, заметили, применяется для ML-оптимизации. Но его можно применять и для нахождения MAP-оценок (об этом далее).
• Если хотите, рассматривайте EM как покоординатный подъем, мне-то что.

Случай независимых наблюдений

А что, если $(x_n, z_n) \sim i. i. d.$? Ну, например, если скрытая переменная отвечают за компоненту смеси, из которой происходит объект, или если скрытая переменная — недостающие переменные в независимых и одинаково распределенных наблюдаемых объектах. О, ну тогда совсем другое дело:

• E-шаг:  $q(Z)=mathbb{P}(Z\mid X,\theta )=p(z_{1}\mid x_{1},\theta )\dots p(z_{n}\mid x_{n},\theta )=q_{1}(z_{1})\dots q_{N}(z_{N}),q_{i}(z_{i})=p(z_{i}\mid x_{i},\theta )$.
• M-шаг:  ${\hat {\theta }}_{new}={\underset {\theta }{\arg \max }}\sum _{Z}q(Z)\log \left(\mathbb {P} (X,Z\mid \theta )\right)={\underset {\theta }{\arg \max }}\sum _{Z}q(Z)sum_{i=1}^{N}\log \left(p(x_{i},z_{i}\mid \theta )\right)={\underset {\theta }{\arg \max }}sum_{i=1}^{N}\sum _{z_{1},\dots ,z_{n}}q(z_{1})\dots q(z_{n})\log \left(p(x_{i},z_{i}\mid \theta )\right)={\underset {\theta }{\arg \max }}sum_{i=1}^{N}sum_{z_{i}}q(z_{i})\log \left(p(x_{i},z_{i}\mid \theta )\right)$.

Регуляризация (тут же про MAP)

Ок, давайте добавим регуляризатор к нашему логарифму правдоподобия: $L(\theta) = \log\left(\mathbb{P}(X \mid \theta)\right) + \lambda R(\theta)$. Давайте не будем переписывать все заново, а методом пристального взгляда осознаем, что нижняя оценка станет такой: $g(\theta) = \sum_{Z} q(Z) \log\left(\frac{\mathbb{P}(X, Z \mid \theta)}{q(Z)}\right) + \lambda R(\theta)$, поэтому E-шаг вообще никак не изменится: мы полагаем оптимальным распределением на скрытых переменных все то же $q(Z) = \mathbb{P}(Z \mid X,\hat{\theta})$ при некотором фиксированном $\hat{\theta}$ (действительно, при фиксированном параметре регуляризационная добавка есть просто добавление константы, поэтому оптимум никуда не сдвинется). На M-шаге мы совершим все те же действия, на каждом этапе таща за собой регуляризационную добавку, поэтому конечный результат: $\hat{\theta}_{new} = \underset{\theta}{\arg\max}\left\{ \mathbb{E}_{Z \sim q(Z) = \mathbb{P}(Z \mid X,\hat{\theta}_{old})} \log\left(\mathbb{P}(X, Z \mid \theta)\right) + \lambda R(\theta)\right\}$. В итоге измененный алгоритм получится таким.

Хорошо, но что насчет MAP-оценки? А в чем суть оценки Maximum A Posteriori? Мы полагаем (Байесовский подход), что параметр $\theta$ — это не что-то строго фиксированное, а тоже в некотором роде случайная величина. Поэтому мы прямо так и говорим: пусть $\theta$ — случайная величина с априорным распределением $p(\theta)$. Тогда по аналогии с методом максимального правдоподобия мы будем максимизировать следующую функцию правдоподобия: $L(\theta) = \log\left(\mathbb{P}(X, \theta)\right) = \log\left(\mathbb{P}(X \mid \theta)\right) + \log\left(p(\theta)\right) \rightarrow \max_{\theta}$. Если приглядеться, то это только что рассмотренный ML с регуляризацией $\lambda R(\theta) = \log\left(p(\theta)\right)$. А посему все (касаемо EM) проделывается в полной аналогии. А вот и результат.

Вывод для гауссиан

Для смеси нормальных распределений полный вывод можно найти по этой ссылке. Затруднения возникают в связи с необходимостью дифференцирования определителей и прочих логарифмов от них, которые появляются в лагранжиане для M-шага; Китов расправляется с этим безобразием при помощи лемм. В остальном всё следует определениям в лоб и требует просто аккуратного расписывания.

P. S.

Сорри за корявые формулы, ох уж эта педивикия...