EM algoritmus

EM algoritmus (z anglického expectation–maximization – očekávaná (střední) hodnota–maximalizace) je ve statistice iterační metoda pro hledání maximálně věrohodného odhadu nebo odhadu parametrů statistického modelu s maximální aposteriorní pravděpodobností (MAP), který závisí na nepozorovaných skrytých proměnných. Při EM iteracích se pravidelně střídají kroky výpočtu střední hodnoty (očekávání, E) s kroky maximalizace (M). V kroku E se vytváří očekávaná logaritmická věrohodnostní funkce na základě aktuálního odhadu parametrů. V kroku M se počítají parametry maximalizující očekávanou logaritmickou věrohodnostní funkci nalezenou v kroku E. Tyto odhady parametrů se pak používají pro určení rozdělení skrytých proměnných v dalším kroku E.

EM clusterování dat o erupcích gejzíru Old Faithful. Náhodný počáteční model (který kvůli různý měřítkům na osách vypadá jako dvě velmi nízké a široké elipsy) se přizpůsobuje pozorovaným datům. V první iteraci se model změní velmi výrazně, ale pak konverguje k rozlišení dvou režimů erupcí gejzíru. Vizualizováno pomocí ELKI.

Historie

EM algoritmus popsali, vysvětlili a pojmenovali Arthur P. Dempster, Nan Laird a Donald Rubin v klasickém článku z roku 1977, ve kterém zmínili, že EM algoritmus byl v minulosti „mnohokrát navržen a popsán pro speciální případy“ jinými autory. Jedním z prvních popisů byla metoda počítání genů pro odhadování frekvence alel, kterou popsal Cedric Smith. Velmi podrobný popis EM metody pro exponenciální funkce publikoval Rolf Sundberg ve své diplomové práci a několika odborných článcích vzniklých při jeho spolupráci s Per Martin-Löfem a Anders Martin-Löfem. Článek Dempstera, Lairda a Rubina z roku 1977 zobecňuje metodu a načrtává analýzu konvergence pro širší třídu problémů. Bez ohledu na dřívější objevy, inovativní článek Dempstera, Lairda a Rubina v Journal of Royal Statistical Society vyvolal nadšenou diskuzi na setkání Royal Statistical Society se Sundberg vedl k označení článku za „brilantní“. Zmiňovaný článek učinil z EM metody důležitý nástroj statistické analýzy.

Později se ukázalo, že analýza konvergence EM algoritmu ve výše zmíněném článku byl chybná; opravenou analýzu konvergence publikoval C. F. Jeff Wu v roce 1983. Jeho důkaz prokázal konvergenci EM metody i mimo rodinu exponenciálních rozdělení, kterou vyžadoval původní článek.

Úvod

EM algoritmus se používá pro hledání (lokální) maximálně věrohodných parametrů statistického modelu v případech, kdy rovnice nelze vyřešit přímo. Tyto modely typicky kromě neznámých parametrů a známých dat pozorování zahrnují skryté proměnné. To znamená, že buď v datech existují chybějící hodnoty nebo lze model formulovat jednodušeji, pokud předpokládáme, že existují další nepozorované datové body. Například smíšený model lze popsat jednodušeji, pokud předpokládáme, že ke každému pozorovanému datovému bodu existuje další nepozorovaný datový bod nebo skrytá proměnná, která popisuje kombinační komponentu, do které každý datový bod patří.

Hledání maximálně věrohodného řešení typicky vyžaduje použití derivace věrohodnostní funkce podle všech neznámých hodnot, parametrů a skrytých proměnných a současně řešení výsledné rovnice. U statistických modelů se skrytými proměnnými to obvykle není možné. Místo toho bývá výsledkem soustava vzájemně propojených rovnic, ve které řešení parametrů vyžaduje hodnoty skrytých proměnných a naopak, přičemž substituce jedné sady rovnic do druhé produkuje neřešitelné rovnice.

