EM-алгоритм — различия между версиями

Материал из Викиконспекты
Перейти к: навигация, поиск
(Оформил докозательство теоремы из шага M)
м (rollbackEdits.php mass rollback)
 
(не показана 21 промежуточная версия 4 участников)
Строка 1: Строка 1:
 
== Определение ==
 
== Определение ==
  
'''Алгоритм EM''' (англ: expectation-maximization) --- итеративный алгоритм поиска оценок максимума правдоподобия модели, в ситуации, когда она зависит от скрытых(ненаблюдаемых) переменных.
+
'''Алгоритм EM''' (англ. ''expectation-maximization'') {{---}} итеративный алгоритм поиска оценок максимума правдоподобия модели, в ситуации, когда она зависит от скрытых (ненаблюдаемых) переменных.
  
 
Алгоритм ищет параметры модели итеративно, каждая итерация состоит из двух шагов:
 
Алгоритм ищет параметры модели итеративно, каждая итерация состоит из двух шагов:
  
'''E(Expectation)''' шаг --- поиск наиболее вероятных значений скрытых переменных.
+
'''E (Expectation)''' шаг {{---}} поиск наиболее вероятных значений скрытых переменных.
  
'''M(Maximisation)''' шаг --- поиск наиболее вероятных значений параметров, для полученных на шаге E значений скрытых переменных.
+
'''M (Maximization)''' шаг {{---}} поиск наиболее вероятных значений параметров, для полученных на шаге E значений скрытых переменных.
  
 
EM алгоритм подходит для решения задач двух типов:
 
EM алгоритм подходит для решения задач двух типов:
  
 
# Задачи с неполными данными.
 
# Задачи с неполными данными.
# Задачи, в которых удобно вводить скрытые переменные для упрощения подсчета функции правдоподобия. Примером такой задачи может служить кластеризация
+
# Задачи, в которых удобно вводить скрытые переменные для упрощения подсчета функции правдоподобия. Примером такой задачи может служить кластеризация.
  
 
== Основной алгоритм ==  
 
== Основной алгоритм ==  
Строка 19: Строка 19:
  
 
Плотность распределения смеси имеет вид:<br/>
 
Плотность распределения смеси имеет вид:<br/>
<tex>p(x) = \sum\limits_{j=1}^k \omega_j p_j(x)</tex><br/>
+
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex>.<br/>
Где <tex> \sum\limits_{j=1}^k w_j = 1;  w_j \geq 0; p_j(x) = \phi(x;\theta_j)</tex> - функция правдоподобия <tex>j</tex>-ой компонеты смеси, <tex>\omega_j</tex> - априорная вероятность <tex>j</tex>-ой компоненты распределения.<br/>
+
Где <tex> \sum\limits_{j=1}^k w_j = 1;  w_j \geq 0; p_j(x) = \phi(x;\theta_j)</tex> {{---}} функция правдоподобия <tex>j</tex>-ой компонеты смеси, <tex>w_j</tex> {{---}} априорная вероятность <tex>j</tex>-ой компоненты смеси.<br/>
  
 
Перед нами стоит две задачи:<br/>
 
Перед нами стоит две задачи:<br/>
  
# По заданной выборке <tex>X^m</tex> случайных и независимых наблюдений полученных из смеси <tex>p(x)</tex>, числу <tex>k</tex> и функцию <tex>\phi</tex>, оценить вектор параметров <tex>\Theta = (\omega_1,..,\omega_k,\theta_1,..,\theta_k)</tex>
+
# По заданной выборке <tex>X^m</tex> случайных и независимых наблюдений полученных из смеси <tex>p(x)</tex>, числу <tex>k</tex> и функции <tex>\phi</tex>, оценить вектор параметров <tex>\Theta = (w_1,..,w_k,\theta_1,..,\theta_k)</tex>.
# Найти <tex>k</tex>
+
# Найти <tex>k</tex>.
  
 
=== Проблема ===  
 
=== Проблема ===  
  
Задачи подобного рода мы умеем решать максимизируя логармиф правдоподобия:<br>
+
Задачи подобного рода мы умеем решать, максимизируя логармиф правдоподобия:<br>
<tex> L(\Theta) = ln \prod\limits_{i=1}^mp(x_i) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow max</tex><br/>
+
<tex> Q(\Theta) = ln \prod\limits_{i=1}^mp(x_i) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}</tex>.<br/>
  
Но пробелма в том, что мы не знаем как аналитически посчитать логарифм суммы. Тут нам и поможет алгоритм EM.
+
Но проблeма в том, что мы не знаем как аналитически посчитать логарифм суммы. Тут нам и поможет алгоритм EM.
  
 
=== Решение ===  
 
