EM-algoritmi

**EM-algoritmi** (EM on  eli *[[odotusarvo]]n maksimointi*) on tilastotieteessä käytetty iteratiivinen menetelmä suurimman uskottavuuden estimaattien löytämiseksi tilastollisten mallien parametreille tilanteessa, jossa osa tiedosta puuttuu.
Puuttuva tieto voi olla esimerkiksi piilevä luokkamuuttuja, josta ei saatu lainkaan havaintoja.

## Kuvaus

Olkoon \theta = \big( \theta_1,\theta_2,... ,\theta_k \big) aineiston Y jakaumaan liittyvien tuntemattomien parametrien muodostama vektori. Täydelliselle aineistolle uskottavuusfunktio voidaan kirjoittaa muodossa
:  L(Y|\theta) = \prod_{i=1}^{n}f(y_i,\theta) .

Hyvin usein osa oleellisista tiedosta jää kuitenkin havaitsematta. Täydellinen data voidaan jakaa kahteen osaan: havaittuun aineistoon X ja puuttuvaan aineistoon Z. Tällöin täydellisen aineiston uskottavuus saadaan kirjoitettua muotoon
:L(X,Z|\theta) = L(X|\theta)L(Z|X,\theta) .
Ottamalla logaritmi puolittain lauseke saadaan muotoon 
:\log L(X|\theta) = \log L(X,Z|\theta) - \log L(Z|X,\theta) .

Oletetaan mallin parametreille tunnettu arvo \theta = \theta^{(t)} iteraatiolla t. Tällöin edellä esitetyn lausekkeen odotusarvo puuttuvien havaintojen suhteen on 
:\operatorname{E}_{Z|X,\theta=\theta^{(t)}}\left[ \log L (X|\theta)  \right]
=\operatorname{E}_{Z|X,\theta=\theta^{(t)}}\left[ \log L (X,Z|\theta) -  \log L (Z|X,\theta)  \right] .

Merkitään nyt täydellisen aineiston logaritmista uskottavuutta seuraavasti:
:\begin{align}
 Q(\theta|\theta^{(t)}) &= \operatorname{E}_{Z|X,\theta=\theta^{(t)}}\left[ \log L (X,Z|\theta)\right]\\
&= \int \log L (X,Z|\theta) L (Z|X ,\theta^{(t)})\, dZ \end{align}

Algoritmissa toistetaan vuorotellen kahta askelta: 
:**E-askel**: Johda termin Q(\theta|\theta^{(t)}) lauseke.
:**M-askel**: Etsi parametrille \boldsymbol\theta sellainen arvo \boldsymbol\theta^{(t+1)}, että uskottavuus maksimoituu.

Aluksi tuntemattomille parametreille asetetaan alkuarvot \theta^{(0)}. Ensimmäinen iteraatio aloitetaan siis laskemalla Q(\theta|\theta^{(0)}).

Varsinainen iteratiivinen algoritmi joudutaan johtamaan erikseen kullekin tilanteelle. 

## Ominaisuuksia

Käytettäessä EM-algoritmia uskottavuusfunktion arvo kasvaa jokaisella iteraatiolla ja parametrin estimaatti lähestyy monotonisesti suurimman uskottavuuden estimaattia.

EM-algoritmi on hyödyllinen uskottavuuden tullessa eksponenttisesta perheestä:
E-askel sievenee tyhjentävien tunnuslukujen odotusarvojen summaksi ja M-askeleessa maksimoidaan lineaarista funktiota. Tällaisessa tapauksessa voidaan usein johtaa suljettu muoto askelten päivittämiseksi Sundbergin kaavalla (Rolf Sundberg julkaisi kaavan, mutta hän hyödynsi  Per Martin-Löfin ja Anders Martin-Löfin julkaisemattomia tuloksia).

## Esimerkkejä
### Gaussinen sekoitus

Olkoon \mathbf{x} = (\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n)  n-kokoinen otos riippumattomia havaintoja kahdesta moniulotteisesta normaalijakaumasta, ulottuvuuksien määrä d>1.
Olkoot \mathbf{z} = (z_1,z_2,\ldots,z_n) latentteja muuttujia, jotka kertovat kummasta ryhmästä kyseinen havainto on peräisin.
:X_i |(Z_i = 1) \sim \mathcal{N}_d(\boldsymbol{\mu}_1,\Sigma_1) \, ja \,  X_i |(Z_i = 2) \sim \mathcal{N}_d(\boldsymbol{\mu}_2,\Sigma_2),
missä
:\operatorname{P} (Z_i = 1 ) = \tau_1 \,  ja \operatorname{P} (Z_i=2) = \tau_2 = 1-\tau_1.