EM algoritmus vychází z pozorování, že tyto dvě sady rovnic lze řešit numericky: Je možné vybrat libovolné hodnoty pro jednu ze obou sad neznámých a pomocí nich odhadnout druhou sadu. Pak pomocí nových hodnot najít lepší odhad první sady, a tak pokračovat ve střídání obou kroků dokud výsledné hodnoty nekonvergují k pevným bodům. Lze dokázat, že v tomto kontextu tento postup konverguje, a že v tomto bodě je derivace věrohodnostní funkce libovolně blízko k nule, což znamená, že nalezený bod je buď maximem nebo sedlovým bodem. Obecně však může existovat více maxim, a není záruka, že algoritmus nalezne globální maximum. Některé věrohodnostní funkce mají v sobě také singularity, tj. nesmyslná maxima. Například jedno z řešení, které může EM algoritmus nalézt ve smíšeném modelu, spočívá v tom, že jedna ze složek má nulový rozptyl a střední parametr pro stejnou složku se rovná jednomu z datových bodů.

Popis

Je-li dán statistický model, který generuje množinu $\mathbf {X}$ pozorovaných dat, množina nepozorovaných skrytých dat nebo chybějících hodnot $\mathbf {Z}$ a vektor neznámých parametrů ${\boldsymbol {\theta }}$, spolu s věrohodnostní funkcí $L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )=p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})$, pak maximálně věrohodný odhad (MLE) neznámých parametrů je určen maximalizací marginální věrohodnosti pozorovaných dat

$L({\boldsymbol {\theta }};\mathbf {X} )=p(\mathbf {X} \mid {\boldsymbol {\theta }})=\int p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})\,d\mathbf {Z}$

Tato hodnota je však často nevyčíslitelná (například pokud $\mathbf {Z}$ je posloupnost událostí, poroste počet hodnot exponenciálně s délkou posloupnosti, takže přesný výpočet sumy bude extrémně obtížný).

EM algoritmus se snaží najít MLE marginální věrohodnosti iterativním používáním následujících dvou kroků:

E krok (Expectation): Definujeme $Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})$ jako očekávanou hodnotu logaritmické věrohodnostní funkce ${\boldsymbol {\theta }}$, vzhledem k aktuálnímu podmíněnému rozdělení pravděpodobnosti $\mathbf {Z}$ pro dané $\mathbf {X}$ a aktuální odhady parametrů ${\boldsymbol {\theta }}^{(t)}$:

$Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)}}\left[\log L({\boldsymbol {\theta }};\mathbf {X} ,\mathbf {Z} )\right]\,$

M krok (Maximalizace): Najdeme parametry, které maximalizují následující výraz:

${\boldsymbol {\theta }}^{(t+1)}={\underset {\boldsymbol {\theta }}{\operatorname {arg\,max} }}\ Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\,$

Typické modely, na které se EM aplikuje, používají $\mathbf {Z}$ jako skrytou proměnnou, což udává příslušnost k jedné z následujících situací:

1. Pozorované datové body $\mathbf {X}$ mohou být diskrétní (hodnoty tvoří konečnou nebo spočetně nekonečnou množinu) nebo spojité (hodnoty tvoří nespočetně nekonečnou množinu). Každému bodu dat může být přiřazen vektor pozorování.
2. Chybějící hodnoty (nazývané také skryté proměnné) $\mathbf {Z}$ jsou diskrétní, vybírané z pevného počtu hodnot a s jednou skrytou proměnnou na pozorovanou jednotku.
3. Parametry jsou spojité a jsou dvou druhů: Parametry, které jsou přiřazené k všem datovým bodům, a parametry přiřazené určité hodnotě skryté proměnné (tj. přiřazené ke všem datovým bodům, jejichž odpovídající skrytá proměnná má tuto hodnotu).

EM je však možné aplikovat i na jiné třídy modelů.

Motivace je následující: Pokud jsou známé hodnoty parametrů ${\boldsymbol {\theta }}$, lze hodnotu skryté proměnné $\mathbf {Z}$ obvykle nalézt maximalizací logaritmické věrohodnostní funkce přes všechny možné hodnoty $\mathbf {Z}$, buď jednoduše iterováním přes $\mathbf {Z}$ nebo pomocí nějakého algoritmu, například Baumova–Welchova algoritmu pro skryté Markovovy modely. Pokud naopak známe hodnoty skryté proměnné $\mathbf {Z}$, můžeme nalézt odhad parametrů ${\boldsymbol {\theta }}$ docela snadno, typicky jednoduchým seskupováním pozorovaných datových bodů podle hodnoty příslušných skrytých proměnných a průměrováním hodnot, nebo nějaké funkce hodnot, bodů z každé skupiny. V případě, když ${\boldsymbol {\theta }}$ i $\mathbf {Z}$ jsou neznámé, nabízí se použít iterativní algoritmus:

