Алгоритам очекивање-максимизација (ЕМ алгоритам)

Из Википедије, слободне енциклопедије

У статистици, алгоритам Очекивање-максимизација (EM - енгл. Expectation-maximization) је итеративна процедура за естимацију параметара на основу критеријума максималне веродостојности (МЛ - енгл. Maximum Likelihood) или оцене апостериорног максимума (МАП - енгл. Maximum a posteriori), код којих су присутне недоступне латентне варијабле (нису видљиве у опсервацијама), које имају особину θ коју не можемо да измеримо, бар не директно.

ЕМ алгоритам се састоји из Е корака, који креира функцију за очекивање од лог веродостојности, израчунату коришћењем тренутне естимације за параметаре, и М корака, који израчунава параметре и максмизује очекивану лог веродостојност, пронађену у Е кораку. Ти процењени пареметри се потом користе за одређивање расподеле латентних варијабли у следећем Е кораку.

Историјат[уреди]

Име ЕМ алгоритма и начин функционисања, дато је у раду из 1977. године, који су написали Артур Демпстер, Нан Лаирд и Доналд Рубин. Њихов рад је генерализовао методу и скицирао анализу конвергенције за ширу класу проблема. Без обзира на раније проналаске, њихова иновативна метода се прославила у јавности и њивов рад окатегорисан као брилијантан. Рад Демпстер-Лаирд-Рубин је основао ЕМ метод као веома важан део статистичке анализе.

Увод[уреди]

ЕМ алгоритам се користи за проналажење параметра максималне веродостојности статистичког модела у случајевима код којих се једначине не могу решити директно. Углавном, ови модели укључују латентне варијабле непознатих параметара и познатих посматраних података. Значи да, или недостају вредности подацима, или се модел може формулисати једноставније претпостављајући да постоје додатне непосматране тачке података.

Проналажење решења максималне веродостојности захтева узимање деривата вероватноће функције у односу на све непознате вредности. Код статистичих модела са латентним варијаблама, ово обично није могуће. Уместо тога, резултат је скуп међусобно повезаних једначина у којој решење за параметре захтева вредности латентних варијабли и обратно, али замена једног скупа једначина у други производи нерешиве једначине. ЕМ алгоритам полази од следећег запажања за решавање ова два сета једначина нумерички. Један једноставно може изабрати произвољне вредности за један од два сета непознатих, користити их да процени други сет, а затим помоћу ових нових вредности нађе бољу процену првоме сету, а онда настави наизменучно између њих, све док резултујуће вредности не конвергирају ка фиксним тачкама. Није очигледно да ће то радити уопште, али у ствари се може доказати у овом контретном контексту да то ради, а да је дериват веродостојности (арбитрарно близу) нула у том тренутку, што значи да је тачка или максимум или седло тачке. У принципу може бити више максимума и никаква гаранција да ће се глобални максимум наћи. Неке веродостојности имају сингуларитете у њима, односно непостојећи максимум.

Опис алгоритма[уреди]

Дати статистички модел који се састоји од скупа \mathbf{X} посматраних података, скуп непримећених латентних података или вредности које недостају \mathbf{Z} и вектр непознатих параметара \boldsymbol\theta, заједно са функцијом веродостојности L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}|\boldsymbol\theta), процена максималне веродостојности (МЛЕ - енгл. Maximum likelihood estimate) од непознатих параметара одређује маргиналне веродостојности посматраних података

L(\boldsymbol\theta; \mathbf{X}) = p(\mathbf{X}|\boldsymbol\theta) = \sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z}|\boldsymbol\theta)

Међутим, ово је често нерешиво (нпр. Ако се деси да број вредности расте експоненцијално што је секвенца дужа, онда ће тачан прорачун суме бити изузетно тежак).

ЕМ алгоритам покушава да пронађе МЛЕ од граничне веродостојности итеративном применом следећа два корака:

Корак очекивања (Е корак): Израчунава очекивану вредност лог веродостојности функције вероватноће, у погледу на условну расподелу \mathbf{Z} датим \mathbf{X} под тренутном проценом параметара \boldsymbol\theta^{(t)}:
Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) = \operatorname{E}_{\mathbf{Z}|\mathbf{X},\boldsymbol\theta^{(t)}}\left[ \log L (\boldsymbol\theta;\mathbf{X},\mathbf{Z})  \right] \,
Корак максимизације (М корак): Проналази параметар који максимизује следеће:
\boldsymbol\theta^{(t+1)} = \underset{\boldsymbol\theta}{\operatorname{arg\,max}} \ Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) \,

