<!-- language-all: lang-r -->

###Short overview about Expectation maximization :

 - **Marginal likelihood**

 Expectation maximization contrasts with 'regular' likelihood maximization by refering to the maximization of a [marginal likelihood][1].

 $$\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 min}} \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 min}} \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][2]) 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][2]) 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][3] (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


  [1]: https://en.wikipedia.org/wiki/Marginal_likelihood
  [2]: https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm#Gaussian_mixture
  [3]: https://math.stackexchange.com/a/734254/466748