Machine Learning: Advanced Probabilistic Methods (2015), P.
Marttinen
[The example is modi…ed from lecture slides for the course "Laskentaintensiiviset
tilastolliset menetelmät (Computational statistics)" by Petri Koistinen]
Suppose that we have N independent observations x = (x1 ; : : : ; xN ) from a
two-component mixture of univariate Gaussian distributions
1 1
p(xn j ) = N (xn j0; 1) + N (xn j ; 1): (1)
2 2
This means that with probability 1/2 the observation xn is generated from the
…rst component N (xn j0; 1); and with probability 1/2 from the second component
N (xn j ; 1). The model (1) has one unknown parameter, , representing the
mean of the second component, and we would like to estimate it using maximum
likelihood
b = arg max flog p(xj )g .
We do this by the EM-algorithm (although direct numerical optimization would
also be straightforward for this simple model).
First we formulate the model using the latent variable representation,
and introduce variables z =(z1 ; : : : ; zN ) which explicitly specify the component
responsible for generating observation xn . In detail:
(1; 0)T ; (xn is from N (xn j0; 1))
zn = (zn1 ; zn2 )T = :
(0; 1)T ; (xn is from N (xn j ; 1))
When we de…ne the distributions for the latent variable model as follows
p(zn1 = 1) = p(zn2 = 1) = 0:5
and
N (xn j0; 1); if zn1 = 1
p(xn jzn ; ) =
N (xn j ; 1); if zn2 = 1
it is easy to see that the marginal distribution of xn obtained by summing over
the latent variables X
p(xn j ) = p(xn jzn ; )p(zn )
zn
is equal to the original distribution (1).
In the EM-algorithm we will maximize the expectation of the log-likelihood
of the complete data (x; z):
(N ) N
Y X
log p(x; zj ) = log p(xn ; zn j ) = log p(xn ; zn j )
n=1 n=1
N
X
= log [0:5 N (xn j0; 1)zn1 N (xn j ; 1)zn2 ]
n=1
N
X
= fzn1 log [N (xn j0; 1)] + zn2 log [N (xn j ; 1)]g + const (2)
n=1
1
E-step 10 : Compute the posterior distribution of the latent variables, given the
current estimate 0 of :
p(zn1 = 1jxn ; 0 ) / p(zn1 = 1)p(xn jzn ; 0)
= 0:5 N (xn j0; 1) (3)
p(zn2 = 1jxn ; 0 ) / p(zn2 = 1)p(xn jzn ; 0)
= 0:5 N (xn j 0 ; 1) (4)
By normalizing (3) and (4) we get
N (xn j 0 ; 1)
(zn2 ) p(zn2 = 1jxn ; 0) = . (5)
N (xn j0; 1) + N (xn j 0 ; 1)
E-step 20 : Evaluate the expectation of the complete data log-likelihood (2)
over the posterior distribution of the latent variables (5):
Q( ; 0) = Ezjx; 0 [log p(x; zj )]
N
X
= fE[zn1 ] log [N (xn j0; 1)] + E[zn2 ] log [N (xn j ; 1)]g
n=1
N
X
= f[1 (zn2 )] log [N (xn j0; 1)] + (zn2 ) log [N (xn j ; 1)]g : (6)
n=1
Note that in (6) we’ve discarded the term not dependent on in equation (2).
As a matter of fact, the …rst term in each sum could also be discarded, but we
retain it here for clarity.
M-step: Maximize Q( ; 0 ) with respect to . To di¤erentiate Q( ; 0 ); we …rst
note the following result, which can be veri…ed by straightforward computation
d
N (xn j ; 1) = N (xn j ; 1)(xn ):
d
With this result at hand, we can write
N
d d X
Q( ; 0) = f[1 (zn2 )] log [N (xn j0; 1)] + (zn2 ) log [N (xn j ; 1)]g
d d n=1
N
X N
X
(zn2 )
= N (xn j ; 1)(xn )= (zn2 )(xn ):
n=1
N (xn j ; 1) n=1
d
Setting d Q( ; 0) = 0; we get
PN
n=1 (zn2 )xn
= PN
n=1 (zn2 )
XN
1
= (zn2 )xn ; (7)
N2 n=1
2
PN
where we have de…ned N2 = n=1 (zn2 ); which can be interpreted as the
e¤ective number of observations assigned to component 2. Equation (7) has
an intuitive interpretation: the mean of component (cluster) 2 is obtained as
a weighted average of all points in the data set, and the weight of data point
xn is equal to the posterior probability (or responsibility) (zn2 ) that the 2nd
component was responsible for generating xn .
Code to run the EM-algorithm: simple_em.m