EM-algoritmi

Wikipedia
Loikkaa: valikkoon, hakuun

EM-algoritmi (EM on lyhenne sanoista Expectation-Maximization eli odotusarvon 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.

Vanhan uskollisen purkautumisiin liittyvän aineiston EM-klusterointi. Jokin aloitusmalli sovitetaan havaittuun aineistoon. (Akseleiden erilaisten mitta-asteikoiden vuoksi jakauma näyttää kahdelta hyvin litteältä ja leveältä soikiolta.) Ensimmäiset iteraatiot muuttavat mallia huomattavasti, minkä jälkeen malli konvergoi kohti geysirin purkausten tyypillisimpiä arvoyhdistelmiä. Visualisointi tehty ELKI:llä.


Kuvaus[muokkaa | muokkaa wikitekstiä]

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. [1][2]

Ominaisuuksia[muokkaa | muokkaa wikitekstiä]

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

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). [4][5][6][7][8][9][10]


Esimerkkejä[muokkaa | muokkaa wikitekstiä]

Gaussinen sekoitus[muokkaa | muokkaa wikitekstiä]

Animaatio, joka havainnollistaa kaksikomponenttisen sekoitusjakauman sovittamista EM-algoritmilla Old Faithful aineistoon. Algoritmin askeleet satunnaisesta aloituksesta konvergenssiin.

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.[2]

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[muokkaa | muokkaa wikitekstiä]

Oletetaan, että meillä on parametrien estimaatit θ(t). Tällöin Zi: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[muokkaa | muokkaa wikitekstiä]

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[muokkaa | muokkaa wikitekstiä]

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[muokkaa | muokkaa wikitekstiä]

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


Historiaa[muokkaa | muokkaa wikitekstiä]

EM-algoritmin historia jaetaan usein kirjoittajien Dempster, Laird ja Rubin vuonna 1977 ilmestynyttä artikkelia[11] 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. [1]

Ennen DLR-artikkelia[muokkaa | muokkaa wikitekstiä]

Kirjallisuudessa ensimmäinen maininta liittyen EM-tyyppiseen algoritmiin esiintyy Newcombin vuoden 1886 artikkelissa [12] kahden yksiulotteisen normaalijakauman sekoituksesta.

Vuonna 1960 Buck [13] 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.[14][15][16]

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

DLR-artikkelin jälkeen[muokkaa | muokkaa wikitekstiä]

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

Lähteet[muokkaa | muokkaa wikitekstiä]

  1. a b c (1997) The EM algorithm and extensions. New York: Wiley. ISBN 0-471-12358-7. 
  2. a b c (2001) "8.5 The EM algorithm", The Elements of Statistical Learning. New York: Springer, 236–243. ISBN 0-387-95284-5. 
  3. Navidi, William (1997). "A Graphical Illustration of the EM Algorithm". The American Statistician 51 (1): 29-31. 
  4. Sundberg, Rolf (1971). "Maximum likelihood theory and applications for distributions generated when observing a function of an exponential family variable". dissertation. 
  5. Sundberg, Rolf (1976). "An iterative method for solution of the likelihood equations for incomplete data from exponential families". Communications in Statistics – Simulation and Computation 5 (1): 55–64. doi:10.1080/03610917608812007. 
  6. Martin-Löf, Anders (1963). "Utvärdering av livslängder i subnanosekundsområdet (Evaluation of sub-nanosecond lifetimes) ("Sundbergin kaava")". 
  7. Martin-Löf, Per. 1966. Statistics from the point of view of statistical mechanics. Lecture notes, Mathematical Institute, Aarhus University. ("Sundberg formula" credited to Anders Martin-Löf).
  8. Martin-Löf, Per. 1970. Statistika Modeller (Statistical Models): Anteckningar från seminarier läsåret 1969–1970 (Notes from seminars in the academic year 1969-1970), with the assistance of Rolf Sundberg. Stockholm University. ("Sundberg formula")
  9. Martin-Löf, P. The notion of redundancy and its use as a quantitative measure of the deviation between a statistical hypothesis and a set of observational data. With a discussion by F. Abildgård, A. P. Dempster, D. Basu, D. R. Cox, A. W. F. Edwards, D. A. Sprott, G. A. Barnard, O. Barndorff-Nielsen, J. D. Kalbfleisch and G. Rasch and a reply by the author. Proceedings of Conference on Foundational Questions in Statistical Inference (Aarhus, 1973), pp. 1–42. Memoirs, No. 1, Dept. Theoret. Statist., Inst. Math., Univ. Aarhus, Aarhus, 1974.
  10. Martin-Löf, Per (1974). "The notion of redundancy and its use as a quantitative measure of the discrepancy between a statistical hypothesis and a set of observational data". Scandinavian Journal of Statistics 1 (1): 3–18. 
  11. "Maximum Likelihood from Incomplete Data via the EM Algorithm" (1977). Journal of the Royal Statistical Society, Series B 39 (1): 1–38. 
  12. Newcomb, S. (1886). "A generalized theory of the combination of observations so as to obtain the best results". American Journal of Mathematics 8: 343-366. 
  13. Buck, S.F (1960). "A method of estimation of missing values in multivariate data suitable for use with an electronic computer". Journal of the Royal Statistical Society B 22: 302-306. 
  14. "Statistical inference for probabilistic functions of finite Markov chains" (1966). Annals of Mathematical Statistics 37: 1554-1563. 
  15. "An inequality with applications to statistical estimation for probabilistic functions of Markov processes and to a model for ecology" (1967). Bulletin of the American Mathematical Society 73: 360-363. 
  16. "A Maximization technique occuring in the statistical analysis of probabilistic functions of Markov chains" (1970). Annals of Mathematical Statistics 41: 164-171. 
  17. "A missing information principle: theory and applications" (1972). Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability 1: 697-715. Berkeley, California: University of California Press.