R/tau_mfa.R

In EMMIXmfa: Mixture Models with Component-Wise Factor Analyzers

tau.mfa <- function(Y, g, q, pivec, B, mu, D, sigma_type, D_type, ...) {

if (!is.matrix(Y))
	Y <- as.matrix(Y)
p <- ncol(Y)
n <- nrow(Y)
Fji <- array(NA, c(n, g))

if (sigma_type == 'common') {

	if (D_type == 'common') {
		inv_D <- diag(1 / diag(D))
    B_inv_D <- B * diag(inv_D)
		inv_O <- try(chol.inv(diag(q) +  t(B_inv_D) %*% B))
		if (any(class(inv_O) %in% "try-error"))
       return(loglike =
                 paste('ill-conditioned or singular Sigma'))

    inv_S <- try(inv_D - B_inv_D %*% inv_O %*% t(B_inv_D))
    if (any(class(inv_S) %in% "try-error"))
       return(loglike =
                   paste('ill-conditioned or singular Sigma'))
    for(i in 1 : g) {
		  logdetS <- sum(log(diag(D))) -  log(det(inv_O))
      mahal_dist <- mahalanobis(x = Y, center = mu[, i, drop=FALSE], 
                    cov = inv_S, inverted = TRUE)
      Fji[, i] <- -0.5 * mahal_dist - (p / 2) * log(2 * pi) - 0.5 * logdetS
    }
  }

} else {

  if (D_type == 'common') {
		
    inv_D <- diag(1/diag(D))
    for (i in 1 : g) {
			B_inv_D <- B[,, i] * diag(inv_D)
      inv_O <- try(chol.inv(diag(q) +  t(B_inv_D) %*% B[,, i]))
      if (any(class(inv_O) %in% "try-error"))
         return(loglike =
                paste('ill-conditioned or singular Sigma[,',i,']'))
      
      inv_S <- try(inv_D - B_inv_D %*% inv_O %*% t(B_inv_D))
      if (any(class(inv_S) %in% "try-error"))
         return(loglike =
                paste('ill-conditioned or singular Sigma[,',i,']'))
      
      logdetS <- sum(log(diag(D))) -  log(det(inv_O))
      mahal_dist <- mahalanobis(x = Y, center = mu[, i, drop=FALSE], 
                      cov = inv_S, inverted = TRUE)
      Fji[,i] <- -0.5 * mahal_dist - (p / 2) * log(2 * pi) - 0.5 * logdetS
    }
  }

  if (D_type == 'unique') {
  
		for (i in 1 : g) {
      inv_D <- diag(1/diag(D[, , i]))
      B_inv_D <- B[,, i] * diag(inv_D)
      inv_O <- try(chol.inv(diag(q) +  t(B_inv_D) %*% B[, , i]))
      if(any(class(inv_O) %in% "try-error"))
         return(loglike =
                 paste('ill-conditioned or singular Sigma[, ', i, ']'))

			inv_S <- try(inv_D - B_inv_D %*% inv_O %*% t(B_inv_D))
      if(any(class(inv_S) %in% "try-error"))
         return(loglike = paste('ill-conditioned or singular Sigma[,',i,']'))

			logdetS <- sum(log(diag(D[,, i]))) - log(det(inv_O))
      mahal_dist <- mahalanobis(Y, mu[, i, drop=FALSE], inv_S, TRUE)
      Fji[,i] <- -0.5 * mahal_dist - (p / 2) * log(2 * pi) - 0.5 * logdetS
    }
  }
}

Fji <- sweep(Fji, 2, log(pivec), '+')
Fjmax <- apply(Fji, 1, max)
Fji <- sweep(Fji, 1, Fjmax, '-')
Fji <- exp(Fji)
tau <- sweep(Fji, 1, rowSums(Fji), '/')
return(tau)
}