Типични модели на које се примењује ЕМ:

  1. Посматране тачке података \mathbf{X} могу бити дискретне (узимајући вредности из коначног или пребројиво бесконачног скупа) или непрекидне (узмањући вредности из непребројиво бесконачног скупа). Могу, у ставри, бити и вектор опсевације повезан са сваком тачком података.
  2. Вредности које недостају (латентне варијабле) \mathbf{Z} су дискретне, извучене из фиксног броја вредности, а постоји једна латентна променљива по посматраној тачки података.
  3. Параметри су непрекидни, и има две врсте: Параметри који су повезани са свим тачкама података, и параметара повезаним са одређеном вредношћу латентне варијабле.

Међутим, могуће је да се ЕМ примени и на друге врсте модела. Ако знамо вредности параметара \boldsymbol\theta, обично се може наћи вредност латентних варијабли \mathbf{Z} повећавањем лог-веродостојности по свим могућим вредностима \mathbf{Z}, или једноставно итеративно преко \mathbf{Z} или преко алгоритма, као што је Витерби алгоритам за скривене Маркове моделе. Насупрот томе, ако знамо вредности латентних варијабли \mathbf{Z}, можемо наћи процену параметрара \boldsymbol\theta прилично лако, једноставним груписањем посматране тачке података на основу вредности придружене латенте варијабле и просека вредности, или нека функција вредности, од тачака у свакој групи. Ово сугерише итеративни алгоритам, у случају када су \boldsymbol\theta и \mathbf{Z}непознати:

  1. Прво, иницијализујте параметре \boldsymbol\theta неким случајним вредностима.
  2. Израчунај најбољу вредност за \mathbf{Z} помоћу ових вредности параметара.
  3. Затим, користите ове израчунате вредности \mathbf{Z} да израчунате бољу процену параметара \boldsymbol\theta. Параметри повезани са одређеном вредношћу \mathbf{Z} ће користити само оне тачке података код којих придружене латентне варијабле имају ту вредност.
  4. Вршите итерацију корака 2 и 3 до конвергенције.

Управо описан алгоритам монотоно прилази локаном минимуму функције, а најчешће се назива хард ЕМ. K-mean алгоритам је пример ове класе алгоритма.

Својства[уреди]

Говорећи о Е кораку, мало је погрешан. Оно што се рачуна у првом кораку су фиксни параметри зависни од података функције Q. Када су параметри Q познати, потпуно су одређени и увећани у другом М кораку ЕМ алогритма. 

Иако ЕМ итерација повећава број посматраних података (тј. маргинали) функције веродостојности, не постоји гаранција да низ конвергира ка процени максималне веродостојности (МЛЕ). За бимодалне дистрибуције, ово значи да ЕМ алгоритам може конвергирати до локалног максимума посматраних података функције веродостојности, зависећи од почетних вредности. Постоји низ хеуристичких или метахеуристичких приступа за "бежање" локалним максимумима, као што су насумични рестарт (кренувши од неколико различитих случајних почетних процена θ(t)), или применом алгоритма симулације жарења.

ЕМ је нарочито користан када је веродостојност у породици експоненцијалних алгоритама: Е корак постаје збир очекивања довољне статистике, а М корак подразумева максимизовање линеарне функције. У том случају обично је могуће извести исправке у затвореној форми за сваки корак, користећи Сундберг формулу (Објавио Ралф Сундберг користећи необјављене резултате Пер Мартин–Лофа и Андерс Мартин-Лофа).

ЕМ метода је модификована за израчунавање максималне апостериорне процене (МАП), са Бајесовом статистиком, у оригиналном раду Демпстер, Лаирд и Рубин.

Постоје и друге методе за проналажење максималне веродостојности, једна од метода је варијација Гаус-Њутновог алгоритма, а постоје и још неке. За разлику од ЕМ, такве методе обично захтевају процену првог или другог деривата функције веродостојности.

Доказ коректности[уреди]

ЕМ алгоритам ради на побољшању Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}), а не побољшава директно \log p(\mathbf{X}|\boldsymbol\theta). Овде смо показали да побољшања првог подразумева побољшање последњг. За било које \mathbf{Z} са не нултом вероватноћом p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta), можемо записати


\log p(\mathbf{X}|\boldsymbol\theta) = \log p(\mathbf{X},\mathbf{Z}|\boldsymbol\theta) - \log p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta) \,.

