2
$\begingroup$

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
$\endgroup$

1 Answer 1

1
$\begingroup$

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
$\endgroup$
1
  • 1
    $\begingroup$ thanks @StupidWolf I think the check for infinite likelihoods can be dropped if alpha_1 is constrained to be greater than 0 $\endgroup$
    – Zeus
    May 14, 2020 at 4:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.