Викиконспекты - Вклад участника [ru]

EM-алгоритм

2020-03-21T11:51:23Z

94.229.111.244: Исправил положение точки

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex>. 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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>. 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex>. 

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

По формуле Байеса справедливо равенство: 
<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>. 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex>. 

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

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex>.</center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>.</center>
|proof=
Посчитаем логарифм правдоподобия: 
<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>. 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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>. 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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>. 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex>: 
<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>. 
А так как <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>. 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
* Число компонент <tex>k</tex> является [[Настройка_гиперпараметров|гиперпараметром]].

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> {{---}} плотность распределения. 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

Также стоит упомянуть алгоритм <tex>c</tex>-means<ref>[https://en.wikipedia.org/wiki/Fuzzy_clustering#Fuzzy_C-means_clustering C-means clustering, Wikipedia]</ref>. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

==Примечания==
<references />

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория: Машинное обучение]]
[[Категория: Кластеризация]]

EM-алгоритм

2020-03-21T11:49:23Z

94.229.111.244: Поправил примечание

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex>. 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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>. 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex>. 

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

По формуле Байеса справедливо равенство: 
<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>. 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex>. 

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

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex>.</center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>.</center>
|proof=
Посчитаем логарифм правдоподобия: 
<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>. 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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>. 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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>. 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex>: 
<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>. 
А так как <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>. 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
* Число компонент <tex>k</tex> является [[Настройка_гиперпараметров|гиперпараметром]].

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> {{---}} плотность распределения. 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

Также стоит упомянуть алгоритм <tex>c</tex>-means<ref>[https://en.wikipedia.org/wiki/Fuzzy_clustering#Fuzzy_C-means_clustering C-means clustering, Wikipedia]</ref>. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

==Примечания==
<references />

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория: Машинное обучение]]
[[Категория: Кластеризация]]

EM-алгоритм

2020-03-21T11:44:51Z

94.229.111.244: Добавил примечаения

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex>. 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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>. 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex>. 

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

По формуле Байеса справедливо равенство: 
<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>. 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex>. 

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

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex>.</center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>.</center>
|proof=
Посчитаем логарифм правдоподобия: 
<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>. 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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>. 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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>. 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex>: 
<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>. 
А так как <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>. 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
* Число компонент <tex>k</tex> является [[Настройка_гиперпараметров|гиперпараметром]].

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> {{---}} плотность распределения. 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

Также стоит упомянуть алгоритм <tex>c</tex>-means<ref>[https://en.wikipedia.org/wiki/Fuzzy_clustering#Fuzzy_C-means_clustering Fuzzy C-means clustering, Wikipedia]</ref>. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

==Примечания==
<references />

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория: Машинное обучение]]
[[Категория: Кластеризация]]

EM-алгоритм

2020-03-21T11:37:58Z

94.229.111.244: Поправил баги

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex>. 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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>. 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex>. 

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

По формуле Байеса справедливо равенство: 
<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>. 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex>. 

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

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex>.</center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>.</center>
|proof=
Посчитаем логарифм правдоподобия: 
<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>. 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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>. 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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>. 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex>: 
<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>. 
А так как <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>. 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

* Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
* Число компонент <tex>k</tex> является [[Настройка_гиперпараметров|гиперпараметром]].

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> {{---}} плотность распределения. 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

Также стоит упомянуть алгоритм [https://ru.wikipedia.org/wiki/Метод_нечёткой_кластеризации_C-средних <tex>c</tex>-means]. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

==Примечания==
<references />

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория: Машинное обучение]]
[[Категория: Кластеризация]]

EM-алгоритм

2020-03-21T11:05:28Z

94.229.111.244: Исправил тире

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex> 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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\limits_{\Theta}</tex> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
|proof=
Посчитаем логарифм правдоподобия: 
<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> 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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> 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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> 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> 
<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> 
А так как <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> 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

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

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория:Машинное обучение]]

EM-алгоритм

2020-03-21T11:03:54Z

94.229.111.244: Поправил опечатки

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex> 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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\limits_{\Theta}</tex> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
|proof=
Посчитаем логарифм правдоподобия: 
<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> 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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> 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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> 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> 
<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> 
А так как <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> 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

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

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория:Машинное обучение]]

EM-алгоритм

2020-03-21T11:00:55Z

