# How to calculate log likelihood for gaussian mixture model

I'm trying to check the calculation of the log likelihood of a 2 component Gaussian Mixture Model using optim, but I get the wrong answer (it should return mu, sigma, alpha actual).

The log likelihood is calculated as:

$$log(p(y_i|\alpha,\mu,\sigma))=\sum_{i=1}^Nlog(\sum_{k=1}^2\alpha_k*N(y_i|\mu_k,\sigma_k))$$

using the method shown on this blog. Thanks for any help.

#generate response
mu_actual <- c(-5,3)
sigma_actual <- c(2,3)
alpha_actual <- c(0.3,0.7)

comp <- rbinom(n=1e3,size=1,prob=alpha_actual[1])
Y <- rnorm(1e3, mean=mu_actual[comp+1], sd=sigma_actual[comp+1])

#likelihood calculation
fn_likelihood_optim <- function(theta,Y){
mu_1 <- theta[1]
mu_2 <- theta[2]
sigma_1 <- theta[3]
sigma_2 <- theta[4]
alpha_1 <- theta[5]
alpha_2 <- 1-theta[5]

like_1 <- log(alpha_1)+dnorm(Y,mu_1,sigma_1,log=TRUE)
like_2 <- log(alpha_2)+dnorm(Y,mu_2,sigma_2,log=TRUE)
log_like <- sum(log(exp(like_1)+exp(like_2)))

if(is.infinite(log_like) & log_like>0) log_like <- 1e6
if(is.infinite(log_like) & log_like<0) log_like <- 1e-6

return(-log_like)
}

#check likelihood with optim - should return theta actual
z <- optim(par = c(0,5,1,2,0.5),
fn = fn_likelihood_optim,
method = "L-BFGS-B",
lower=c(-10,-10,0.001,0.001,0),
upper=c(10,10,10,10,0.99),
Y=Y)

z$par  ## 1 Answer You are getting infinity if alpha_1 every touches zero, so set the alpha_1 lower boundary to 0.01 or something not zero. I think there is a typo here if(is.infinite(log_like) & log_like>0) log_like <- 1e6, should be 1e-6. set.seed(111) comp <- rbinom(n=1e3,size=1,prob=alpha_actual[1]) Y <- rnorm(1e3, mean=mu_actual[comp+1], sd=sigma_actual[comp+1]) fn_likelihood_optim <- function(theta,Y){ mu_1 <- theta[1] mu_2 <- theta[2] sigma_1 <- theta[3] sigma_2 <- theta[4] alpha_1 <- theta[5] alpha_2 <- 1-theta[5] like_1 <- log(alpha_1)+dnorm(Y,mu_1,sigma_1,log=TRUE) like_2 <- log(alpha_2)+dnorm(Y,mu_2,sigma_2,log=TRUE) log_like <- sum(log(exp(like_1)+exp(like_2))) if(is.infinite(log_like) & log_like>0) log_like <- 1e-6 if(is.infinite(log_like) & log_like<0) log_like <- 1e-6 return(-log_like) } #check likelihood with optim - should return theta actual z <- optim(par = c(0,5,1,2,0.5), fn = fn_likelihood_optim, method = "L-BFGS-B", lower=c(-10,-10,0.001,0.001,0.01), upper=c(10,10,10,10,0.99), Y=Y) z$par
[1]  3.0831575 -5.0493771  2.7413426  1.9876056  0.2740447

• thanks @StupidWolf I think the check for infinite likelihoods can be dropped if alpha_1 is constrained to be greater than 0
– Zeus
Commented May 14, 2020 at 4:25