1. Nejdříve inicializovat parametry ${\boldsymbol {\theta }}$ na nějaké náhodné hodnoty.
2. Spočítat pravděpodobnost každé možné hodnoty $\mathbf {Z}$ pro dané ${\boldsymbol {\theta }}$.
3. Pak, použít právě vypočtené hodnoty $\mathbf {Z}$ pro vypočítání lepšího odhadu parametrů ${\boldsymbol {\theta }}$.
4. Opakovat kroky 2 a 3 dokud nekonvergují.

Právě popsaný algoritmus se monotónně přibližuje lokálnímu minimu nákladové funkce.

Vlastnosti

Označení krok očekávání (E) je poněkud nevhodné pojmenování. V prvním kroku se počítají pevné parametry funkce Q (závislé na datech). Jakmile jsou parametry Q známé, je Q plně určena a ve druhém (M) kroku EM algoritmu se provádí její maximalizace.

Ačkoli EM iterace skutečně zvyšuje (marginální) hodnotu věrohodnostní funkce na pozorovaných datech, není záruka, že posloupnost konverguje k odhadu maximální věrohodnosti. Pro multimodální distribuce to znamená, že EM algoritmus může konvergovat k lokálnímu maximu pozorované věrohodnostní funkce, v závislosti na počátečních hodnotách. Pro únik z lokálního maxima existuje množství heuristik a metaheuristik, například gradientní algoritmus s opakovaným náhodným startem (začíná s několika různými náhodnými počátečními odhady θ^(t)) nebo použití metod simulovaného žíhání.

EM algoritmus je zvláště užitečný, když věrohodnost patří do rodiny exponenciálních rozdělení: v kroku E se použije suma očekávání dostačujících statistik a v kroku M se maximalizuje lineární funkce. V takovém případě je obvykle možné odvodit tvaru analytické řešení aktualizací v každém kroku pomocí Sundbergova vzorce (který publikoval Rolf Sundberg na základě nepublikovaných výsledků Per Martin-Löfa a Anders Martin-Löfa).

V původním článku od Dempstera, Lairda a Rubina byla uvedena úprava EM metody pro výpočet maximálního aposteriorního odhadu (MAP) pro bayesovské vyvozování.

Existují i jiné metody pro hledání maximálně věrohodných odhadů, jako například metoda gradientního spádu, metoda sdružených gradientů nebo varianty Gaussova–Newtonova algoritmu. Na rozdíl od EM však tyto metody typicky vyžadují vyhodnocování první nebo druhé derivace věrohodnostní funkce.

Důkaz korektnosti

EM algoritmus se snaží zlepšit $Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})$, místo toho, aby přímo zlepšoval $\log p(\mathbf {X} \mid {\boldsymbol {\theta }})$. Ukážeme, že vylepšení prvního znamená také vylepšení druhého.

Pro jakékoli $\mathbf {Z}$ s nenulovou pravděpodobností $p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})$, můžeme psát

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

Vezměme očekávání přes možné hodnoty neznámých dat $\mathbf {Z}$ podle aktuálního odhadu parametru $\theta ^{(t)}$, znásobíme obě strany $p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})$ a sečteme (nebo zintegrujeme) podle $\mathbf {Z}$. Levá strana je očekávaná hodnota konstanty, takže dostáváme:

${\begin{aligned}\log p(\mathbf {X} \mid {\boldsymbol {\theta }})&=\sum _{\mathbf {Z} }p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})\log p(\mathbf {X} ,\mathbf {Z} \mid {\boldsymbol {\theta }})-\sum _{\mathbf {Z} }p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }}^{(t)})\log p(\mathbf {Z} \mid \mathbf {X} ,{\boldsymbol {\theta }})\\&=Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})+H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\,,\end{aligned}}$

kde $H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})$ je definované jako opačná hodnota součtu, který nahrazuje. Tato poslední rovnice platí pro každou hodnotu ${\boldsymbol {\theta }}$, včetně ${\boldsymbol {\theta }}={\boldsymbol {\theta }}^{(t)}$,

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

a odečtením této poslední rovnice od předchozí rovnice dává

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

