- # Update gam[,]
- for (r in 1:k)
- {
- gam[i, r] <- pi[r] * exp(-0.5
- * sum((Y[i, ] %*% rho[, , r] - X[i, ] %*% phi[, , r])^2)) * detRho[r]
- }
+ # Update gam[,]; use log to avoid numerical problems
+ logGam <- sapply(1:k, function(r) {
+ log(pi[r]) + log(detRho[r]) - 0.5 *
+ sum((Y[i, ] %*% rho[, , r] - X[i, ] %*% phi[, , r])^2)
+ })
+
+ logGam <- logGam - max(logGam) #adjust without changing proportions
+ gam[i, ] <- exp(logGam)
+ norm_fact <- sum(gam[i, ])
+ gam[i, ] <- gam[i, ] / norm_fact
+ sumLogLLH <- sumLogLLH + log(norm_fact) - log((2 * base::pi)^(m/2))