Az EM-algoritmust (expectation-maximization algorithm) először általánosságban tárgyaljuk, majd konkretizáljuk egy olyan feladatra, amely kevert minták szétválasztását célozza^[9].

Általánosságban az EM-algoritmust olyan feladatok megoldására fejlesztették ki, amelyekben ismeretlen, de konkrét, elvileg rögzített paraméterek becslését szeretnénk meghatározni valamilyen adathalmaz alapján úgy, hogy bizonyos adatelemeket nem ismerünk (pl. nem tudjuk ezeket mérni), csak valószínűség-eloszlásukról van információnk [Dem77], [Mit97].

Az EM-algoritmus szemléltetésére több példa is felhozható. Első példánk legyen a következő: meg akarjuk becsülni egy sokdimenziós adathalmaz átlagát úgy, hogy az adataink hiányosak, bizonyos adatvektorok egy vagy több komponense hiányzik. (Pl. egy – valamilyen elektromos vagy elektronikus berendezésekből álló – termékhalmazon mérték a termék áramfelvételét, két meghibásodás közötti időt stb., de néhány terméknél egyik-másik adat hiányzik, elveszett stb.) Egy másik lehetséges példa, amelyet későbbiekben részletesebben is tárgyalunk a következő: van egy adatsorunk, melyben az adatok több eloszlásból származhatnak. Az adatok alapján szeretnénk megbecsülni az egyes eloszlások várható értékét, de anélkül, hogy tudnánk, hogy az egyes adatok melyik eloszlásból származnak. Egy hétköznapi példaként gondoljunk arra, hogy az amerikaiak magasság, vérnyomás stb. adatai alapján azt szeretnénk megmondani, hogy az amerikaiakon belül az egyes népcsoportok (fehér, néger, mexikói, ázsiai) átlagmagassága, vérnyomásának átlaga és szórása stb. mekkora. Mindezt azonban úgy szeretnénk meghatározni, hogy az egyes méréseknél nem lett feljegyezve, hogy ki melyik népcsoportba tartozik. Van tehát egy olyan adathalmazunk, mely pl. a magasságértékeket tartalmazza, de az egyes értékek mellett nincs meg az az információ, hogy melyik magasságérték melyik népcsoportba tartozó ember magassága, vagyis hiányos adatokkal van dolgunk. Itt az ismeretlen paraméter, a hiányzó adat az, hogy melyik adat melyik csoportból származik, a keresett paramétervektor pedig a csoportok átlagos magasságából képzett vektor. (Természetesen ennek a demográfiai illetve népegészségügyi példának számos műszaki analógiáját lehet találni, de ez talán egy könnyen érthető példa.)

Az EM-algoritmust általános formájában a következőképpen szokták megfogalmazni (pl. [Mit97]). Egy valószínűség eloszlás valamilyen $\theta$ paramétervektorát akarjuk becsülni csupán részlegesen megfigyelhető adatok alapján. Legyen X={x₁, x₂, …, x_L} a megfigyelt adatkomponensek halmaza és Z={z₁, z₂, …, z_L} azon adatkomponensek halmaza, amelyeket nem tudtunk megfigyelni. Tehát az $Y=X\cup Z$ lenne a teljes adathalmaz, és ha mind X, mind Z ismert volna a feladatot könnyen meg lehetne oldani. Minthogy Z-t nem ismerjük, ezt hiányzó ismeretnek nevezzük. A nem ismert Z komponenseket valószínűségi változóként kezelhetjük, amelyek eloszlása az ismeretlen $\theta$ paramétervektortól függ. Hasonlóképpen Y is valószínűségi változó, hiszen X és Z segítségével definiáltuk. A Z hiányzó adatokat jellemző valószínűség eloszlás függ a keresett paramétertől. Ha ismernénk Z valószínűség eloszlását, akkor az X megfigyelt adatok alapján a $\theta$ paramétervektor becslését meg tudnánk határozni. Ugyanis a $\theta$ paramétervektor becsléséhez Z teljes ismerete nem szükséges Z várható értéke alapján a becslés meghatározható (Meg tudnánk fogalmazni egy maximum likelihood becslési feladatot, melynek megoldásaként $\theta$ olyan becslését határozhatnánk meg, mely mellett a megfigyléseink a legnagyobb valószínűségűek.) A $\theta$ becslés segítségével a hiányzó adatok valószínűség eloszlását pontosabban meg tudjuk adni, ami (a maximum likelihood becslésen keresztül) a $\theta$ paraméter jobb becslését eredményezi. Az eljárást folytatva a paramétervektor egyre jobb becslését kapjuk. Megfelelő feltételek mellett ez az iteratív eljárás konvergens és a $\theta$ becslések sorozata a valódi $\theta$ paraméterhez konvergál [Dem77].