=== Решение ===  
  
Основная идея алгоритма EM заключается в том, что мы добавляме скрытые переменные такие, что:<br/>
+
Основная идея алгоритма EM заключается в том, что мы добавляем скрытые переменные такие, что:<br/>
  
# Они могут быть выражены через <tex>\Theta</tex>
+
# Они могут быть выражены через <tex>\Theta</tex>.
# Они помогают разбить сумму так: <tex>p (X, H|\Theta) = \prod\limits_{i=1}^k p (X|H, \Theta) p(H|\Theta)</tex>, где <tex>H</tex> - матрица скрытых переменных.
+
# Они помогают разбить сумму так: <tex>p (X, H|\Theta) = \prod\limits_{i=1}^k p (X|H, \Theta) p(H|\Theta)</tex>, где <tex>H</tex> {{---}} матрица скрытых переменных.
  
 
Тогда алгоритм EM сводится к повторению шагов, указанных в [[#Определение|Определении]].
 
Тогда алгоритм EM сводится к повторению шагов, указанных в [[#Определение|Определении]].
Строка 45: Строка 45:
 
=== E-шаг ===  
 
=== E-шаг ===  
  
<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> <br />
+
<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex>.<br />
  
Скрытые переменные представляют из себя матрицу <tex>H = (h_{ij})_{m \times k}</tex><br/>
+
Скрытые переменные представляют из себя матрицу <tex>H = (h_{ij})_{m \times k}</tex>,<br/>
где <tex>h_{ij} = P(\theta_j | x_i)</tex> - вероятность того, что <tex>x_i</tex> пренадлежит <tex>j</tex>-ой компоненте.<br/>
+
где <tex>h_{ij} = P(\theta_j | x_i)</tex> {{---}} вероятность того, что <tex>x_i</tex> пренадлежит <tex>j</tex>-ой компоненте.<br/>
  
 
По формуле Байеса справедливо равенство:<br />
 
По формуле Байеса справедливо равенство:<br />
<tex> h_{ij} = \frac{w_jp_j(x_i)}{p (x_i)} = \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^k w_s p_s(x_i)}</tex><br/>
+
<tex> h_{ij} = \frac{w_jp_j(x_i)}{p (x_i)} = \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^k w_s p_s(x_i)}</tex>.<br/>
  
Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex><br/>
+
Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex>.<br/>
  
Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных<br/>
+
Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных.<br/>
  
 
=== M-шаг ===  
 
=== M-шаг ===  
Строка 62: Строка 62:
 
|statement=
 
|statement=
 
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам:<br/>
 
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам:<br/>
<center><tex>\theta_j = \arg\max\limits{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
+
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex>.</center>
 
Оптимальные же веса считаются как:<br/>
 
Оптимальные же веса считаются как:<br/>
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
+
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>.</center>
 
|proof=
 
|proof=
 
Посчитаем логарифм правдоподобия:<br>
 
Посчитаем логарифм правдоподобия:<br>
<tex> Q(\Theta) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow max</tex><br/>
+
<tex> Q(\Theta) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}</tex>.<br/>
 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи:<br/>
 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи:<br/>
<tex> L(\Theta, X^m) = \sum\limits_{i=1}^m ln \biggl( \sum\limits_{j=1}^k w_j p_j(x_i) \biggr)  - \lambda \biggl(\sum\limits_{j=1}^k w_j - 1 \biggr) </tex><br/>
+
<tex> L(\Theta, X^m) = \sum\limits_{i=1}^m ln \biggl( \sum\limits_{j=1}^k w_j p_j(x_i) \biggr)  - \lambda \biggl(\sum\limits_{j=1}^k w_j - 1 \biggr) </tex>.<br/>
 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим:<br/>
 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим:<br/>
<tex>\frac{\partial L} {\partial w_j} = \sum\limits_{i=1}^m \frac{p_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} - \lambda = 0, j = 1..k</tex><br />
+
<tex>\frac{\partial L} {\partial w_j} = \sum\limits_{i=1}^m \frac{p_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} - \lambda = 0, j = 1..k</tex>.<br/>
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> <br />
+
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex>:<br />
<tex>\sum\limits_{j=1}^k \sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} = \lambda \sum\limits_{j=1}^kw_j</tex> <br />
+
<tex>\sum\limits_{j=1}^k \sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} = \lambda \sum\limits_{j=1}^kw_j</tex>.<br />
А так как <tex>\sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = 1</tex> и <tex>\sum\limits_{j=1}^kw_j = 1</tex>, из чего следует <tex>\lambda = m</tex> <br />
+
А так как <tex>\sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = 1</tex> и <tex>\sum\limits_{j=1}^kw_j = 1</tex>, из чего следует <tex>\lambda = m</tex>.<br />
  