Jensenova nerovnost nám však říká, že $H({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})\geq H({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})$, z čehož můžeme vyvodit, že

$\log p(\mathbf {X} \mid {\boldsymbol {\theta }})-\log p(\mathbf {X} \mid {\boldsymbol {\theta }}^{(t)})\geq Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})-Q({\boldsymbol {\theta }}^{(t)}\mid {\boldsymbol {\theta }}^{(t)})\,.$

Neboli výběr ${\boldsymbol {\theta }}$ pro zlepšení $Q({\boldsymbol {\theta }}\mid {\boldsymbol {\theta }}^{(t)})$ způsobí, že se $\log p(\mathbf {X} \mid {\boldsymbol {\theta }})$ zlepší nejméně o totéž.

Jako postup maximalizace maximalizace

EM algoritmus lze chápat jako dva střídající se kroky maximalizace, tj. jako příklad vzestupu podle souřadnice. Uvažujme funkci:

$F(q,\theta ):=\operatorname {E} _{q}[\log L(\theta ;x,Z)]+H(q),$

kde q je libovolné rozdělení pravděpodobnosti nepozorovaných dat z a H(q) je entropie rozdělení q. Tuto funkci lze zapsat jako

$F(q,\theta )=-D_{\mathrm {KL} }{\big (}q\parallel p_{Z\mid X}(\cdot \mid x;\theta ){\big )}+\log L(\theta ;x),$

kde $p_{Z\mid X}(\cdot \mid x;\theta )$ je podmíněné rozdělení pravděpodobnosti nepozorovaných dat, jsou-li dána pozorovaná data $x$, a $D_{KL}$ je Kullbackova–Leiblerova divergence.

Pak na kroky v EM algoritmu můžeme pohlížet jako:

Expectation krok: Zvolíme $q$, které maximalizuje $F$:

$q^{(t)}=\operatorname {arg\,max} _{q}\ F(q,\theta ^{(t)})$

Maximalizační krok: Zvolíme $\theta$, které maximalizuje $F$:

$\theta ^{(t+1)}=\operatorname {arg\,max} _{\theta }\ F(q^{(t)},\theta )$

Aplikace

EM se často používá pro shlukovou analýzu ve strojovém učení a počítačovém vidění. Při zpracování přirozeného jazyka se velmi často používají dvě varianty EM algoritmu – Baumův–Welchův algoritmus pro skryté Markovovy modely a inside-outside algoritmus pro učení pravděpodobnostních bezkontextových gramatik bez učitele.

EM se často používá pro odhad parametrů smíšených modelů, především v kvantitativní genetice.

V psychometrice je EM algoritmus téměř nepostradatelný pro odhadování parametrů položek a skrytých schopností modelů teorie odpovědi na položku.

Díky schopnosti pracovat s chybějícími daty a pozorovat neidentifikované proměnné se EM stává užitečným nástrojem pro oceňování a správu rizik portfolia.

EM algoritmus (a jeho rychlejší varianta OSEM Ordered subset expectation maximization) se často používá při vytváření obrazu v zobrazovacích metodách v lékařství, zvláště pro pozitronovou emisní tomografii a jednofotonovou emisní výpočetní tomografii. Rychlejší varianty EM jsou popsány dále.

Ve strukturálním inženýrství je algoritmus STRIDE Structural Identification using Expectation Maximization pouze výstupní metoda pro identifikaci přirozených vibračních vlastností strukturálního systému pomocí senzorových dat (viz operační modální analýza).

EM algoritmy pro filtrování a vyhlazování

Kalmanův filtr se typicky používá pro online odhad stavu a vyhlazování s minimálním rozptylem může být použito pro offline nebo dávkový odhad stavu. Nicméně tato řešení s minimálním rozptylem vyžadují odhady parametrů modelu stavového prostoru. EM algoritmy lze používat pro řešení problémů sdružených stavů a odhady parametrů.

Filtrovací a vyhlazovací EM algoritmy vznikají opakováním následujícího dvoukrokového postupu:

krok E

Použít Kalmanův filtr nebo vyhlazování s minimálním rozptylem navržený s odhady aktuálních parametrů pro získání odhadů aktualizovaných stavů.

krok M

Použití filtrovaných nebo vyhlazených odhadů stavu pro výpočet maximální věrohodnosti pro získání odhadů aktualizovaných parametrů.