Általánosabban az eljárást a következőképpen fogalmazhatjuk meg. Van egy kiinduló hipozézisünk az ismeretlen paraméterről. Jelöljük ezt h-val. A h aktuális hipotézis megad egy aktuálisan igaznak tartott paramétervektort. Egy iteráció során egy új hipotézist, h’-t állítunk fel, mely a paramétervektor egy jobb becslését adja. Olyan új hipotézist keresünk, mely mellett a (teljes) megfigyeléseinkre vonatkozó $P\left(Y|h^{\prime }\right)$ feltételes sűrűségfüggvény maximumát a tényleges adatainknál veszi fel. Ez egy maximum likelihood becslési feladat lenne, ahol a feltételes sűrűségfüggvény a likelihood függvény. Azonban mivel az Y teljes adathalmazt nem ismerjük, a $P\left(Y|h^{\prime }\right)$ likelihood függvény maximumát biztosító hipotézist sem tudjuk meghatározni. Meghatározhatjuk azonban a $P\left(Y|h^{\prime }\right)$ likelihood függvény várhatóértékét $E\left\{P\left(Y|h^{\prime }\right)\right\}$ -t, és a várhatóérték alapján fogalmazhatunk meg egy maximum likelihood becslési feladatot. Az $E\left\{P\left(Y|h^{\prime }\right)\right\}$ meghatározásához az Y teljes adatok ismeretére nincs szükségünk, elegendő, ha Y eloszlása ismert, ami azt jelenti, hogy a hiányzó adatok ismerete nélkül, de a hiányzó adatok eloszásának ismeretében a várhatóérték likelihood függvény megfogalmazható. A várhatóértéket az Y-t generáló valószínűségi eloszlás felett képezzük, amelyet viszont a $\theta$ paramétervektor határoz meg. Az aktuális hipoztézis (az ehhez tartozó aktuális paraméterbecslés) alapján meghatározható a várhatóérték likelihood függvény, $E\left\{P\left(Y|h^{\prime }\right)\right\}$ , majd az ennek maximumát biztosító új hipotézis. (Legtöbbször praktikus okokból e helyett az $E\left\{\text{log}\left(P\left(Y|h^{\prime }\right)\right)\right\}$ -t maximáljuk, de ez a lényegen nem változtat, a két kifejezés maximuma egybeesik.)

Definiáljunk egy $Q\left(h^{\prime }|h\right)$ függvényt, ami $E\left\{\text{log}\left(P\left(Y|h^{\prime }\right)\right)\right\}$ -t mint h’ függvényét adja meg azzal a feltételezéssel, hogy az aktuális hipotézisünk, h szerint a paramétervektor értéke $\theta$ , valamint, hogy az adataink megfigyelhető része X:

$Q\left(h^{\prime }|h\right)=E\left\{P\left(Y|h^{\prime }\right)|h,X\right\}$ (13.14)

Ezekután az EM-algoritmus az alábbi két lépést ismétli, amíg a konvergencia be nem következik:

E-lépés: várhatóérték-képzés (expectation). Meghatározzuk a $Q\left(h^{\prime }|h\right)$ függvényt a jelenlegi h hipotézis ( $\theta$ paramétervektor) és a megfigyelt X adathalmaz alapján, az Y felett nyert valószínűségi eloszlás várható értékét képezve:

$Q\left(h^{\prime }|h\right)\leftarrow E\left\{P\left(Y|h^{\prime }\right)|h,X\right\}$ (13.15)

M-lépés: maximumkeresés (maximization). A h hipotézist javítjuk, megkeressük azt a h’-t (a paramétervektor ${\theta }^{\prime }$ új becslését), amely maximálja az E-lépés során nyert Q-t:

$h\leftarrow \underset{h\text{'}}{\text{argmax}}Q\left(h^{\prime }|h\right)$ (13.16)

Az M-lépés után az új paramétervektorral újra meghatározzuk az Y feletti eloszlás (amely változott a paramétervektor változása miatt) várható értékét, majd újra maximáljuk a Q-t stb.

Az általános bemutatás után szemléltessük az EM-algoritmust a fent már említett, legismertebb alkalmazási példájával, a kevert eloszlások paraméterbecslésével. Tegyük fel, hogy L darab többdimenziós mintánk van X={x₁, x₂, …, x_L}, amelyek K különböző eloszlásból származnak a következő módon. Először valamilyen valószínűséggel eldől, hogy melyik eloszlásból jön a következő minta (bár mi ezeket a valószínűségeket nem ismerjük). Azt az eseményt, hogy a k-dik eloszlásból származik a következő minta, jelölje ω_k, ennek a priori bekövetkezési valószínűsége P(ω_k). Ezt követően az adott, k-dik eloszlásra jellemző $p\left(x|{\omega }_k,{\theta }_k\right)$ valószínűségeloszlás alapján előáll az adott érték. A ${\theta }_k$ k=1, 2, …, K ismeretlen paramétervektorok együtt adják az összetett eloszlás $\theta$ ismeretlen eredő paramétervektorát, amely meghatározza a minták eredő eloszlását.

Egy adott − bár ismeretlen − paramétersor (paramétervektor) esetén a kevert eloszlást a következő összefüggés írja le:

$p\left(x|\theta \right)={\displaystyle \sum _{k=1}^{K}p\left(x|{\omega }_k,{\theta }_k\right)}P\left({\omega }_k\right)$ (13.17)

A rögzített, de ismeretlen $\theta$ mellett a mért adathalmazunk valószínűségét a következő összefüggés, az úgynevezett likelihood függvény adja meg, ami nem más, mint az egyes értékek valószínűségeinek szorzata:

$p\left(X|\theta \right)={\displaystyle \prod _{j=1}^{L}p\left(x_j|\theta \right)}$ (13.18)

A $\theta$ ismeretlen paramétervektor maximum likelihood (ML) becsléseaz a $\widehat{\theta }$ becslő, amely a $p\left(X|\theta \right)$ -t maximálja.

$\widehat{\theta }=\underset{\theta }{\text{argmax}}p\left(X|\theta \right)$ (13.19)

Ez az összefüggés azt fejezi ki, hogy a ténylegesen mért adatsor a különböző $\theta$ -k esetén különböző valószínűségű, és érthető módon azt a $\theta$ -t választjuk ezek közül, amely mellett a mért adatsor a legvalószínűbb. Tegyük fel, hogy $p\left(X|\theta \right)$ differenciálható $\theta$ szerint, jelöljük a likelihood logaritmusát $L$ -lel, a logaritmikus likelihood függvény ${\theta }_i$ szerinti gradiensét pedig ${\nabla }_{{\theta }_i}\left(L\right)$ -lel, ekkor:

$L={\displaystyle \sum _{j=1}^L\text{log}\left(p\left(x_j|\theta \right)\right)}$ (13.20)

és

${\nabla }_{{\theta }_i}\left(L\right)={\displaystyle \sum _{j=1}^L\frac{1}{p\left(x_j|\theta \right)}}{\nabla }_{{\theta }_i}\left[{\displaystyle \sum _{k=1}^{K}p\left(x_j|{\omega }_k,{\theta }_k\right)P\left({\omega }_k\right)}\right]$ (13.21)

Tegyük fel, hogy az összes $\theta$ és ${\theta }_k$ paramétervektorok elemei függetlenek egymástól, amennyiben $i\ne k$ . Mivel $P({\omega }_i)$ nem függ ${\theta }_i$ -től:

$\begin{array}{l}{\nabla }_{{\theta }_i}\left(L\right)={\displaystyle \sum _{j=1}^L\frac{1}{p\left(x_j|\theta \right)}}{\nabla }_{{\theta }_i}p\left(x_j|{\omega }_i,{\theta }_i\right)P\left({\omega }_i\right)=\\ {\displaystyle \sum _{j=1}^L\frac{p\left(x_j|{\omega }_i,{\theta }_i\right)P\left({\omega }_i\right)}{p\left(x_j|\theta \right)}}\frac{\text{1}}{p\left(x_j|{\omega }_i,{\theta }_i\right)}{\nabla }_{\theta {}_i}p\left(x_j|{\omega }_i,{\theta }_i\right)\end{array}$ (13.22)

A kifejezést bővítettük $\frac{p\left(x_j|{\omega }_i,{\theta }_i\right)}{p\left(x_j|{\omega }_i,{\theta }_i\right)}=1$ - gyel. A Bayes tétel alapján a szumma utáni első tényező nem más, mint az ${\omega }_i$ esemény a posteriori valószínűsége. Tehát

$P\left({\omega }_i|x_j,{\theta }_i\right)=\frac{p\left(x_j|{\omega }_i,{\theta }_i\right)P\left({\omega }_i\right)}{p\left(x_j|\theta \right)}$ (13.23)

Ez annak valószínűségét adja meg, hogy a megfigyelt x_j minta az i-edik eloszlásból származik adott $\theta$ (és ezen belül ${\theta }_i$ ) mellett. Ezek után az (13.22) egyenlet a következő alakra hozható:

${\nabla }_{{\theta }_i}\left(L\right)={\displaystyle \sum _{j=1}^{L}P\left({\omega }_i|x_j,{\theta }_i\right)}{\nabla }_{{\theta }_i}\text{log}\left(p\left(x_j|{\omega }_i,{\theta }_i\right)\right)$ $$ (13.24)

Mikor elérjük a maximumot, akkor a gradiens 0-vá válik, tehát az ML becslés ki kell elégítse a következő összefüggéseket:

${\displaystyle \sum _{j=1}^{L}P\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)}{\nabla }_{{\theta }_i}\text{log}\left(p\left(x_j|{\omega }_i,{\widehat{\theta }}_i\right)\right)=0i=\mathrm{1,...,}K$ (13.25)

Tehát ezen egyenletrendszer megoldásai között kell keressük az ML becslőt.

Ha tudnánk, hogy melyik minta melyik eloszlásból származik, akkor az ismeretlen paramétereket (pl. az eloszlások várható értékét, szórását stb.) viszonylag egyszerűen tudnánk becsülni. Ezek viszont a nem megfigyelt, nem megfigyelhető változók. Így viszont az alábbi nemlineáris egyenletrendszert kéne megoldanunk (amennyiben sem az a priori eloszlások, sem az egyes ${\theta }_i$ eloszlásparaméterek nem ismertek):

$$ $\widehat{P}\left({\omega }_i\right)=\frac{1}{L}{\displaystyle \sum _{j=1}^L\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)}$ (13.26)

