###Short overview about Expectation maximization :
- Marginal likelihood
Expectation maximization contrasts with 'regular' likelihood maximization by refering to the maximization of a marginal likelihood.
$$\underbrace{p(X\vert \theta)}_{\substack{\text{marginal likelihood}\\\text{ $\mathcal{L}(\theta \vert X)$}}} = \int_z \underbrace{p(X, z \vert \theta)}_{\substack{\text{likelihood}\\\text{ $\mathcal{L}(\theta \vert X,\underset{\uparrow \\ \substack{\llap{\text{This $z$ is }\rlap{\text{missing data}}}}}{z})$}}} \text{d}z = \int_z p(X \vert \theta,z) p(z\vert X,\theta) \text{d}z $$
So this relates to an integral over some likelihood with an additional parameter $z$ (e.g. missing data).
- EM algorithm
In the EM algorithm this integral is not maximized directly:
$$\hat\theta = \underset{\theta}{\text{arg max}} \left( \int_z p(X \vert \theta,z) p(z\vert X,\theta) \text{d}z \right)$$
but instead it is computed in an iterative way:
$$\hat\theta_{k+1} = \underset{\theta}{\text{arg max}} \left( \int_z p(X \vert \theta,z) p(z\vert X,\hat\theta_k) \text{d}z \right)$$
This is done by picking an initial $\theta_1$ and updating repetitively. Note that now the optimization keeps the term $p(z\vert X,\theta)$ fixed.
###Not all problems are like that.
So this marginal likelihood is only the case for problems with unobserved data. For instance in finding a Gaussian Mixture (example on wikipedia) one may consider an observed variable $z$ that refers the class (which component in the mixture) that a measurement belongs to.
There are many problems that do not consider a marginal likelihood and evaluate parameters directly in order to optimize some likelihood function (or some other cost function but not a marginalized/expectation one).
###Code example
The use of marginal likelihood is not always about explicitly missing data.
In the example below the (example on wikipedia) is worked out numerically in R.
In this case you do not have explicitly a case with missing data (the likelihood can be defined directly and does not need to be a marginal likelihood integrating over missing data), but one has a mixture of two multivariate Gaussian distributions. The problem with that is that one can not do as normally and compute the sum of the logarithm of the terms (which has computational advantages), because not the terms are not a product, but they involve a sum as well. Although, it would be possible to compute those logarithms of sums using an approximation (this is not done in the code below, instead the parameters are created in order to be advantageous and do not generate infinite values).
library(MASS)
# data example
set.seed(1)
x1 <- mvrnorm(100, mu=c(1,1), Sigma=diag(1,2))
x2 <- mvrnorm(30, mu=c(3,3), Sigma=diag(sqrt(0.5),2))
x <- rbind(x1,x2)
col <- c(rep(1,100),rep(2,30))
plot(x,col=col)
# Likelihood without integrating over z
Lsimple <- function(par,X=x) {
tau <- par[1]
mu1 <- c(par[2],par[3])
mu2 <- c(par[4],par[5])
sigma_1 <- par[6]
sigma_2 <- par[7]
likterms <- tau*dnorm(X[,1],mean = mu1[1], sd = sigma_1)*
dnorm(X[,2],mean = mu1[2], sd = sigma_1)+
(1-tau)*dnorm(X[,1],mean = mu2[1], sd = sigma_2)*
dnorm(X[,2],mean = mu2[2], sd = sigma_2)
logLik = sum(log(likterms))
-logLik
}
# Marginal likelihood integrating over z
LEM <- function(par,X=x,oldp=oldpar) {
tau <- par[1]
mu1 <- c(par[2],par[3])
mu2 <- c(par[4],par[5])
sigma_1 <- par[6]
sigma_2 <- par[7]
oldtau <- oldp[1]
oldmu1 <- c(oldp[2],oldp[3])
oldmu2 <- c(oldp[4],oldp[5])
oldsigma_1 <- oldp[6]
oldsigma_2 <- oldp[7]
f1 <- oldtau*dnorm(X[,1],mean = oldmu1[1], sd = oldsigma_1)*
dnorm(X[,2],mean = oldmu1[2], sd = oldsigma_1)
f2 <- (1-oldtau)*dnorm(X[,1],mean = oldmu2[1], sd = oldsigma_2)*
dnorm(X[,2],mean = oldmu2[2], sd = oldsigma_2)
pclass <- f1/(f1+f2)
### note that now the terms are a product and can be replaced by a sum of the log
#likterms <- tau*dnorm(X[,1],mean = mu1[1], sd = sigma_1)*
# dnorm(X[,2],mean = mu1[2], sd = sigma_1)*(pclass)*
# (1-tau)*dnorm(X[,1],mean = mu2[1], sd = sigma_2)*
# dnorm(X[,2],mean = mu2[2], sd = sigma_2)*(1-pclass)
loglikterms <- (log(tau)+dnorm(X[,1],mean = mu1[1], sd = sigma_1, log = TRUE)+
dnorm(X[,2],mean = mu1[2], sd = sigma_1, log = TRUE))*(pclass)+
(log(1-tau)+dnorm(X[,1],mean = mu2[1], sd = sigma_2, log = TRUE)+
dnorm(X[,2],mean = mu2[2], sd = sigma_2, log = TRUE))*(1-pclass)
logLik = sum(loglikterms)
-logLik
}
# solving with direct likelihood
par <- c(0.5,1,1,3,3,1,0.5)
p1 <- optim(par, Lsimple,
method="L-BFGS-B",
lower = c(0.1,0,0,0,0,0.1,0.1),
upper = c(0.9,5,5,5,5,3,3),
control = list(trace=3, maxit=10^3))
p1
# solving with direct LEM
oldpar <- c(0.5,1,1,3,3,1,0.5)
for (i in 1:100) {
p2 <- optim(oldpar, LEM,
method="L-BFGS-B",
lower = c(0.1,0,0,0,0,0.1,0.1),
upper = c(0.9,5,5,5,5,3,3),
control = list(trace=1, maxit=10^3))
oldpar <- p2$par
print(i)
}
p2
# the result is the same:
p$par
p2$par