.nr PS 10
.nr VS 12
.\"fspecial TR CTR
.\"fspecial TI CTI
.\"fspecial TB CTB
.\"kern 1 
.special SS S SA
.fspecial SS SA
.ds ABSTRACT РЕФЕРАТ
.TL 
Практическая реализация EM-алгоритма для смеси гауссовских распределений
.AU
Г.А. Ситкарев, <sitkarev@unixkomi.ru>
.ND Декабрь, 2011
.AI
Сыктывкарский Государственный Университет
.EQ
delim ##
define vdots '\[u22EE]'
define ddots '\[u22F1]'
define rdtos '\[u22F0]'
.EN
.NH 1
Условия задачи
.PP
Дано #N# точек
#( bold x sub 1 , bold x sub 2 , ... , bold x sub N )#
размерности #M#.
Известно,
что все точки принадлежат #K#
многомерным гауссовским распределениям вида
#omega sub k^ eta ( bold x ~|~ bold mu sub k,~SIGMA sub k )#
с неизвестными параметрами
#bold mu sub k#,
#SIGMA sub k#,
#omega sub k#:
.RS
.IP #bold mu sub k#
вектор средних значений для #k#-го распределения;
.IP #SIGMA sub k#
ковариационная матрица #M times M# для #k#-го распределения;
.IP #omega sub k#
вес #k#-го распределения, #sum from k omega sub k = 1#.
.RE
.PP
Каждая точка #bold x sub n# имеет некоторую вероятность
присутствия в распределении #k#.
Эту вероятность мы будем обозначать как #P sub {nk}#.
Все параметры #P sub {nk}# для всех точек удобно свести
в матрицу #P sub {N times K}#:
.EQ
P sub {N times K} ~=~
 left [
  pile {
   P sub {11} above
   P sub {21} above
   vdots above 
   P sub {N1}
  } ~
  pile {
   P sub {12} above
   P sub {22} above
   vdots above
   P sub {N2}
  } ~
  pile {
   P sub {13} above
   P sub {23} above
   vdots above
   P sub {N3}
  } ~ 
  pile {
   ... above
   ... above
   ddots above
   ...
  } ~
  pile {
   P sub {1K} above
   P sub {2K} above
   vdots above
   P sub {NK}
  }
 right ] .
.EN
.NH 1
Правдоподобность и вероятности #P sub {nk} #
.PP
Вероятностную плотность для точки с координатами 
#bold x# будем обозначать как\ #P( bold x )#.
Это вероятность попасть в точку по этим координатам,
если случайным образом выбирать её из всей выборки.
Правдоподобность параметров распределения будем обозначать
как \[Lscript].
Наша задача максимизировать \[Lscript],
подобрав параметры #bold mu sub k#, #SIGMA sub k#,
#omega sub k#.
Так как мы полагаем точки с данными независимыми,
их \[Lscript] будет произведением вероятностей присутствия точки из
выборки в #bold x sub n#:
.EQ
 \[Lscript]~=~ prod from n=1 to N P( bold x sub n ) =
  P( bold x sub 1 ) ^times^
  P( bold x sub 2 ) ^times^ ... ^times^
  P( bold x sub N ).
.EN
Для того,
чтобы посчитать \[Lscript],
нам нужно найти плотности распределения в точках
#bold x sub n# и перемножить их.
Для каждой точки #bold x sub n# плотность определяется
как сумма плотностей из #K# распределений,
пропорционально весу #k#-го распределения:
.EQ
P( bold x sub n ) = sum from {k=1} to K
  omega sub k ^ p sub j ^ ( bold x).
.EN
Так как значения #P sub {nk}#,
то есть значения строки #n# матрицы #P sub {N times K}#,
есть вероятности вхождения #bold x sub n# в
распределение #k#,
их сумма должна равняться единице:
.EQ
sum from k=1 to k P sub {nk} = 1.
.EN
Значит #P( bold x sub n )# это сумма
.EQ
P( bold x sub n ) =
 P sub {n1} cdot P( bold x sub n ) ~+~
 P sub {n2} cdot P( bold x sub n ) ~+~ ... ~+~
 P sub {NK} cdot P( bold x sub n ),
.EN
но нам также известно,
что
.EQ
P( bold x sub n ) =
 omega sub 1 ^ eta ( bold x sub n ~|~ bold mu sub 1 ,~SIGMA sub 1 )
 ~+~
 omega sub 2 ^ eta ( bold x sub n ~|~ bold mu sub 2 ,~SIGMA sub 2 )
 ~+~ ... ~+~
 omega sub K ^ eta ( bold x sub n ~|~ bold mu sub K ,~SIGMA sub K ).