Узимамо очекивање над вредностима \mathbf{Z} множењем обе стране са p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta^{(t)}). Сабирањем и (или интегрисањем) преко \mathbf{Z}. Лева страна је константа очекивања, па добијамо:


\begin{align}
\log p(\mathbf{X}|\boldsymbol\theta) &
= \sum_{\mathbf{Z}} p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta^{(t)}) \log p(\mathbf{X},\mathbf{Z}|\boldsymbol\theta)
- \sum_{\mathbf{Z}} p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta^{(t)}) \log p(\mathbf{Z}|\mathbf{X},\boldsymbol\theta) \\
& = Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta|\boldsymbol\theta^{(t)}) \,,
\end{align}

где је H(\boldsymbol\theta|\boldsymbol\theta^{(t)}) одређен негирањем суме коју је заменио. Ова последња једначина важи за било коју вредност \boldsymbol\theta, uključujući \boldsymbol\theta = \boldsymbol\theta^{(t)},


\log p(\mathbf{X}|\boldsymbol\theta^{(t)})
= Q(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}) \,,

и одузимањем ове последње једначине са оном из претходне, добијамо


\log p(\mathbf{X}|\boldsymbol\theta) - \log p(\mathbf{X}|\boldsymbol\theta^{(t)})
= Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)})
 + H(\boldsymbol\theta|\boldsymbol\theta^{(t)}) - H(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}) \,,

Међутим, Гисова неједнакост нам говори да H(\boldsymbol\theta|\boldsymbol\theta^{(t)}) \ge H(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}), па можемо закључити да


\log p(\mathbf{X}|\boldsymbol\theta) - \log p(\mathbf{X}|\boldsymbol\theta^{(t)})
\ge Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}) \,.

Другим речима, бирајући \boldsymbol\theta да унапредимо Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) изван Q(\boldsymbol\theta^{(t)}|\boldsymbol\theta^{(t)}) ће унапредити \log p(\mathbf{X}|\boldsymbol\theta) preko \log p(\mathbf{X}|\boldsymbol\theta^{(t)}) најмање толико.

Апликације[уреди]

Под неким околностима, врло је згодно гледати на ЕМ алгоритам као два наизменична корака максимизирања. Размотримо функцију:

F(q,\theta) = \operatorname{E}_q [ \log L (\theta ; x,Z)] + H(q) = -D_{\mathrm{KL}}\big(q \big\| p_{Z|X}(\cdot|x;\theta ) \big) + \log L(\theta;x)

где је q произвољна расподела вероватноћа над непосматраним подацима z, pZ|X(· |x;θ) је условна расподела непосматраних података где су дати посматрани подаци x, H је ентропија и DKL је Кулбак-Лајблер дивергенција.

Онда се кораци ЕМ алгоритма могу посматрати као:

Корак очекивања: Изабери q да максимизујеш F:
 q^{(t)} = \operatorname*{arg\,max}_q \ F(q,\theta^{(t)})
Корак максимизације: Изабери θ да максимизујеш F:
 \theta^{(t+1)} = \operatorname*{arg\,max}_\theta \ F(q^{(t)},\theta)

Пример[уреди]

Гаусова расподела[уреди]

Анимација која демонстрира ЕМ алгоритам користећи модел двокомпонентне Гаусове расподеле над подацима гејзира Стари Верни. Алгоритам иде од насумичне иницијализације до конвергенције.

Нека је x = (x1,x2,…,xn) пример од n независних опсервација из расподеле две мултиваријационе нормалне расподеле димензије d, и нека је z=(z1,z2,…,zn) латентна варијабла којом се одређује компонента из које потиче посматрање.

X_i |(Z_i = 1) \sim \mathcal{N}_d(\boldsymbol{\mu}_1,\Sigma_1) i X_i |(Z_i = 2) \sim \mathcal{N}_d(\boldsymbol{\mu}_2,\Sigma_2)

где је

\operatorname{P} (Z_i = 1 ) = \tau_1 \, i \operatorname{P} (Z_i=2) = \tau_2 = 1-\tau_1

Циљ је да се процене непознати параметари који представљају “мешање” вредности између Гаусових и начини и коваријансе сваког.

\theta = \big(\boldsymbol{\tau},\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma_1,\Sigma_2 \big)

где је функција веродостојности:

L(\theta;\mathbf{x},\mathbf{z}) = P(\mathbf{x},\mathbf{z} \vert \theta) = \prod_{i=1}^n  \sum_{j=1}^2  \mathbb{I}(z_i=j) \ \tau_j \ f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j)

