R/normalmixEM2comp2.R

In PheNorm: Unsupervised Gold-Standard Label Free Phenotyping Algorithm for EHR Data

normalmixEM2comp2 <- function (x, lambda, mu, sigsqrd, eps = 1e-08, maxit = 1000, verb = FALSE) {
  arbvar <- (length(sigsqrd) == 2)
  mu1 <- mu[1]
  mu2 <- mu[2]
  sigsqrd1 <- sigsqrd[1]
  sigsqrd2 <- sigsqrd[arbvar + 1]
  mx <- mean(x)
  const <- length(x) * 0.918938533204673
  dl <- 1 + eps
  iter <- 0
  ll <- rep(0, maxit + 1)
  a1 <- (x - mu1)^2
  b1 <- (lambda/sqrt(sigsqrd1)) * exp(-a1/2/sigsqrd1)
  a2 <- (x - mu2)^2
  b2 <- ((1 - lambda)/sqrt(sigsqrd2)) * exp(-a2/2/sigsqrd2)
  l <- sum(log(b1 + b2))
  while (dl > eps && iter < maxit) {
    iter <- iter + 1
    ll[iter] <- l
    postprobs <- b1/(b1 + b2)
    lambda <- mean(postprobs)
    mu1 <- mean(postprobs * x)/lambda
    mu2 <- (mx - lambda * mu1)/(1 - lambda)
    if (arbvar) {
      sigsqrd1 <- mean(postprobs * a1)/lambda
      sigsqrd2 <- mean((1 - postprobs) * a2)/(1 - lambda)
    }
    else {
      sigsqrd1 <- sigsqrd2 <- mean(postprobs * a1 + (1 - postprobs) * a2)
    }
    a1 <- (x - mu1)^2
    b1 <- (lambda/sqrt(sigsqrd1)) * exp(-a1/2/sigsqrd1)
    a2 <- (x - mu2)^2
    b2 <- ((1 - lambda)/sqrt(sigsqrd2)) * exp(-a2/2/sigsqrd2)
    oldl <- l
    l <- sum(log(b1 + b2))
    dl <- l - oldl
    if (verb) {
      cat("iteration =", iter, " log-lik diff =", dl, " log-lik =",
          l - const, "\n")
    }
  }
  ##cat("number of iterations=", iter, "\n")
  iter <- iter + 1
  ll[iter] <- l
  postprobs <- cbind(postprobs, 1 - postprobs)
  colnames(postprobs) <- c(paste("comp", ".", 1:2, sep = ""))
  out <- list(x = x, lambda = c(lambda, 1 - lambda), mu = c(mu1,mu2),
              sigma = sqrt(c(sigsqrd1, sigsqrd2)[1:(1 + arbvar)]),
              loglik = l - const, posterior = postprobs, all.loglik = ll[1:iter] - const,
              restarts = 0, ft = "normalmixEM")
  class(out) <- "mixEM"
  out
}