<tex>w_i = \frac{1}{m}\sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = \frac{1}{m}\sum\limits_{i=1}^m h_{ij}</tex><br />
+
<tex>w_j = \frac{1}{m}\sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = \frac{1}{m}\sum\limits_{i=1}^m h_{ij}</tex>.<br />
  
Приравняв к нулю лагранжиан по <tex>\theta_j</tex> схожим способом найдем:<br />
+
Приравняв к нулю производную Лагранжиана по <tex>\theta_j</tex>, схожим способом найдем:<br />
  
<tex> \theta_j = \arg\max\limits{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</tex><br />
+
<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</tex><br />
  
  
Строка 87: Строка 87:
 
=== Критерий остановки ===  
 
=== Критерий остановки ===  
  
Алгоритм EM вы полняется до сходимости, но как нам определить, что сходимость наступила?. Мы можем останавливаться, когда либо <tex>Q(\Theta)</tex>, либо <tex>H</tex> перестают сильно меняться. Но, обычно, удобней контролировать изменения значений скрытых переменных, так как они имеют смысл вероятностей и принимают значения из отрезка <tex>[0,1]</tex>. Поэтому один из возможных критериев остановки будет выглядеть так: <tex>max_{i,j} |h_{ij} - h_{ij}^{(0)}| > \delta</tex>
+
Алгоритм EM выполняется до сходимости, но как нам определить, что сходимость наступила? Мы можем останавливаться, когда либо <tex>Q(\Theta)</tex>, либо <tex>H</tex> перестают сильно меняться. Но, обычно, удобней контролировать изменения значений скрытых переменных, так как они имеют смысл вероятностей и принимают значения из отрезка <tex>[0,1]</tex>. Поэтому один из возможных критериев остановки будет выглядеть так: <tex>\max\limits_{i,j} |h_{ij} - h_{ij}^{(0)}| > \delta</tex>.
  
 
=== Псевдокод ===  
 
=== Псевдокод ===  
  
  Input:<tex>X^m, k, \theta^{(0)}</tex>
+
  Input:<tex>X^m, k, \Theta^{(0)}</tex>
 
  Repeat
 
  Repeat
 
     '''E-step''': for all i = 1..m; j = 1..k:
 
     '''E-step''': for all i = 1..m; j = 1..k:
         <tex>h_{ij} = \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^k w_s p_s(x_i)}</tex>
+
         <tex>h_{ij} = \frac{w_j \phi(x_i; \theta_j)}{\sum\limits_{s=1}^k w_s \phi(x_i; \theta_j)}</tex>
 
     '''M-step''': for all j = 1..k:
 
     '''M-step''': for all j = 1..k:
         <tex>\theta_j = argmax_{\theta} \sum\limits_{i=1}^m h_{ij}*ln \phi (x_i, \theta)</tex>
+
         <tex>\theta_j = \arg\max\limits_{\theta} \sum\limits_{i=1}^m h_{ij}*ln \phi (x_i, \theta)</tex>
 
         <tex>w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>
 
         <tex>w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>
  Until stopping criteria is satisfied
+
  Until a stopping criterion is satisfied
 +
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>
  
 
=== Плюсы и минусы ===  
 
=== Плюсы и минусы ===  
Строка 104: Строка 105:
 
Плюсы:<br/>
 
Плюсы:<br/>
  
* Сходится в большинтсве ситуаций
+
* Сходится в большинтсве случаев.
* Наиболее гибкое решение
+
* Наиболее гибкое решение.
* Легко может добавить нечуствительность к шуму
+
* Существуют простые модификации, позволяющие уменьшить чуствительность алгоритма к шуму в данных.
  
 
Минусы:<br/>
 
Минусы:<br/>
  
* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму
+
* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
* Число компонент <tex>k</tex> является гиперпараметром.
+
* Число компонент <tex>k</tex> является [[Настройка_гиперпараметров|гиперпараметром]].
  
 
== Модификации ==  
 
== Модификации ==  
  
Базовый алгоритм EM является очень гибким для модификаций, позволяющих улучшить его работу. В этом разделе мы приведем краткое описаниен некоторых из них.
+
Базовый алгоритм EM является очень гибким для модификаций, позволяющих улучшить его работу. В этом разделе мы приведем краткое описание некоторых из них.
  
 
=== Generalized EM-algorithm ===  
 
