Dlaczego algorytm EM działa (ang. Expectation–maximization)?

Algorytm EM ma wiele zastosowań, w sytuacji gdy bezpośrednia maksymalizacji funkcji log-wiarogodności l(θ,Y)l(\theta, Y) jest trudna. W wielu sytuacjach, problem maksymalizacji można uprościć, jeżeli przestrzeń parametrów rozszerzy się o dodatkowe - nieobserwowane zmienne.

Typowe przykłady takich sytuacji to uzupełnianie braków w danych oraz badanie mieszanin rozkładów.

Poniższy opis rozpoczniemy od ogólnego sformułowania algorytmu EM, a następnie przedstawimy konkretny przykład. Ostatni punkt to przedstawienie zastosowania algorytmu EM do maksymalizacji rozkładu a-posteriori.

Ogólne sformułowanie problemu

Oznaczmy przez YY obserwowaną zmienną losową (potencjalnie wielowymiarową) o rozkładzie FF z rodziny indeksowanej parametrem θ\theta. Bezpośrednia estymacja θ\theta metodą największej wiarygodności jest trudna, więc rozszerzamy przestrzeń o zmienne ukryte/nieobserwowane ZZ.

Łączny wektor zmiennych losowych (Y,Z)(Y, Z) pochodzi z rodziny indeksowanej parametrem θ∈Θ\theta \in \Theta. Gęstość tego rozkładu oznaczmy przez f(y,z;θ)f(y, z; \theta).

Algorytm EM

Algorytm EM ma następującą postać:

  1. Wybierz początkową wartość wektora parametrów θ^(1)\hat\theta^{(1)},

  2. Krok E (Expectation). Wyznacz warunkową wartość oczekiwaną dla aktualnego θ^(j)\hat\theta^{(j)}, gdzie jj to numer iteracji. Q(θ,θ^(j))=E(l(θ;Z,Y)∣Y,θ^(j)). Q\left(\theta, \hat\theta^{(j)}\right) = E\left(l(\theta; Z, Y) | Y, \hat\theta^{(j)}\right). Będziemy ją traktować jako funkcję po θ\theta.

  3. Krok M (Maximization). Wyznacz kolejną wartość θ^(j+1)\hat\theta^{(j+1)}, taką która maksymalizuje Q(θ,θ^(j))Q\left(\theta, \hat\theta^{(j)}\right) po θ\theta. θ^(j+1)=argmaxθQ(θ,θ^(j)). \hat\theta^{(j+1)} = arg \max_\theta Q\left(\theta, \hat\theta^{(j)}\right).

  4. Powtarzaj kroki 2-3 do spełnienia określonego warunku stop.

Dlaczego to działa?

Pytanie, na które musimy odpowiedzieć, to dlaczego optymalizacja Q(θ,θ^(j))Q\left(\theta, \hat\theta^{(j)}\right) pomaga w maksymalizacji funkcji wiarygodności.

Zacznijmy od rozkładu warunkowego P(Z∣Y;θ)=P(Z,Y;θ)P(Y;θ), P(Z|Y; \theta) = \frac{ P(Z, Y; \theta) }{ P(Y; \theta) }, a więc P(Y;θ)=P(Z,Y;θ)P(Z∣Y;θ). P(Y; \theta) = \frac{P(Z,Y; \theta)}{P(Z|Y; \theta)}.

Interesuje nas funkcja log-wiarogodności. Korzystając z powyższego, można ją wyrazić jako l(θ;Y)=l0(θ;Z,Y)−l1(θ;Z∣Y). l(\theta; Y) = l_0(\theta; Z, Y) - l_1(\theta; Z|Y). Zauważmy, że funkcje wiarygodności l()l(), l0()l_0() i l1()l_1() dotyczą różnych zmiennych, stąd różne oznaczenia.

