Восстановление смеси плотностей распределений с помощью EM-алгоритма.

Евгений Борисов

четверг, 26 декабря 2013 г.

В этой статье мы поговорим о задаче восстановления смеси плотностей распределений точек и о методе, который называется EM-алгоритм (expectation–maximization algorithm). Этот метод можно применять для решения разных практических задач, например для решения задачи кластеризации.

1 Введение

Рассмотрим множество точек X, поэтому набору точек требуется восстановить плотность их распределения. Пусть плотность имеет следующий вид

∑ k ∑k p(x) = wjφ (x;𝜃j) ; wj = 1 ; wj ≥ 0 j=1 j=1

где 𝜃_j - параметры, φ - заданная функция, k - количество компонент смеси, w_j - весовой коэффициент (априорная вероятность компоненты)

Задача: найти оптимальные 𝜃_j,w_j по простой выборке X и заданным k и φ.

2 Базовый вариант EM-алгоритма

Для решения описанной выше задачи можно использовать EM-алгоритм, который состоит из двух основных процедур: E-шаг - оценка параметров и М-Шаг - оптимизация параметров. Базовый вариант EM-алгоритма выглядит следующим образом.

EM(X,k,w,𝜃,δ)

выбрать начальные значения w и 𝜃
stepE: вычислить оценку параметров $g′ij := gij$

$---wj-φ(xi,𝜃j)--- gij := ∑k w φ(x ,𝜃 ) s=1 s i s$
i = 1,…,m - номер объекта в выборке,
j = 1,…,k - номер компоненты смеси распределения
stepM: оптимизировать параметры

$m -1 ∑ wj := m gij i=1$
$∑m 𝜃j := arg max gij ln φ (xi;𝜃) 𝜃 i=1$ (1)
если max _i,j|g_ij − g′_ij| > δ
то переход на п.2
где 0 < δ < 1 – максимальное значение изменения параметров (условие остановки EM)
конец

Это базовый вариант метода. Результат его работы зависит от начальных значений параметров, здесь не вполне ясно как выбирать количество компонент k и начальные параметры 𝜃_j. Для решения этих задач существует модификация базового алгоритма - EM-алгоритм с последовательным добавлением компонент.

3 EM-алгоритм с последовательным добавлением компонент

Идея этой модификации ЕМ-алгоритма вполне описывается в названии. Мы последовательно добавляем компоненты к смеси, на каждом шаге выполняем базовый ЕМ и оцениваем качество результата по определённой процедуре. Этот вариант EM-алгоритма выглядит следующим образом.

выбрать начальные значения параметров
k = 1, w = 1
параметры 𝜃 выбираются в зависимости от вида функции φ,
например, если φ это многомерная нормальная (гауссовская) плотность
то 𝜃₁ = (w₁,μ₁, Σ₁), где μ₁ - мат.ожидание, Σ₁ - матрица ковариаций на X.
найти максимальное значение плотности p на X

$pmax = mxa∈xX p (x )$ $∑ k p(x) = wjφ (x;𝜃j) j=1$
выбрать из X точки с плотностью ниже пороговой R ⋅ p_max,
где 0 < R < 1 – коэффициент правдоподобия

$U = {x ∈ X |p(x) < R ⋅ pmax }$
если |U| < m₀
то переход на п.8
где 0 < m₀ < |X| – минимальная длинна выборки для обработки
добавляем компоненту смеси

$k = k + 1$

$mk = |U|$

$w = mk- ; w = w (1 − w ) ; j = 1,...k − 1 k m j j k$
параметры для нормальной (гауссовской) плотности φ

$w ∑ μk = --k xi mk xi∈U$ $Σk = cov(U )$
выполняем базовую процедуру оптимизации

$EM (X, k,w, 𝜃,δ)$
где 0 < δ < 1 - максимальное значение изменения параметров(условие остановки EM)
переход на п.2
конец

4 Реализация

В этом разделе мы рассмотрим реализации базового и последовательного EM-алгоритмов для точек на плоскости (n = 2).

Определим φ(x; 𝜃) как многомерную нормальную (гауссовская) плотность

Рис.1: функция Гаусса

( ) exp--−-12(x-−-μ)TΣ-−1(x-−-μ)-- φ (x;𝜃) = N (x;μ,Σ ) = ∘ ----n------ (2π) detΣ

где μ ∈ ℝⁿ - мат.ожидание (центр) x ∈ X, Σ ∈ ℝ^n×n матрица ковариаций X. Таким образом 𝜃_j это совокупность (μ_j, Σ_j) мат.ожидания и матрицы ковариаций для компоненты смеси j.

В этом случае решение задачи (1) по нахождению оптимальных 𝜃_j = (μ_j, Σ_j) выглядит следующим образом (j = 1,…,k)

( | ∑m |||| Nj = gij ||| i=1 |{ 1 ∑m μj = --- gijxi ||| Nj i=1 ||| ∑m |||( Σj = -1- gij(xi − mj )(xi − mj )T Nj i=1

4.1 Реализация базового EM-алгоритма

В этом примере будем искать смесь из двух компонент (k = 2).

Начальные значения μ₁ - мат.ожидание на X, μ₂ генерируется случайным образом, начальные матрицы ковариаций Σ₁ = Σ₁ = cov(X)

Отметим, что результат зависит от начальных значений параметров, и выбранный способ инициализации является не самым удачным, в данном случае он был выбран из-за его простоты.

На рисунках ниже проиллюстрирована работа ЕМ-алгоритма, который восстанавливает смесь плотностей.


Рис.2: начальное состояние	Рис.3: результат работы ЕМ


Рис.4: изменение параметра g	Рис.5: анимация процесса (кликабельно)

Реализация в системе Octave [ здесь ].

4.2 Реализация последовательного EM-алгоритма

На рисунках ниже проиллюстрирована работа последовательного ЕМ-алгоритма, который восстанавливает смесь плотностей.


Рис.7: начальное состояние	Рис.8: промежуточный результат работы ЕМ (две компоненты)


Рис.9: промежуточный результат работы ЕМ (три компоненты)	Рис.10: окончательный результат работы ЕМ (четыре компоненты)

На рисунках ниже проиллюстрировано изменение количества точек с низким значением плотности (отмечены желтым цветом) на разных этапах выполнения алгоритма.


Рис.11: начальное состояние	Рис.12: промежуточный результат работы ЕМ (две компоненты)


Рис.13: промежуточный результат работы ЕМ (три компоненты)	Рис.14: окончательный результат работы ЕМ (четыре компоненты)

На рисунках ниже проиллюстрировано изменение параметра g на разных этапах выполнения алгоритма.


Рис.15: промежуточный результат работы ЕМ (две компоненты)	Рис.16: промежуточный результат работы ЕМ (три компоненты)


Рис.17: промежуточный результат работы ЕМ (четыре компоненты)

Реализация в системе Octave [ здесь ].

Список литературы

[1] Воронцов К.В. Статистические методы классификации – http://shad.yandex.ru/lectures/machine_learning.xml

[2] GNU Octave – http://www.gnu.org/software/octave/

При использовании материалов этого сайта, пожалуйста вставляйте в свой текст ссылку на мою статью.