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 (which helps to minimize by expressing derivatives).
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
# (this is done here indirectly/computationally with optim,
# but could be done analytically be expressing te derivative and solving)
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