R/logL_mcfa.R

In EMMIXmfa: Mixture Models with Component-Wise Factor Analyzers

loglike.mcfa <- function(Y, g, q, pivec, A, xi, omega, D, ...) {

p <- ncol(Y)
if (is.null(p))
  p <- 1
n <- nrow(Y)

Fji <- array(NA, c(n, g))
inv_D <- diag(1 / diag(D))
#  inv_D %*% A
inv_D_A <- sweep(A, 1, diag(inv_D), '*')
for (i in 1 : g) {

  inv_O <- try(chol.inv(chol.inv(omega[,, i]) + 
                        t(inv_D_A) %*% A), silent = TRUE) 
  if (any(class(inv_O) %in% "try-error"))
    return(loglike <-
        paste('Ill-conditioned or singular Sigma[,', i, ']'))

  inv_S <- try(inv_D - inv_D_A %*% inv_O %*% t(inv_D_A), silent = TRUE)
  if (any(class(inv_S) %in% "try-error")) {
    return(loglike <- paste('Ill-conditioned or singular Sigma_',
                      i, sep = ''))
  }

  # |D + A \Om A^T| using matrix determinant lemma
  # |D| + |\Om| + |\Om^{-1} + A^T D^{-1} A| 
  logdetD <- log(det(as.matrix(omega[,, i]))) +
                sum(log(diag(D))) -
                #+ log(det(chol.inv(omega[,, i]) + t(A) %*% inv_D %*% A))
                log(det(inv_O))
                
  mahal_dist <- mahalanobis(Y, t(A %*% xi[, i, drop = FALSE]), inv_S, TRUE)

  Fji[, i] <- -0.5 * mahal_dist - (p/2) * log(2 * pi) - 0.5 * logdetD
}

Fji <- sweep (Fji, 2, log(pivec), '+')
Fjmax <- apply (Fji, 1, max)
Fji <- sweep (Fji, 1, Fjmax, '-')
loglike <- sum (Fjmax, log(rowSums(exp(Fji))))

return(loglike)
}