Předpokládejme, že Kalmanův filtr nebo vyhlazování s minimálním rozptylem pracuje s měřeními systému s jediným vstupem a jediným výstupem, který je ovlivněn aditivním bílým šumem. Aktualizovaný odhad měření rozptylu šumu lze získat z výpočtu maximální věrohodnosti

${\widehat {\sigma }}_{v}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{(z_{k}-{\widehat {x}}_{k})}^{2}$

kde ${\widehat {x}}_{k}$ jsou odhady skalárního výstupu vypočítané filtrem nebo vyhlazováním z N skalárních měření $z_{k}$. Výše uvedená aktualizace může také být aplikována na aktualizaci Poissonovo měření intenzity šumu. Podobně pro autoregresivní proces prvního řádu lze vypočítat odhad rozptylu aktualizovaného šumu procesu podle vzorce

${\widehat {\sigma }}_{w}^{2}={\frac {1}{N}}\sum _{k=1}^{N}{({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {X}}_{k})}^{2}$

kde ${\widehat {x}}_{k}$ a ${\widehat {x}}_{k+1}$ jsou skalární odhady stavu vypočítané filtrem nebo vyhlazování. Aktualizovaný odhad koeficientů modelu je získaný podle vzorce

${\widehat {F}}={\frac {\sum _{k=1}^{N}({\widehat {x}}_{k+1}-{\widehat {F}}{\widehat {x}}_{k})}{\sum _{k=1}^{N}{\widehat {x}}_{k}^{2}}}.$

Konvergence odhadů parametrů obdobných jako jsou uvedeny výše jsou dobře studované.

Varianty

Bylo navrženo několik metod pro zrychlení někdy pomalé konvergence EM algoritmu, například metody používající sdružený gradient a modifikace Newtonovy metody (Newtonova–Raphsonova metoda). EM lze také používat pro metody omezeného odhadu.

