- # For each lambda, computation of the parameters
- out =
- if (ncores > 1)
- parLapply(cl, 1:length(S), computeAtLambda)
- else
- lapply(1:length(S), computeAtLambda)
+ ## Computation of the loglikelihood
+ # Precompute det(rhoLambda[,,r]) for r in 1...k
+ detRho <- sapply(1:k, function(r) gdet(rhoLambda[, , r]))
+ sumLogLLH <- 0
+ for (i in 1:n)
+ {
+ # Update gam[,]; use log to avoid numerical problems
+ logGam <- sapply(1:k, function(r) {
+ log(piLambda[r]) + log(detRho[r]) - 0.5 *
+ sum((Y[i, ] %*% rhoLambda[, , r] - X[i, ] %*% phiLambda[, , r])^2)
+ })
+
+ logGam <- logGam - max(logGam) #adjust without changing proportions
+ gam <- exp(logGam)
+ norm_fact <- sum(gam)
+ sumLogLLH <- sumLogLLH + log(norm_fact) - log((2 * base::pi)^(m/2))
+ }
+ llhLambda <- c(sumLogLLH/n, (dimension + m + 1) * k - 1)
+ list(phi = phiLambda, rho = rhoLambda, pi = piLambda, llh = llhLambda)
+ }