94.229.111.244: Добавил точки

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex> 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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\limits_{\Theta}</tex> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
|proof=
Посчитаем логарифм правдоподобия: 
<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> 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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> 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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> 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> 
<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> 
А так как <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> 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

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

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория:Машинное обучение]]

EM-алгоритм

2020-03-20T20:45:47Z

94.229.111.244: Поправил нерабочую ссылку на Определение

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex> 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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\limits_{\Theta}</tex> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
|proof=
Посчитаем логарифм правдоподобия: 
<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> 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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> 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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> 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> 
<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> 
А так как <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> 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

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

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

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

[[Файл:em_clustering.png|thumb|600px|Результат выполнения программы]]
'''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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория:Машинное обучение]]

EM-алгоритм

2020-03-19T19:10:53Z

94.229.111.244: Поправил кучу мелких багов

== Определение ==

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

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

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

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

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

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

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

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

Плотность распределения смеси имеет вид: 
<tex>p(x) = \sum\limits_{j=1}^k w_j p_j(x)</tex> 
Где <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>-ой компоненты смеси. 

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

# По заданной выборке <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> 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\limits_{\Theta}</tex> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
|proof=
Посчитаем логарифм правдоподобия: 
<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> 
При условии, что<tex> \sum\limits_{j=1}^k w_j = 1; w_j \geq 0</tex> имеет смысл рассматривать Лагранжиан задачи: 
<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> 
Приравняв нулю производную Лагранжиана по <tex>w_j</tex>, получим: 
<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> 
Умножим на <tex>w_j</tex> и просуммируем уравнения для всех <tex>j</tex> 
<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> 
А так как <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> 

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

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

<tex> \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).</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>
Repeat
'''E-step''': for all i = 1..m; j = 1..k:
<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:
<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>
Until a stopping criterion is satisfied
Return <tex>\Theta = (\theta_j, w_j)_{j=1}^k</tex>

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

Плюсы: 

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

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

[[Файл:Gaussians2.png|right|thumb|400px|Несколько итераций алгоритма]]

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> {{---}} вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

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

[[Файл:kmeans.jpg|right|thumb|200px|Пример работы k-means]]

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

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

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

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

[[Файл:Prog.png|thumb|400px|Результат выполнения программы]]
'''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 = (
('GaussianMixture', gmm),
('KMeans', two_means)
)
for name, algorithm in clustering_algorithms:

# Этап обучения
algorithm.fit(X)

# Применяем алгоритм
y_pred = 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, однако есть данные, с которыми он не справляется без дополнительных преобразований.

== См. также ==
*[[Кластеризация]]
*[[Алгоритм_k-Means|Алгоритм k-Means]]

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

# Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
# [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]

[[Категория:Машинное обучение]]

EM-алгоритм

2020-03-17T09:00:58Z

94.229.111.244:

== Определение ==

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

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

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

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

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

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

== Проблема восстановления распределения смеси ==

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

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

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

# По заданной выборке <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>k</tex>

=== Проблема ===

Задачи подобного рода мы умеем решать максимизируя логармиф правдоподобия: 
<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> 

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

=== Решение ===

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

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

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

=== E-шаг ===

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 

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

По формуле Байеса справедливо равенство: 
<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> 

Также <tex>\sum\limits_{j=1}^k h_{ij} = 1</tex> 

Таким образом, зная значения вектора параметров <tex>\Theta</tex>, мы легко можем пересчитать значения скрытых переменных 

=== M-шаг ===

{{Теорема
|statement=
Если известны скрытые переменные, то задача минимизации <tex>Q(\Theta)</tex> сводится к <tex>k</tex> независимым подзадачам: 
<center><tex>\theta_j = argmax_{\theta} \sum\limits_{i=1}^m h_{ij}*ln \phi (x_i, \theta)</tex></center>
Оптимальные же веса считаются как: 
<center><tex> w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex></center>
}}

=== Псевдокод ===

Input:<tex>X^m, k, \theta^{(0)}</tex>
Repeat
'''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>
'''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>w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}</tex>
Until stopping criteria is satisfied

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

Плюсы: 

* Сходится в большинтсве ситуаций
* Наиболее гибкое решение
* Легко может добавить нечуствительность к шуму

Минусы: 

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

== Модификации ==

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

=== Generalized EM-algorithm ===

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

=== Stochastic EM-algorithm ===

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

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

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

<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> — вектор параметров, 
<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> - плотность распределения 

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

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>

M-шаг:

: <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> \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k</tex>

== Задача разделения смеси распределений ==

=== Общий алгоритм ===