=== Generalized EM-algorithm ===  
Строка 123: Строка 124:
 
=== Stochastic EM-algorithm ===
 
=== Stochastic EM-algorithm ===
  
Как уже было отмечено в [[#Плюсы_и_минусы|Плюсы и минусы]], базовый алгоритм чувствителен к начальному приближению и могут быть ситуации, когда алгоритм "застрянет" в лоакльном экстремуме. Для того, чтобы предотвратить это, будем на каждой итерации алгоритма случайно "встряхивать" выборку. В этой модификации у нас добавляется шаг S, на котором мы и будем "встряхивать" выборку. И на шаге M мы будем решать уже задачу максимуму невзвешенного правдоподобия. Эта модификация хороша тем, что нечуствиетльная к начальном приблежению.  
+
Как уже было отмечено в [[#Плюсы_и_минусы|Плюсы и минусы]], базовый алгоритм чувствителен к начальному приближению и могут быть ситуации, когда алгоритм "застрянет" в локальном экстремуме. Для того, чтобы предотвратить это, будем на каждой итерации алгоритма случайно "встряхивать" выборку. В этой модификации у нас добавляется шаг S, на котором мы и будем "встряхивать" выборку. И на шаге M мы будем решать уже задачу максимуму невзвешенного правдоподобия. Эта модификация хороша тем, что нечуствиетльная к начальном приблежению.  
  
 
== Пример. Разделение смеси Гауссиана ==  
 
== Пример. Разделение смеси Гауссиана ==  
Строка 131: Строка 132:
 
Каноническим примером использования EM алгоритма является задача разделения смеси гауссиана. Данные у нас получены из нормального распределения. В этом случае параметрами функций ялвяются матожидание и дисперсия.<br/>
 
Каноническим примером использования EM алгоритма является задача разделения смеси гауссиана. Данные у нас получены из нормального распределения. В этом случае параметрами функций ялвяются матожидание и дисперсия.<br/>
  
<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> вектор параметров, <br/>
+
<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, <br/>
<tex>p_j(x) = N(x;\mu_j, \sigma_j) = \frac1{\sqrt{2\pi}\sigma_j} \exp \biggl(-\frac{(x - \mu_j)^2}{2\sigma_j^2}\biggr) </tex> - плотность распределения <br/>
+
<tex>p_j(x) = N(x;\mu_j, \sigma_j) = \frac1{\sqrt{2\pi}\sigma_j} \exp \biggl(-\frac{(x - \mu_j)^2}{2\sigma_j^2}\biggr) </tex> {{---}} плотность распределения.<br/>
  
 
Посчитаем значения для каждого шага. <br/>
 
Посчитаем значения для каждого шага. <br/>
Строка 138: Строка 139:
 
E-шаг:
 
E-шаг:
  
: <tex> h_{ij} = \frac{w_j N(x_i, \mu_j, \sigma_j)}{\sum\limits_{s=1}^k w_s N(x_i, \mu_s, \sigma_s)} </tex>
+
: <tex> h_{ij} = \frac{w_j N(x_i, \mu_j, \sigma_j)}{\sum\limits_{s=1}^k w_s N(x_i, \mu_s, \sigma_s)}.</tex>
  
 
M-шаг:
 
M-шаг:
  
: <tex>w_j = \frac{1}{m} \sum\limits_{i=1}^m h_{ij}</tex>
+
: <tex>w_j = \frac{1}{m} \sum\limits_{i=1}^m h_{ij}.</tex>
: <tex> \mu_j = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}x_i</tex>
+
: <tex> \mu_j = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}x_i.</tex>
: <tex> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>
+
: <tex> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.</tex>
  
 
== Использование в задаче кластеризации ==
 
== Использование в задаче кластеризации ==
Строка 150: Строка 151:
 
[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]
 
[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]
  
Как уже упоминалось в [[#Опредиление|Определении]], алгоритм EM подходит для решения задачи кластеризации. И одной из его имплементаций для это задачи является алгоритм <tex>k</tex>-means. В этом алгоритме в качестве скрытых переменных выступают метки классов объектов. Параметрами же являются центроиды искомых классов. Тогда на шаге E мы относим объекты к какому то одному классу на основе расстояний до центроид. А на шаге M мы пересчитываем центроиды кластеров, исходя из полученной на шаге E разметке.<br/>
+
Как уже упоминалось в [[#Определение|Определении]], алгоритм EM подходит для решения задачи кластеризации. И одной из его имплементаций для этой задачи является алгоритм [[Алгоритм_k-Means|<tex>k</tex>-Means]]. В этом алгоритме в качестве скрытых переменных выступают метки классов объектов. Параметрами же являются центроиды искомых классов. Тогда на шаге E мы относим объекты к какому-то одному классу на основе расстояний до центроид. А на шаге M мы пересчитываем центроиды кластеров, исходя из полученной на шаге E разметке.<br/>
  
Также стоит упомянуть алгоритм <tex>c</tex>-means. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.
+
Также стоит упомянуть алгоритм <tex>c</tex>-means<ref>[https://en.wikipedia.org/wiki/Fuzzy_clustering#Fuzzy_C-means_clustering C-means clustering, Wikipedia]</ref>. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.
  
  
Строка 160: Строка 161:
 
== Реализация на python ==
 
== Реализация на python ==
  
 +
В пакете sklearn алгоритм EM представлен объектом GaussianMixture. Проиллюстрируем его работу на примере задачи кластеризации и сравним его с алгоритмом <tex>k</tex>-means:
 +
 +
[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
 
  '''import''' numpy as np
 
  '''import''' numpy as np
 
  '''import''' matplotlib.pyplot as plt
 
  '''import''' matplotlib.pyplot as plt
Строка 167: Строка 171:
 
  np.random.seed(12)
 
  np.random.seed(12)
 
   
 
   
  # Создаем datasets с использованием стандартных sklearn.datasets
+
  <font color="green"># Создаем datasets с использованием стандартных sklearn.datasets</font>
 
  n_samples = 2000
 
  n_samples = 2000
 
  random_state = 170
 
  random_state = 170
Строка 175: Строка 179:
 
  varied = datasets.make_blobs(n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state)
 
  varied = datasets.make_blobs(n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state)
 
   
 
   
  # Создаем анизатропно разделенные данные
+
  <font color="green"># Создаем анизатропно разделенные данные</font>
 
  X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
 
  X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
 
  transformation = [[0.6, -0.6], [-0.4, 0.8]]
 
  transformation = [[0.6, -0.6], [-0.4, 0.8]]
Строка 181: Строка 185:
 
  aniso = (X_aniso, y)
 
  aniso = (X_aniso, y)
 
   
 
   
  # Выставляем параметры для matplotlib.pyplot
+
  <font color="green"># Выставляем параметры для matplotlib.pyplot</font>
 
  plt.figure(figsize=(9 * 2 + 3, 12.5))
 
  plt.figure(figsize=(9 * 2 + 3, 12.5))
 
  plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
 
  plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
Строка 187: Строка 191:
 
  defaul_n = 3
 
  defaul_n = 3
 
   
 
   
  # Варьируем значение количества классов в зависимости от данных, ведь для нас это гиперпараметр
+
  <font color="green"># Варьируем значение количества классов в зависимости от данных, ведь для нас это гиперпараметр</font>
 
  datasets = [
 
  datasets = [
 
     (varied, defaul_n),
 
     (varied, defaul_n),
Строка 196: Строка 200:
 
     X, y = dataset
 
     X, y = dataset
 
   
 
   
     # Нормализация данных
+
     <font color="green"># Нормализация данных</font>
 
     X = StandardScaler().fit_transform(X)
 
     X = StandardScaler().fit_transform(X)
 
   
 
   
     # Непосредственно наш алгоритм - Gaussian Mixture
+
     <font color="green"># Непосредственно наш алгоритм - Gaussian Mixture</font>
 
     gmm = mixture.GaussianMixture(n_components=n_cluster, covariance_type='full')
 
     gmm = mixture.GaussianMixture(n_components=n_cluster, covariance_type='full')
 
      
 
      
     # Для сравнения берем алгоритм - K-means
+
     <font color="green"># Для сравнения берем алгоритм - k-means</font>
 
     two_means = cluster.KMeans(n_clusters=n_cluster)
 
     two_means = cluster.KMeans(n_clusters=n_cluster)
     clustering_algorithms = (
+
     clustering_algorithms = {
         ('GaussianMixture', gmm),
+
         'Real distribution': None,
         ('KMeans', two_means)
+
        'Gaussian Mixture': gmm,
     )
+
         'k-Means': two_means
 +
     }
 
     for name, algorithm in clustering_algorithms:
 
     for name, algorithm in clustering_algorithms:
 
   
 
   
 
         # Этап обучения
 
         # Этап обучения
         algorithm.fit(X)
+
         if algorithm is not None:
 +
            algorithm.fit(X)
 
          
 
          
 
         # Применяем алгоритм
 
         # Применяем алгоритм
        y_pred = algorithm.predict(X)
+
        y_pred = y if algorithm is None else algorithm.predict(X)
 
          
 
          
 
         # Рисуем результаты
 
         # Рисуем результаты
Строка 229: Строка 235:
 
  plt.show()
 
  plt.show()
  
[[Файл:Prog.png|thumb|250px|Результат программы]]
+
Как и следовало ожидать, алгоритм EM работает на некоторых данных лучше чем k-means, однако есть данные, с которыми он не справляется без дополнительных преобразований.
 
 
Как и следовало ожидать, алгоритм работает на некоторых данных лучше чем k-means, однако есть данные, с которыми он не справляется без дополнительных преобразований.
 
  
 
== См. также ==
 
== См. также ==
 
*[[Кластеризация]]
 
*[[Кластеризация]]
 +
*[[Алгоритм_k-Means|Алгоритм k-Means]]
 +
 +
==Примечания==
 +
<references />
  
 
== Источники информации ==
 
== Источники информации ==
  
 +
# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
 
# [http://www.machinelearning.ru/wiki/images/6/6d/Voron-ML-1.pdf Математические методы обучения по прецедентам К. В. Воронцов]
 
# [http://www.machinelearning.ru/wiki/images/6/6d/Voron-ML-1.pdf Математические методы обучения по прецедентам К. В. Воронцов]
 +
# [http://www.machinelearning.ru/wiki/index.php?title=EM-%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC Статья про EM-алгоритм на machinelearning.ru]
 +
# [https://machinelearningmastery.com/expectation-maximization-em-algorithm/ A Gentle Introduction to Expectation-Maximization]
 
# [http://dendroid.sk/2011/05/09/k-means-clustering/ k-means]
 
# [http://dendroid.sk/2011/05/09/k-means-clustering/ k-means]
  
[[Категория:Машинное обучение]]
+
[[Категория: Машинное обучение]]
 +
[[Категория: Кластеризация]]

Текущая версия на 19:36, 4 сентября 2022

Определение

Алгоритм EM (англ. expectation-maximization) — итеративный алгоритм поиска оценок максимума правдоподобия модели, в ситуации, когда она зависит от скрытых (ненаблюдаемых) переменных.

Алгоритм ищет параметры модели итеративно, каждая итерация состоит из двух шагов:

E (Expectation) шаг — поиск наиболее вероятных значений скрытых переменных.

M (Maximization) шаг — поиск наиболее вероятных значений параметров, для полученных на шаге E значений скрытых переменных.

EM алгоритм подходит для решения задач двух типов:

  1. Задачи с неполными данными.
  2. Задачи, в которых удобно вводить скрытые переменные для упрощения подсчета функции правдоподобия. Примером такой задачи может служить кластеризация.

Основной алгоритм

Постановка задачи

Плотность распределения смеси имеет вид:
[math]p(x) = \sum\limits_{j=1}^k w_j p_j(x)[/math].
Где [math] \sum\limits_{j=1}^k w_j = 1; w_j \geq 0; p_j(x) = \phi(x;\theta_j)[/math] — функция правдоподобия [math]j[/math]-ой компонеты смеси, [math]w_j[/math] — априорная вероятность [math]j[/math]-ой компоненты смеси.

Перед нами стоит две задачи:

  1. По заданной выборке [math]X^m[/math] случайных и независимых наблюдений полученных из смеси [math]p(x)[/math], числу [math]k[/math] и функции [math]\phi[/math], оценить вектор параметров [math]\Theta = (w_1,..,w_k,\theta_1,..,\theta_k)[/math].
  2. Найти [math]k[/math].

Проблема

Задачи подобного рода мы умеем решать, максимизируя логармиф правдоподобия:
[math] Q(\Theta) = ln \prod\limits_{i=1}^mp(x_i) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}[/math].

Но проблeма в том, что мы не знаем как аналитически посчитать логарифм суммы. Тут нам и поможет алгоритм EM.

Решение

Основная идея алгоритма EM заключается в том, что мы добавляем скрытые переменные такие, что:

  1. Они могут быть выражены через [math]\Theta[/math].
  2. Они помогают разбить сумму так: [math]p (X, H|\Theta) = \prod\limits_{i=1}^k p (X|H, \Theta) p(H|\Theta)[/math], где [math]H[/math] — матрица скрытых переменных.

Тогда алгоритм EM сводится к повторению шагов, указанных в Определении.

E-шаг

[math]p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)[/math].

Скрытые переменные представляют из себя матрицу [math]H = (h_{ij})_{m \times k}[/math],
где [math]h_{ij} = P(\theta_j | x_i)[/math] — вероятность того, что [math]x_i[/math] пренадлежит [math]j[/math]-ой компоненте.

По формуле Байеса справедливо равенство:
[math] h_{ij} = \frac{w_jp_j(x_i)}{p (x_i)} = \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^k w_s p_s(x_i)}[/math].

Также [math]\sum\limits_{j=1}^k h_{ij} = 1[/math].

Таким образом, зная значения вектора параметров [math]\Theta[/math], мы легко можем пересчитать значения скрытых переменных.

M-шаг

Теорема:
Если известны скрытые переменные, то задача минимизации [math]Q(\Theta)[/math] сводится к [math]k[/math] независимым подзадачам:
[math]\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)[/math].

Оптимальные же веса считаются как:

[math] w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}[/math].
Доказательство:
[math]\triangleright[/math]

Посчитаем логарифм правдоподобия:
[math] Q(\Theta) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}[/math].
При условии, что[math] \sum\limits_{j=1}^k w_j = 1; w_j \geq 0[/math] имеет смысл рассматривать Лагранжиан задачи:
[math] L(\Theta, X^m) = \sum\limits_{i=1}^m ln \biggl( \sum\limits_{j=1}^k w_j p_j(x_i) \biggr) - \lambda \biggl(\sum\limits_{j=1}^k w_j - 1 \biggr) [/math].
Приравняв нулю производную Лагранжиана по [math]w_j[/math], получим:
[math]\frac{\partial L} {\partial w_j} = \sum\limits_{i=1}^m \frac{p_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} - \lambda = 0, j = 1..k[/math].
Умножим на [math]w_j[/math] и просуммируем уравнения для всех [math]j[/math]:
[math]\sum\limits_{j=1}^k \sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} = \lambda \sum\limits_{j=1}^kw_j[/math].
А так как [math]\sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = 1[/math] и [math]\sum\limits_{j=1}^kw_j = 1[/math], из чего следует [math]\lambda = m[/math].

[math]w_j = \frac{1}{m}\sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = \frac{1}{m}\sum\limits_{i=1}^m h_{ij}[/math].

Приравняв к нулю производную Лагранжиана по [math]\theta_j[/math], схожим способом найдем:

[math] \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).[/math]
[math]\triangleleft[/math]

Критерий остановки

Алгоритм EM выполняется до сходимости, но как нам определить, что сходимость наступила? Мы можем останавливаться, когда либо [math]Q(\Theta)[/math], либо [math]H[/math] перестают сильно меняться. Но, обычно, удобней контролировать изменения значений скрытых переменных, так как они имеют смысл вероятностей и принимают значения из отрезка [math][0,1][/math]. Поэтому один из возможных критериев остановки будет выглядеть так: [math]\max\limits_{i,j} |h_{ij} - h_{ij}^{(0)}| \gt \delta[/math].

Псевдокод

Input:[math]X^m, k, \Theta^{(0)}[/math]
Repeat
   E-step: for all i = 1..m; j = 1..k:
       [math]h_{ij} = \frac{w_j \phi(x_i; \theta_j)}{\sum\limits_{s=1}^k w_s \phi(x_i; \theta_j)}[/math]
   M-step: for all j = 1..k:
       [math]\theta_j = \arg\max\limits_{\theta} \sum\limits_{i=1}^m h_{ij}*ln \phi (x_i, \theta)[/math]
       [math]w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}[/math]
Until a stopping criterion is satisfied
Return [math]\Theta = (\theta_j, w_j)_{j=1}^k[/math]

Плюсы и минусы

Плюсы:

  • Сходится в большинтсве случаев.
  • Наиболее гибкое решение.
  • Существуют простые модификации, позволяющие уменьшить чуствительность алгоритма к шуму в данных.

Минусы:

  • Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
  • Число компонент [math]k[/math] является гиперпараметром.

Модификации

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

Generalized EM-algorithm

Осоновная идея этой модификации заключается в том, что на шаге M мы не будем пытаться найти наилучшее решение. Это применимо в случаях, когда максимизация [math]Q(\Theta)[/math] является сликшом дорогой, поэтому нам достаточно сделать лишь несколько итераций, для того, чтобы сместиться в сторону максимума значения [math]Q(\Theta)[/math]. Эта модификация имеет неплохую сходимость.

Stochastic EM-algorithm

Как уже было отмечено в Плюсы и минусы, базовый алгоритм чувствителен к начальному приближению и могут быть ситуации, когда алгоритм "застрянет" в локальном экстремуме. Для того, чтобы предотвратить это, будем на каждой итерации алгоритма случайно "встряхивать" выборку. В этой модификации у нас добавляется шаг S, на котором мы и будем "встряхивать" выборку. И на шаге M мы будем решать уже задачу максимуму невзвешенного правдоподобия. Эта модификация хороша тем, что нечуствиетльная к начальном приблежению.

Пример. Разделение смеси Гауссиана

Несколько итераций алгоритма

Каноническим примером использования EM алгоритма является задача разделения смеси гауссиана. Данные у нас получены из нормального распределения. В этом случае параметрами функций ялвяются матожидание и дисперсия.

[math]\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)[/math] — вектор параметров,
[math]p_j(x) = N(x;\mu_j, \sigma_j) = \frac1{\sqrt{2\pi}\sigma_j} \exp \biggl(-\frac{(x - \mu_j)^2}{2\sigma_j^2}\biggr) [/math] — плотность распределения.