где је \mathbb{I} индикатор функције, а f је расподела вероватноће од више варијанта. Ово може бити преписано фамилији експоненцијалних форми:

L(\theta;\mathbf{x},\mathbf{z}) = \exp \left\{ \sum_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big] \right\}.

Да би се видела последња једнакост, имајте на уму да за свако i сви индикатори \mathbb{I}(z_i=j) су једнаки нули, осим једног који је једнак један. Унутрашња сума се на тај начин смањује на један члан.

Корак Е[уреди]

Тренутна процена параметара θ(t) условна расподела Zi је детерминисана са Бајесовом теоремом, да буде пропорционалне висине од нормалне расподеле вероватноће, са тежином τ:

T_{j,i}^{(t)} := \operatorname{P}(Z_i=j | X_i=\mathbf{x}_i ;\theta^{(t)}) = \frac{\tau_j^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)})}{\tau_1^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_1^{(t)},\Sigma_1^{(t)}) + \tau_2^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_2^{(t)},\Sigma_2^{(t)})} .

Резултат Е корака у функцији:

\begin{align}Q(\theta|\theta^{(t)})
&= \operatorname{E} [\log L(\theta;\mathbf{x},\mathbf{Z})] \\
&= \sum_{i=1}^n \sum_{j=1}^2 T_{j,i}^{(t)} \big[ \log \tau_j  -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big]
\end{align}

Да бисмо видели последњу једнакост, имајте на уму да смо сумирањем по свим могућим вредностима од z где је вероватноћа сваког z производ \prod_{i=1}^n T_{j,i}^{(t)}. Сада погледајмо коефицијенте сваког члана унутар суме за.L(\theta;\mathbf{x},\mathbf{z}). Биће два члана T_{1,i}^{(t)}(\sum_{j=1,2}\prod_{k\ne i}T_{j,k}^{(t)}) and T_{2,i}^{(t)}(\sum_{j=1,2}\prod_{k\ne i}T_{j,k}^{(t)}). Будући да термин у загради маргинализује по свим могућим случајевима кад k\ne i, је једнако 1. Тако су коефицијенти сваког члана T_{1,i}^{(t)} и T_{2,i}^{(t)} који теже једнакости.

Корак М[уреди]

Квадратни облик Q(θ|θ(t)) значи да одређивање максималне вредности θ је релативно једноставно. Имајте на уму да τ, (μ1,Σ1) и (μ2,Σ2) могу бити све маскимизовано независно једни од других, јер се сви они појављују у одвојеним линеарним члановима. За почетак, узмимо у обзир τ, које има ограничење τ1 + τ2=1:

\begin{align}\boldsymbol{\tau}^{(t+1)}
&= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}}\  Q(\theta | \theta^{(t)} ) \\
&= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}} \ \left\{ \left[  \sum_{i=1}^n T_{1,i}^{(t)} \right] \log \tau_1 + \left[  \sum_{i=1}^n T_{2,i}^{(t)} \right] \log \tau_2  \right\}
\end{align}

Ово има исти облик као МЛЕ за биномну расподелу, тако да:

\tau^{(t+1)}_j = \frac{\sum_{i=1}^n T_{j,i}^{(t)}}{\sum_{i=1}^n (T_{1,i}^{(t)} + T_{2,i}^{(t)} ) } = \frac{1}{n} \sum_{i=1}^n T_{j,i}^{(t)}

У наредним проценама (μ1,Σ1):

\begin{align}(\boldsymbol{\mu}_1^{(t+1)},\Sigma_1^{(t+1)})
&= \underset{\boldsymbol{\mu}_1,\Sigma_1} {\operatorname{arg\,max}}\  Q(\theta | \theta^{(t)} ) \\
&= \underset{\boldsymbol{\mu}_1,\Sigma_1} {\operatorname{arg\,max}}\  \sum_{i=1}^n T_{1,i}^{(t)} \left\{ -\tfrac{1}{2} \log |\Sigma_1| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_1)^\top\Sigma_1^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_1) \right\}
\end{align}

Ово има исту форму као и тежак МЛЕ за нормалну расподелу, тако да

\boldsymbol{\mu}_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{1,i}^{(t)}} and \Sigma_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)})^\top }{\sum_{i=1}^n T_{1,i}^{(t)}}

и, преко симетрије:

\boldsymbol{\mu}_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{2,i}^{(t)}} and \Sigma_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)})^\top }{\sum_{i=1}^n T_{2,i}^{(t)}} .

Литература[уреди]

Спољашње везе[уреди]