$\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)=\frac{p\left(x_j|{\omega }_i,{\widehat{\theta }}_i\right)\widehat{P}\left({\omega }_i\right)}{{\displaystyle \sum _{k=1}^{K}p\left(x_j{|{\omega }_k,\theta }_k\right)\widehat{P}\left({\omega }_k\right)}}$ (13.27)

${\displaystyle \sum _{j=1}^L\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)}{\nabla }_{{\theta }_i}\text{log}\left(p\left(x_j|{\omega }_i,{\widehat{\theta }}_i\right)\right)=0i=\mathrm{1,...,}K$ (13.28)

Foglaljuk össze, hogy mit ismerünk, és mit kívánunk megbecsülni ebben az egyenletrendszerben. Két dolgot ismerünk:

- az X={x₁, x₂, …, x_L} mintahalmazt,

az egyes eloszlások struktúráját (legtöbbször pl. azt tudjuk vagy fetételezzük, hogy normál eloszlásról van szó). Ezt fejezi ki az egyenletrendszerben a $p\left(x_j|{\omega }_i,{\widehat{\theta }}_i\right)$ szerepeltetése, amely ${\theta }_i$ paraméterekkel jellemezhető.

Nem ismerjük:

- az egyes eloszlások a priori valószínűségeit (ezért használunk az összefüggésekben $\widehat{P}\left({\omega }_i\right)$ -t, ami a valószínűség becsült voltára utal),