Algoritmus Parameter-expanded expectation maximization (PX-EM) vede často ke zrychlení díky „použití `nastavení kovariance' pro opravu analýzy kroku M, přičemž využívá zvláštní informace zachycené v přisuzovaných úplných datech“.

Expectation conditional maximization (ECM) nahrazuje každý M krok v posloupností kroků podmíněné maximalizace (CM), ve kterých se každý parametr θ_i maximalizuje samostatně, když se ostatní parametry nemění. Rozšířením tohoto algoritmu je Expectation conditional maximization either (ECME).

Uvedený přístup lze dále rozšířit na algoritmus generalized expectation maximization (GEM), u kterého se hledají pouze přírůstky účelové funkce F jak v kroku E tak v kroku M, jak je popsané v části Jako postup maximalizace maximalizace. GEM se dále vyvíjí v distribuovaném prostředí a vykazuje slibné výsledky.

Je také možné považovat EM algoritmus za podtřídu MM algoritmu (Majorize/Minimize nebo Minorize/Maximize, podle kontextu), a používat pro něj všechny techniky vyvinuté pro obecnější případ.

α-EM algoritmus

Funkce Q používaná v EM algoritmu je založena na logaritmické věrohodnostní funkci. Proto je považována za log-EM algoritmus. Použití logaritmické věrohodnostní funkce lze zobecnit na použití α-logaritmu poměru věrohodností. Potom lze α-log poměru věrohodností pozorovaných dat přesně vyjádřit jako rovnost použitím funkce Q aplikované na α-logaritmus poměru věrohodností a α-divergence. Získání této funkce Q je zobecněný krok E. Jeho maximalizace je zobecněný krok M. Výsledný algoritmus se nazývá α-EM algoritmus. α-EM algoritmus, jehož autorem je Yasuo Matsuyama je zobecněním log-EM algoritmu. Není potřeba žádný výpočet gradientu nebo matice Hessianů. Při výběru vhodného α vykazuje α-EM algoritmus rychlejší konvergenci než log-EM. α-EM algoritmus vede k rychlejší verzi algoritmu α-HMM odhadu skrytých Markovových modelů.

Vztah k variačním bayesovským metodám

EM je částečně nebayesovská metoda založená na maximální věrohodnosti. Její konečný výsledek dává rozdělení pravděpodobnosti přes skryté proměnné (v bayesovském stylu) spolu s bodovým odhadem θ (buď maximálně věrohodný odhad anebo v a posterior režimu). Může být požadována jeho plně bayesovská verze, která dává rozdělení pravděpodobnosti podle θ a skryté proměnné. Při bayesovském přístupu k inference bychom jednoduše považovali θ za další skrytou proměnnou. Při tomto přístupu mizí rozdíl mezi kroky E a M. Při použití faktorizované aproximace Q popsané výše (variační bayesovská metoda) může řešení iterovat přes každou skrytou proměnnou (včetně θ) a optimalizovat je po jedné. Nyní je potřeba k kroků pro každou iteraci, kde k je počet skrytých proměnných. Pro grafické modely je to snadné udělat, protože nové Q pro každou proměnnou závisí pouze na jeho Markovově blanketu, takže pro efektivní inference lze používat lokální předávání variačních zpráv (Variational message passing, VMP).

Geometrická interpretace

V informační geometrii jsou kroky E a M interpretovány jako projekce pod duální afinní konexe, nazývané e-spojení a m-spojení; Pomocí těchto termínů lze pohlížet na Kullbackovu–Leiblerovu divergenci.

Příklady

Smíšené normální rozdělení

Porovnání k-means clusteringu a EM na umělá data vizualizované pomocí ELKI. Při použití variancí může EM algoritmus popsat normální rozdělení přesně, zatímco k-means rozděluje data na Voroného buňky. Střed clusteru je znázorněn světlejším větším symbolem.

Animovaná demonstrace, jak se EM algoritmus přizpůsobuje dvěma komponentám gaussovského smíšeného modelu pro datový soubor Old Faithful. Algoritmus prochází od náhodné inicializace po konvergenci.

Nechť $\mathbf {x} =(\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{n})$ jsou vzorky $n$ nezávislých pozorování z kombinace dvou vícerozměrných normálních rozdělení dimenze $d$ a nechť $\mathbf {z} =(z_{1},z_{2},\ldots ,z_{n})$ jsou skryté proměnné, které určují komponentu, z nichž pozorování pochází.

$X_{i}\mid (Z_{i}=1)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{1},\Sigma _{1})$ a $X_{i}\mid (Z_{i}=2)\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{2},\Sigma _{2})$

kde

$\operatorname {P} (Z_{i}=1)=\tau _{1}\,$ a $\operatorname {P} (Z_{i}=2)=\tau _{2}=1-\tau _{1}$

Cílem je odhadnout neznámé parametry reprezentující mísení hodnot dvou normálních rozdělení a střední hodnoty a kovariance obou:

$\theta ={\big (}{\boldsymbol {\tau }},{\boldsymbol {\mu }}_{1},{\boldsymbol {\mu }}_{2},\Sigma _{1},\Sigma _{2}{\big )}$

kde Věrohodnostní funkce věrohodnosti pro neúplná data je

$L(\theta ;\mathbf {x} )=\prod _{i=1}^{n}\sum _{j=1}^{2}\tau _{j}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{j})$,

a Věrohodnostní funkce věrohodnosti pro úplná data je

$L(\theta ;\mathbf {x} ,\mathbf {z} )=p(\mathbf {x} ,\mathbf {z} \mid \theta )=\prod _{i=1}^{n}\prod _{j=1}^{2}\ [f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{j})\tau _{j}]^{\mathbb {I} (z_{i}=j)}$

nebo

$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\}.$

kde $\mathbb {I}$ je charakteristická funkce a $f$ je hustota pravděpodobnosti vícerozměrného normálního rozdělení.

V poslední rovnosti se pro každé i jedna hodnota charakteristické funkce $\mathbb {I} (z_{i}=j)$ rovná nule a jedna jedničce. Vnitřní součet tedy má vždy pouze jeden člen.

Krok E

Máme-li aktuální odhad parametrů θ^(t), bude podmíněné rozdělení pravděpodobnosti Z_i určeno Bayesovou větou, aby to byla proporcionální výška normální hustoty vážená parametrem τ:

$T_{j,i}^{(t)}:=\operatorname {P} (Z_{i}=j\mid 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)})}}$.

Tyto hodnoty, které jsou obvykle považovány za výstup kroku E (přestože to není funkce Q popsaná níže), se nazývají „členské pravděpodobnosti“.

Tento E krok odpovídá nastavení této funkce pro Q:

${\begin{aligned}Q(\theta \mid \theta ^{(t)})&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,\mathbf {\theta } ^{(t)}}[\log L(\theta ;\mathbf {x} ,\mathbf {Z} )]\\&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,\mathbf {\theta } ^{(t)}}[\log \prod _{i=1}^{n}L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\operatorname {E} _{\mathbf {Z} \mid \mathbf {X} ,\mathbf {\theta } ^{(t)}}[\sum _{i=1}^{n}\log L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\sum _{i=1}^{n}\operatorname {E} _{Z_{i}\mid \mathbf {X} ;\mathbf {\theta } ^{(t)}}[\log L(\theta ;\mathbf {x} _{i},Z_{i})]\\&=\sum _{i=1}^{n}\sum _{j=1}^{2}P(Z_{i}=j\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})\log L(\theta _{j};\mathbf {x} _{i},j)\\&=\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{aligned}}$

Střední hodnota („očekávání“) $\log L(\theta ;\mathbf {x} _{i},Z_{i})$ v sumě se bere vzhledem k hustotě pravděpodobnosti $P(Z_{i}\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)})$, která může být pro každé $\mathbf {x} _{i}$ z trénovací množiny jiná. Před provedením kroku E jsou známé všechny potřebné hodnoty kromě $T_{j,i}$, které se počítají podle rovnice na začátku kroku E.

Tuto plně podmíněnou střední hodnotu není třeba počítat v jednom kroku, protože τ a μ/Σ se vyskytují ve zvláštních lineárních členech a mohou být tedy maximalizovány nezávisle.

Krok M

Protože Q(θ | θ^(t)) je kvadratická funkce, je určování maximalizujících hodnot θ relativně přímočaré. Také τ, (μ ₁,Σ₁) a (μ ₂,Σ₂) mohou být vesměs maximalizovány nezávisle, protože se všechny vyskytují v samostatných lineárních členech.

Pro začátek uvažujme τ s omezením τ₁ + τ₂=1:

${\begin{aligned}{\boldsymbol {\tau }}^{(t+1)}&={\underset {\boldsymbol {\tau }}{\operatorname {arg\,max} }}\ Q(\theta \mid \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{aligned}}$

což má stejný tvar jako MLE pro binomické rozdělení, takže

$\tau _{j}^{(t+1)}={\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)}$.

Pro další odhady (μ ₁,Σ₁):

${\begin{aligned}({\boldsymbol {\mu }}_{1}^{(t+1)},\Sigma _{1}^{(t+1)})&={\underset {{\boldsymbol {\mu }}_{1},\Sigma _{1}}{\operatorname {arg\,max} }}\ Q(\theta \mid \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{aligned}}$.

což má stejný tvar jako vážené MLE pro normální rozdělení, takže

${\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)}}}$ a $\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)}}}$

a díky symetrii

${\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)}}}$ a $\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)}}}$.

Ukončení

Iterativní proces ukončíme, jakmile $E_{Z\mid \theta ^{(t)},\mathbf {x} }[\log L(\theta ^{(t)};\mathbf {x} ,\mathbf {Z} )]\leq E_{Z\mid \theta ^{(t-1)},\mathbf {x} }[\log L(\theta ^{(t-1)};\mathbf {x} ,\mathbf {Z} )]+\varepsilon$ pro předem zvolenou hodnotu $\varepsilon$.

Zobecnění

Výše popsaný algoritmus lze zobecnit pro kombinaci více než dvou vícerozměrných normálních rozdělení.

Zkrácená a cenzurovaná regrese

EM algoritmus byl implementován pro případ, kdy existuje podkladový model lineární regrese, který vysvětluje variace určité velikosti, ale skutečně pozorované hodnoty jsou cenzurovanou nebo zkrácenou verzí hodnot reprezentovaných v modelu. Speciální případy tohoto modelu zahrnují cenzurované nebo zkrácené pozorování z jednoho normální rozdělení.

Alternativy

EM algoritmus typicky konverguje k lokální optimu. Toto lokální optimum však nemusí být globálním optimem, a také nemáme jakýhokoli odhad rychlosti konvergence. V mnoharozměrném prostoru může být konvergence libovolně špatná a může existovat exponenciální počet lokálních optim. Proto existuje potřeba pro alternativní metody pro zaručené učení, zvláště v mnohadimanzionálních prostředích. Vedle algoritmu EM existují i jiné algoritmy, které lépe zaručují konzistenci. Tyto algoritmy využívají přístupy založené na momentech nebo tak zvané spektrální techniky. Zájem o algoritmy pro učení parametrů pravděpodobnostního modelu používající momenty se v poslední době stále zvyšuje, protože za určitých podmínek mohou na rozdíl od EM algoritmu, jehož velkým problémem je uváznutí v lokálním optimu, zaručit například globální konvergenci. Algoritmy se zaručeným učením mohou být odvozeny pro několik důležitých modelů, jako jsou smíšené modely, skryté Markovovy modely atd. V těchto spektrálních metodách, se neobjevují žádná nepravá lokální optima, a skutečné parametry lze při splnění určitých podmínek regularity konzistentně odhadnout.

Odkazy

Reference

V tomto článku byl použit překlad textu z článku Expectation–maximization algorithm na anglické Wikipedii.

Související články

• Kombinační rozdělení
• Složené rozdělení
• Odhad hustoty rozdělení
• Spektroskopie s úplnou absorpcí
• MM algoritmus

Literatura

• HOGG, Robert; MCKEAN, Joseph; CRAIG, Allen. Introduction to Mathematical Statistics. Upper Saddle River, NJ: Pearson Prentice Hall, 2005. 
• DELLAERT, Frank. The Expectation Maximization Algorithm. [s.l.]: [s.n.], 2002.  dává jednodušší vysvětlení EM algoritmu jako maximalizace se spodní mezí.
• BISHOP, Christopher M., 2006. Pattern Recognition and Machine Learning. [s.l.]: Springer. Dostupné online. ISBN 978-0-387-31073-2. 
• GUPTA, M. R.; CHEN, Y. Theory and Use of the EM Algorithm. Foundations and Trends in Signal Processing. 2010. DOI 10.1561/2000000034.  Dobře napsaná krátká kniha o EM včetně podrobného odvození EM pro GMMs, skryté Markovovy modely a Dirichlet.
• BILMES, Jeff. A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models. [s.l.]: [s.n.], 1998.  obsahuje zjednodušené odvození EM rovnice pro gaussovské Kombinace a gaussovské Kombinace skrytých Markovových Modelů.
• Einicke, G.A. Smoothing, Filtering and Prediction: Estimating the Past, Present and Future. Rijeka, Croatia: Intech, 2012. Dostupné online. ISBN 978-953-307-752-9. DOI 10.5772/2706. 
• MCLACHLAN, Geoffrey J.; KRISHNAN, Thriyambakam. The EM Algorithm and Extensions. 2. vyd. Hoboken: Wiley, 2008. ISBN 978-0-471-20170-0. 

Externí odkazy

• Různé 1D, 2D a 3D ukázky EM se smíšenými modely poskytnutý jako materiál pro SOCR Statistics Online Computational Resource činnosti a aplety. Tyto aplety a aktivity ukazují empirické vlastnosti EM algoritmu pro odhad parametrů při různých nastaveních.
• k-MLE: Rychlý algoritmus učení statistickým smíšeným modelům
• Hierarchie tříd v C++ (GPL) včetně gaussovských kombinací
• The on-line textbook: Information Theory, Inference, and Learning Algorithms Online učebnice: informační teorie, inferenční a učicí algoritmy, od Davida MacKaye obsahuje jednoduché příklady EM algoritmu jako například clustering pomocí měkkého k-means algoritmu a zdůrazňuje variační pohled na EM algoritmus, jak je popsané v kapitole 33.7 verze 7.2 (čtvrté vydání).
• Variační Algoritmy pro přibližné bayesovské vyvozování, M. J. Beala obsahuje porovnání EM na variačních bayesovských EM a odvození několika modelů včetně variačních bayesovských skryté Markovovy modely (kapitoly).
• Expectation Maximization Algoritmus: Krátký tutorial, kompletní odvození EM algoritmu od Seana Bormana.
• EM Algoritmus, autor: Xiaojin Zhu.
• EM algoritmus a varianty: neformální tutoriál Alexise Roche. Podrobný a velmi jasný popis EM algoritmu a mnoha zajímavých variant.

Portály: Matematika