Необходимо описать плотность распределения функции на X как сумму k функций, которые можно рассматривать как элементы параметрического семейства функций <tex> p_j(x) = \phi(x;\theta_j)</tex>. Плотность распределения будет выглядеть как 
<tex>p(x) = \sum\limits_{i=1}^k \omega_j p_j(x); \sum\limits_{i=1}^k w_j = 1; w_j >= 0 </tex>
 где <tex>\omega_j</tex>- априорная вероятность j компоненты распределения.
Задача разделения смеси заключается в том, чтобы, имея выборку <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>

E-шаг:

<tex>p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)</tex> 
Введем обозначение: <tex> g_{ij} = P(\theta_j | x_i) </tex> это и будут скрытые параметры данной задачи - апостериорная вероятность того, что обучающий объект <tex> x_i </tex> получен из <tex>j</tex>-й компоненты

По формуле Байеса справедливо равенство: 
<tex> g_{ij} = \frac{w_jp_j(x_i)}{\sum\limits_{t=1}^k w_t p_t(x_i)}</tex> 
Таким образом при зная значение параметров легко найти скрытые переменные.

Перейдем к M-шагу.

Посчитаем для аддитивности логарифм правдоподобия: 
<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) \longrightarrow max</tex> 
при условии <tex>\sum\limits_{i=1}^k w_j = 1; w_j >= 0</tex> имеет смысл рассматривать лагранжиан задачи: 
<tex>\frac{\partial L} {\partial w_j} = \sum\limits_{i=1}^m \frac{p_j(x_i)}{\sum\limits_{t=1}^kw_tp_t(x_i)} - \lambda = 0.</tex> 

Умножим на <tex>\omega_j</tex> и просуммируем уравнения для всех <tex>j</tex> 

<tex>\sum\limits_{j=1}^k \sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{t=1}^kw_tp_t(x_i)} = \lambda \sum\limits_{j=1}^kw_j</tex> 

Так как можно заменить порядок суммы: <tex> \sum\limits_{i=1}^m \sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{t=1}^kw_tp_t(x_i)} = \lambda \sum\limits_{j=1}^kw_j</tex>.

А так как <tex>\sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{t=1}^kw_tp_t(x_i)} = 1</tex> и <tex>\sum\limits_{j=1}^kw_j = 1</tex>, из чего следует <tex>\lambda = m</tex> 

<tex>\omega_j = \frac{1}{m}\sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{t=1}^kw_tp_t(x_i)} = \frac{1}{m}\sum\limits_{i=1}^mg_{ij}</tex> 

Приравняв к нулю лагранжиан по <tex>\theta_j</tex> схожим способом найдем: 

<tex> \theta_j = \arg\max\limits{\theta}\sum\limits_{i=1}^mg_{ij}\ln(\phi(x_i;\theta)).</tex> 

Таким образом на M-шаге необходимо взять среднее значение <tex>g_{ij}</tex> и решить k независимых оптимизационных задач.

=== Разделение смеси гауссиан ===
[[Файл:Gaussians.png|right|250px| Несколько итераций алгоритма]]
Важным на практике примером является случай, когда параметрическое семейство - нормальные распределения. Параметрами функций будут являться матожидание и дисперсия. 
<tex>\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)</tex> — вектор параметров, 
<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>

== k-means как EM алгоритм ==
[[Файл:kmeans.jpg|right|250px|K-means]]
Скрытыми переменными в данной задаче являются классы, к которым относятся объекты для кластеризации. Сами же параметры это центры масс классов. На шаге E - распределяются все объекты по классам исходя из расстояния от центра, на шаге M находится оптимальное месторасположение центра.

Аналогично рассматривается и алгоритм c-means. Скрытые переменные здесь будут вероятности принадлежности к классам, которые находятся на E-шаге по расстоянию от центра. Центр так же рассчитывается на M-шаге исходя из скрытых переменных.

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

'''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 = (
('GaussianMixture', gmm),
('KMeans', two_means)
)
for name, algorithm in clustering_algorithms:

# Этап обучения
algorithm.fit(X)

# Применяем алгоритм
y_pred = 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()

[[Файл:Prog.png|thumb|250px|Результат программы]]

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

== См. также ==
*[[Кластеризация]]

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

# [http://www.machinelearning.ru/wiki/images/6/6d/Voron-ML-1.pdf Математические методы обучения по прецедентам К. В. Воронцов]
# [http://dendroid.sk/2011/05/09/k-means-clustering/ k-means]

[[Категория:Машинное обучение]]