Tavoite on estimoida jakaumien tuntemattomat keskiarvot ja kovarianssit, sekä jakaumien sekoittumista kuvaava arvo \tau:
:\theta = \big( \boldsymbol{\tau},\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma_1,\Sigma_2 \big),
missä uskottavuusfunktio on:
: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) ,
missä \mathbb{I} on indikaattorifunktio ja f on moniulotteisen normaalijakauman tiheysfunktio. Tämä voidaan kirjoittaa uudelleen eksponenttisen perheen muotoon:
: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\}. 
Voidaan huomata, että kullekin *i* indikaattori \mathbb{I}(z_i=j) saa arvon yksi vain yhdellä *j*, ja toisella *j* indikaattorin arvo on nolla. Sisempi summa siis supistuu yhdeksi lausekkeeksi eikä summausta tarvita.

#### E-askel

Oletetaan, että meillä on parametrien estimaatit *θ*(*t*). Tällöin *Z**i*:n ehdollinen jakauma voidaan kirjoittaa todennäköisyytenä Bayesin kaavan mukaisesti:
: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)})} .

Siten E-askel johtaa seuraavaan funktioon:
:\begin{align}Q(\theta|\theta^{(t)})
&= \operatorname{E} [\log L(\theta;\mathbf{x},\mathbf{z}) ] \\
&= \operatorname{E} [\log \prod_{i=1}^{n}L(\theta;\mathbf{x}_i,\mathbf{z}_i) ] \\
&= \operatorname{E} [\sum_{i=1}^n \log L(\theta;\mathbf{x}_i,\mathbf{z}_i) ] \\
&= \sum_{i=1}^n\operatorname{E} [\log L(\theta;\mathbf{x}_i,\mathbf{z}_i) ] \\
&= \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}

#### M-askel
Huomataan, että  \tau, (\mu_1,\Sigma_1) ja (\mu_2,\Sigma_2) voidaan kukin maksimoida toisistaan riippumatta, sillä ne ovat edellä esitetyssä lausekkeessa eri termeissä.

Tarkastellaan aluksi parametria *τ*, jolla on rajoite *τ*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}
Tämä on samaa muotoa kuin binomijakauman suurimman uskottavuuden estimaatti. Siten
:\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)}.

Tarkastellaan seuraavaksi parametrien (\mu_1,\Sigma_1) estimaatteja:
:\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}
Tämä on samaa muotoa normaalijakauman painotetun SU-estimaatin kanssa, joten
:\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)}}  ja \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)}} .
Vastaavasti
:\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)}}  ja \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)}} .

#### Lopettaminen

Lopeta iterointi, jos \log L(\theta^{t};\mathbf{x},\mathbf{z}) ja \log L(\theta^{(t-1)};\mathbf{x},\mathbf{z}) ovat riittävän lähellä toisiaan (erotus alle jonkin ennalta asetetun kynnysarvon).

#### Yleistäminen

Yllä esitetty algoritmi voidaan yleistää useampien kuin kahden monimuuttujaisen normaalijakauman sekoituksille.

## Historiaa
EM-algoritmin historia jaetaan usein kirjoittajien Dempster, Laird ja Rubin vuonna 1977 ilmestynyttä artikkelia edeltävään ja sitä seuraavaan aikaan. Kyseisessä artikkelissa annettiin runsaasti esimerkkejä algoritmin sovelluksista, ja kuvailtiin sen konvergenssiä ja muita perusominaisuuksia. Tätä artikkelia kutsutaan usein DLR-artikkeliksi. 

### Ennen DLR-artikkelia
Kirjallisuudessa ensimmäinen maininta liittyen EM-tyyppiseen algoritmiin  esiintyy Newcombin vuoden 1886 artikkelissa  kahden yksiulotteisen normaalijakauman sekoituksesta.

Vuonna 1960 Buck  esitteli p-ulotteisen populaation keskiarvovektorin ja kovarianssimatriisin estimointia tilanteessa, jossa osa aineistosta puuttui. Hän käytti regressiota ja puuttuvien havaintojen selittämistä havaitulla aineistolla. Hänen menetelmässään tarvitut regressiokertoimet ja kovarianssimatriisin kerrointen korjaukset saadaan yhdellä täydellisten havaintojen informaatiomatriisin kääntämisellä ja sopivilla matriisilaskuilla. EM-algoritmin peruselementit esiintyvät Buckin menetelmässä.

EM-algoritmin soveltamista Markov-malleille käsiteltiin sarjassa artikkeleita: Baum ja Petrie (1966), Baum ja Eagon (1967) ja Baum, Petrie, Soules ja Weiss (1970). Nämä artikkelit sisältävät helposti yleistettävissä olevia konvergenssituloksia. Näissä artikkeleissa kehitetty algoritmi toimii myös perustana nykyään käytetyille piilo-Markov-mallien EM-algoritmeille.

Vuonna 1972 Orchard ja Woodbury esittelivät täydellisen ja ei-täydellisen aineiston logaritmisten uskottavuusfunktioiden välisen suhteen.

### DLR-artikkelin jälkeen

Rajapyykkinä toimivan artikkelin jälkeen EM-algoritmia on sovellettu muun muassa neuroverkkoihin,
koneoppimisessa, psykometriikassa ja lääketieteellisessä kuvantamisessa (esimerkiksi PET-kuvauksissa).

## Lähteet