.EN
Последнее означает,
что каждое значение #P sub {nk}# можно вычислить по формуле:
.EQ
P sub {nk} =
 {omega sub k ^ eta ( bold x sub n ~|~ bold mu sub k ,~SIGMA sub k )} over 
 {P( bold x sub n )}.
.EN
При заполнении матрицы #P sub {N times K}#
значениями #P sub {nk}#,
вычисления рационально построить по следующей схеме:
.RS
Для каждой строки #n# матрицы #P sub {N times K}#:
.IP а)
Заполнить строку значениями
#P sub {nk} cdot P( bold x sub n ) =
omega sub k ^ eta ( bold x sub n ~|~ bold mu sub k ,~SIGMA sub k ).#
.IP б)
Подсчитать сумму всех элементов в строке #sum P sub {nk}#.
.IP в)
Поделить все элементы строки на сумму,
полученную на предыдущем шаге.
.RE
.NH 1
Формулы максимизации
.PP
Значения #bold mu bar sub k ,~SIGMA bar sub k ,~omega bar sub k#
на шаге максимизации вычисляются по следующим формулам:
.DS
.EQ (1)
omega bar sub k ~=~ 1 over N sum from n P sub {nk}
.EN
.EQ (2)
bold mu bar sub k ~=~
 {sum from n P sub {nk} cdot bold x sub n } over {sum from n to nothing P sub {nk} }
.EN
.EQ (3)
SIGMA bar sub k ~=~
 {
  sum from n
   P sub {nk} cdot ( bold x sub n - italic {bold mu} bar sub k )
    ( bold x sub n - bold mu bar sub k ) sup T 
 } over
 {
  sum from n to nothing
   P sub {nk}
 }
.EN
.DE
.NH 1
Пошаговое выполнение алгоритма
.PP
.nr step 1 1
Пользователь предоставляет данные #bold x sub n#,
размерности #N#, #M# и #K#,
значение #epsilon# и,
возможно,
вектор начальных значений для #bold mu sub k#.
Реализация алгоритма может состоять из следующих шагов:
.IP \n[step]. 3
На первом шаге:
.RS
.IP а)
Установить веса #omega sub k#, как #omega sub k = 1 over k#.
.IP б)
Установить #SIGMA sub k# в
#I = size -2 
{
  left [ matrix { 
  ccol { 1 above 0 above size -2 {vdots} above 0 }
  ccol { 0 above 1 above size -2 {vdots} above 0 }
  ccol { size -1 {...} above size -1 {...} above size -2 {ddots} above size -1 {...} }
  ccol { 0 above 0 above size -2 {vdots} above 1 }
 } right ]
}.#
.IP в)
Выбрать #bold mu sub k#,
как #K# случайно взятых значения из набора #bold x sub n#,
или же взять начальные значения,
предоставленные пользователем.
.RE
.IP \n+[step].
На втором шаге:
.RS
.IP а)
Выполнить шаг \fIE\fP:
.RS
.IP а.1)
вычислить #P sub {nk}# и #P( bold x sub n )#;
.IP а.2)
вычислить сумму всех логарифмов #P( bold x sub n )#;
.IP а.3)
сохранить #size -2 {sum from n} log P( bold x sub n )#
в переменную #loglike_prev#.
.RE
.RE
.RS
.IP б)
Выполнить шаг \fIM\fP:
.RS
.IP б.1)
пользуясь формулой (1) и весами #P sub {N times K}#,
вычислить #omega bar sub k#;
.IP б.2)
пользуясь формулой (2) и весами #P sub {N times K}#,
вычислить #bold mu bar sub k#;
.IP б.3)
пользуясь формулой (3) и весами #P sub {N times K}#,
вычислить #SIGMA bar sub k#.
.RE
.RE
.IP \n+[step].
На третьем шаге:
.RS
.IP а)
Выполнить шаг \fIE\fP,
как в 2а,
но сохраняя результат 2а.3 в переменную \fIloglike\fP.
.IP б)
Вычислить #abs(loglike - loglike_prev)#,
и сравнить с #epsilon#,
заданным пользователем.
Если #abs(loglike - loglike_prev) < epsilon#,
тогда перейти на шаг 4.
.IP в)
Выполнить шаг \fIM\fP, как в 2б.
.IP г)
Перейти на шаг 3.
.RE
.IP \n+[step].
Завершить выполнение алгоритма с успешным статусом,
а значения\ #bold mu bar sub k,~SIGMA bar sub k,~omega bar sub k#
вернуть пользователю.