- az egyes eloszlások paramétereit (erre utal a $\widehat{\theta }$ használata),

nem tudjuk, hogy melyik minta melyik eloszlásból származik − ez tulajdonképpen az az adat, amit nem tudunk megfigyelni −, és természetesen nem ismerjük az a posteriori valószínűsége- ket sem ( $\widehat{P}\text{(}{\omega }_i|x_j,{\widehat{\theta }}_i\text{)}$ ).

Az egyenletrendszer gyakran nem oldható meg vagy a megoldás nagyon nagy számításigényű. Ennek megfelelően az EM-algoritmus iteratív megoldást kínál, a következő két lépés ismételt alkalmazásával. Első lépésben a (13.26) és a (13.27) összefüggések alapján, a pillanatnyilag használt $\theta$ segítségével kiszámítjuk az a posteriori valószínűség eloszlást, és a várható értékét képezzük − $E\left\{\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)\right\}$ (expectation-lépés). Ezek után az így kiszámított várható értékeket írjuk be a (13.28) logaritmikus likelihood egyenletbe, és megkeressük azt az új $\theta$ értéket, amely maximálja a logaritmikus likelihood-ot (maximization-lépés). Ezt a két lépést ciklikusan valamilyen kezdeti értékekből kiindulva, valamilyen leállási feltétel eléréséig folytatjuk. A kezdetiérték választás, a konvergencia, a leállási feltételek kérdéseinek tárgyalása meghaladják a jelen könyv kereteit.

A legtöbbször vizsgált eset, amikor az egyes részeloszlások Gauss-jellegűek. Ilyenkor a következő − nagyon logikus − eredményre jutunk az ismeretlen átlagok (és szórások) becslésénél. Az eredményeket − a jobb áttekinthetőség kedvéért − skalár részeloszlások esetére közöljük, amikor csak az eloszlások várható értéke ( $m_i$ ) ismeretlen, a szórásuk azonos és ismert. (Természetesen az is ismeretlen, hogy a keverékből melyik adat melyik eloszláshoz tartozik.)

$E\left\{\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)\right\}=\frac{e^{\frac{-{\left(x_j-m_j\right)}^2}{2{\sigma }^2}}}{{\displaystyle \sum _r{e}^{\frac{-{\left(x_j-m_r\right)}^2}{2{\sigma }^2}}}}$ (13.29)

$m_i=\frac{{\displaystyle \sum _{j}E\left\{\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)\right\}{x}_j}}{{\displaystyle \sum _{j}E\left\{\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)\right\}}}$ (13.30)

Látható, hogy ha tudnánk, hogy melyik minta melyik részeloszlásból származik, akkor visszakapnánk a várható érték becslésére a jól ismert átlagértéket – azon értékek átlagát, amelyek az adott részeloszláshoz tartoznak. (Az $E\left\{\widehat{P}\left({\omega }_i|x_j,{\widehat{\theta }}_i\right)\right\}$ várhatóértékek 1, illetve 0 értéket vennének fel, attól függően, hogy az adott mérés a kiválasztott eloszláshoz tartozik-e vagy sem.)