Посчитаем значения для каждого шага.

E-шаг:

[math] h_{ij} = \frac{w_j N(x_i, \mu_j, \sigma_j)}{\sum\limits_{s=1}^k w_s N(x_i, \mu_s, \sigma_s)}.[/math]

M-шаг:

[math]w_j = \frac{1}{m} \sum\limits_{i=1}^m h_{ij}.[/math]
[math] \mu_j = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}x_i.[/math]
[math] \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.[/math]

Использование в задаче кластеризации

Пример работы k-means

Как уже упоминалось в Определении, алгоритм EM подходит для решения задачи кластеризации. И одной из его имплементаций для этой задачи является алгоритм [math]k[/math]-Means. В этом алгоритме в качестве скрытых переменных выступают метки классов объектов. Параметрами же являются центроиды искомых классов. Тогда на шаге E мы относим объекты к какому-то одному классу на основе расстояний до центроид. А на шаге M мы пересчитываем центроиды кластеров, исходя из полученной на шаге E разметке.

Также стоит упомянуть алгоритм [math]c[/math]-means[1]. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.



Реализация на python

В пакете sklearn алгоритм EM представлен объектом GaussianMixture. Проиллюстрируем его работу на примере задачи кластеризации и сравним его с алгоритмом [math]k[/math]-means:

Результат выполнения программы
import numpy as np
import matplotlib.pyplot as plt
from sklearn import cluster, datasets, mixture
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice
np.random.seed(12)

# Создаем datasets с использованием стандартных sklearn.datasets
n_samples = 2000
random_state = 170
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5, noise=.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)
varied = datasets.make_blobs(n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state)

# Создаем анизатропно разделенные данные
X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
aniso = (X_aniso, y)

# Выставляем параметры для matplotlib.pyplot
plt.figure(figsize=(9 * 2 + 3, 12.5))
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
plot_num = 1
defaul_n = 3

# Варьируем значение количества классов в зависимости от данных, ведь для нас это гиперпараметр
datasets = [
   (varied, defaul_n),
   (aniso, defaul_n),
   (blobs, defaul_n),
   (noisy_circles, 2)]
for i_dataset, (dataset, n_cluster) in enumerate(datasets):
   X, y = dataset

   # Нормализация данных
   X = StandardScaler().fit_transform(X)

   # Непосредственно наш алгоритм - Gaussian Mixture
   gmm = mixture.GaussianMixture(n_components=n_cluster, covariance_type='full')
   
   # Для сравнения берем алгоритм - k-means
   two_means = cluster.KMeans(n_clusters=n_cluster)
   clustering_algorithms = {
       'Real distribution': None,
       'Gaussian Mixture': gmm,
       'k-Means': two_means
   }
   for name, algorithm in clustering_algorithms:

       # Этап обучения
       if algorithm is not None:
           algorithm.fit(X)
       
       # Применяем алгоритм
        y_pred = y if algorithm is None else algorithm.predict(X)
       
       # Рисуем результаты
       plt.subplot(len(datasets), len(clustering_algorithms), plot_num)
       if i_dataset == 0:
           plt.title(name, size=18)
       colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a']), int(max(y_pred) + 1))))
       plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])
       plt.xlim(-2.5, 2.5)
       plt.ylim(-2.5, 2.5)
       plt.xticks(())
       plt.yticks(())
       plot_num += 1
plt.show()

Как и следовало ожидать, алгоритм EM работает на некоторых данных лучше чем k-means, однако есть данные, с которыми он не справляется без дополнительных преобразований.

См. также

Примечания

Источники информации

  1. Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
  2. Математические методы обучения по прецедентам К. В. Воронцов
  3. Статья про EM-алгоритм на machinelearning.ru
  4. A Gentle Introduction to Expectation-Maximization
  5. k-means