Nałóżmy teraz obustronnie warunkową wartość oczekiwaną E(.∣Y,θ′)E(. | Y, \theta'). Parametr θ′\theta' będziemy w przyszłości wykorzystywali jako ustalony, roboczy, suboptymalny wektor parametrów. Otrzymujemy l(θ;Y)=E(l0(θ;Z,Y)∣Y,θ′)−E(l1(θ;Z∣Y)∣Y,θ′) l(\theta; Y) = E(l_0(\theta; Z, Y) | Y, \theta') - E(l_1(\theta; Z|Y) | Y, \theta')

Aby uprościć zapis poniżej, przyjmijmy takie standardowe oznaczenia dla prawej strony tego równania. l(θ;Y)=Q(θ,θ′)−R(θ,θ′) l(\theta; Y) = Q(\theta, \theta') - R(\theta, \theta')

Rozłożyliśmy funkcję log-wiarogodności, którą chcielibyśmy maksymalizować po θ\theta na dwa człony. Algorytm EM będzie bezpośrednio maksymalizował Q(θ,θ′)Q(\theta, \theta').

Funkcja R(θ,θ′)R(\theta, \theta') ma swoje maksimum gdy wartość oczekiwana jest warunkowana tym samym rozkładem co będąca pod nią funkcja wiarygodności. Jest tak dlatego, że funkcja wiarygodności jest wypukła i można wykorzystać nierówności Jensena. Pośrednie kroki uzasadnienia tej właściwości są opisane w ćwiczeniu 8.2 Elements of Statistical Learning.

A więc R(θ,θ′)≤R(θ′,θ′)R(\theta, \theta') \leq R(\theta', \theta'), co z kolei oznacza, że jeżeli jesteśmy w stanie znaleźć θ\theta, które zwiększy Q(θ,θ′)Q(\theta, \theta') to z pewnością zwiększymy też l(θ;Y)l(\theta; Y).

A może tak policzyć coś na palcach?

Ogólne sformułowanie algorytmu EM umożliwia stosowanie go do rozmaitych sytuacji. Ale aby wyrobić sobie intuicję przeliczmy jeden prosty przykład, jednowymiarowej mieszaniny dwóch rozkładów normlanych różniących się tylko średnią. Dla większej liczby składowych, wymiarów i parametrów obliczenia przeprowadza się według tego samego schematu, jest jedynie więcej symboli.

Model mieszaniny. Dwie składowe Y1∼N(μ1,σ2), Y_1 \sim \mathcal N(\mu_1, \sigma^2), Y2∼N(μ2,σ2), Y_2 \sim \mathcal N(\mu_2, \sigma^2), rozkład mieszający Z∼B(1,π), Z \sim B(1, \pi), i mieszanina Y=(1−Z)Y1+ZY2. Y = (1-Z) Y_1 + Z Y_2.

Mieszanina to w π\pi rozkład zmiennej Y1Y_1 a w 1−π1-\pi rozkład zmiennej Y2Y_2. Gęstość tego rozkładu można opisać jako fY(x)=(1−π)ϕ1(x)+πϕ2(x). f_Y(x) = (1-\pi) \phi_1(x) + \pi \phi_2(x).

Nasz wektor parametrów to θ=(π,μ1,μ2,σ)\theta=(\pi, \mu_1, \mu_2, \sigma) a funkcja log-wiarogodności ma postać l(θ,y)=∑ilog[(1−π)ϕ1(yi)+πϕ2(yi)]. l(\theta, y) = \sum_i \log [(1-\pi)\phi_1(y_i) + \pi \phi_2(y_i)]. Optymalizacja takiej sumy logarytmów sum nie jest prosta.

Zgodnie z duchem EM ułatwimy sobie zadanie dodając do modelu zmienne ukryte ZZ. Funkcja wiarogodności w takim modelu ma postać l(θ,y,z)=∑i[(1−zi)logϕ1(yi)+zilogϕ2(yi)]+∑i[(1−zi)logπ+zilog(1−π)]. l(\theta, y, z) = \sum_i [(1-z_i) \log \phi_1 (y_i) + z_i \log \phi_2 (y_i)] + \sum_i [(1-z_i) \log \pi + z_i \log (1-\pi) ].

Tę funkcję będzie nam łatwiej maksymalizować. Przyjrzyjmy się krokom w algorytmie EM.

Krok E

Chcemy wyznaczyć warunkową wartość oczekiwaną funkcji l(θ,y,z)l(\theta, y, z) po y,θ(i)y, \theta^{(i)}. Człony funkcji wiarogodności zawierające yiy_i się nie zmienią. Musimy jedynie policzyć co stanie się z członem ziz_i.

E(Zi∣Y;θ(i))=Pr(Zi=1∣Y;θ(i))=Pr(Zi=1,Y;θ(i))Pr(Y;θ(i)) E(Z_i | Y; \theta^{(i)}) = Pr(Z_i = 1 | Y; \theta^{(i)}) = \frac{Pr(Z_i = 1, Y; \theta^{(i)})}{Pr(Y; \theta^{(i)})} a więc E(Zi∣Y;θ(i))=π^(i)ϕ2(yi)(1−π^(i))ϕ1(yi)+π^(i)ϕ2(yi)=ηi. E(Z_i | Y; \theta^{(i)}) = \frac{\hat \pi^{(i)} \phi_2(y_i) }{(1-\hat \pi^{(i)}) \phi_1(y_i) + \hat \pi^{(i)} \phi_2(y_i) } = \eta_i.

Wyliczone wartości oczekiwane wstawiamy w miejsce ziz_i.

Krok M

Funkcja Q(θ,θ^(i))Q(\theta, \hat \theta^{(i)}) jest już funkcją parametrów (μ1,μ2,σ,π)(\mu_1, \mu_2, \sigma, \pi). Z uwagi na jej postać, można każdy z parametrów maksymalizować niezależnie. Wyznaczamy pochodną i przyrównujemy ją do zera.

Otrzymamy

μ^1(j+1)=∑i(1−ηi)yi∑i(1−ηi), {\hat\mu_1}^{(j+1)} = \frac{\sum_i (1-\eta_i)y_i}{\sum_i (1-\eta_i)}, μ^2(j+1)=∑iηiyi∑iηi, {\hat\mu_2}^{(j+1)} = \frac{\sum_i \eta_i y_i}{\sum_i \eta_i}, σ^22,(j+1)=[∑i(1−ηi)(yi−μ^1(j))2+∑iηi(yi−μ^2(j))2]/n, \hat \sigma_2^{2, (j+1)} = \left[\sum_i (1-\eta_i) \left(y_i - {\hat\mu_1}^{(j)}\right)^2 + \sum_i \eta_i \left(y_i - {\hat\mu_2}^{(j)}\right)^2\right]/n, π^(j+1)=∑iηi/n. \hat \pi^{(j+1)} = \sum_i \eta_i / n.

Kroki E i M należy powtarzać aż nie uzyska się przyzwoitej zbieżności.

Estymacja MAP

Przedstawiony powyżej algorytm EM maksymalizuje funkcje wiarogodności. Rosnąca popularność metod Bayesowskich spowodowała, że jego odmiana jest też wykorzystywana do znajdowania estymatorów maksimum a-posteriori (MAP).

Punktem wyjścia jest wzór na rozkład a-posteriori p(θ∣Y)≈l(θ,y)+log(p(θ)). p(\theta | Y) \approx l(\theta, y) + \log(p(\theta)).

Algorytm EM jest ogólnym algorytmem maksymalizacji. Możemy wykorzystać go również do szukania maksimum rozkładu aposteriori, dodając w kroku E człon log(p(θ))\log(p(\theta))

QMAP(θ,θ′)=Q(θ,θ′)+log(p(θ)). Q^{MAP}(\theta, \theta') = Q(\theta, \theta') + \log(p(\theta)).

Inne materiały

Alternatywne wyprowadzenie EM przez Andrew Ng http://cs229.stanford.edu/notes/cs229-notes8.pdf.