Mixture of Normal and Uniform to Construct 1D Gate

cutpoint.r

cut_mixture <- function(x, range = NULL, tol = 1e-6, maxit = 100) {
  n <- length(x)
  
  if (is.null(range)) {
    range <- c(min(x), max(x))
  }
  
  # Initialization
  # Initial value of pi1 are the proportion of x within 2 standard deviations of mu
  huber_out <- huber(x)
  mu <- huber_out$mu
  sigma2 <- huber_out$s^2
  pi1 <- mean(findInterval(x, mu + c(-2, 2) * sqrt(sigma2)) == 1)
  pi2 <- 1 - pi1
  
  Q_prev <- 0
  
  for (iter in seq_len(maxit)) {
    # E-Step
    r1 <- pi1 * dnorm(x, mean = mu, sd = sqrt(sigma2))
    r1 <- r1 / (r1 + pi2 * dunif(x, min = range[1], max = range[2]))
    r2 <- 1 - r1
    
    # M-Step
    sum_r1 <- sum(r1)
    pi1 <- sum_r1 / n
    pi2 <- 1 - pi1
    mu <- sum(r1 * x) / sum_r1
    sigma2 <- sum(r1 * (x - mu)^2) / sum_r1
    
    # Calculate the complete likelihood
    Q <- sum(matrix(log(c(pi1, pi2)), nrow = 1) %*% rbind(r1, r2))
    Q <- Q + r1 %*% dnorm(x, mean = mu, sd = sqrt(sigma2), log = TRUE)
    Q <- Q + r2 %*% dunif(x, min = range[1], max = range[2], log = TRUE)
    Q <- drop(Q)
    message(Q)
    
    if (abs(diff(c(Q_prev, Q))) < tol) {
      break
    } else {
      Q_prev <- Q
    }
  }
  
  probs <- cbind(r1, r2)
  x_range <- seq(mu, range[2], length = 1000)
  density_1 <- pi1 * dnorm(x_range, mean = mu, sd = sqrt(sigma2))
  density_2 <- pi2 * dunif(x_range, min = range[1], max = range[2])
  density_x <- density_1 + density_2
  cutpoint <- x_range[which(density_1 < density_2)[1]]
  
  list(probs = probs, density = density_x, cutpoint = cutpoint, mean = mu, var = sigma2)
}

sim.r

library(flowClust)
set.seed(42)
x <- SimulateMixture(N = 1000, w = c(0.6, 0.4), mu = matrix(c(2,2.5), ncol = 1),
                     sigma = array(c(0.05, 0.05, dim = c(2, 1, 1))))
x <- as.vector(x)
x <- x[x > 0]

mixture_out <- cut_mixture(x = x)

plot(density(x))
abline(v = mixture_out$cutpoint, col = "red")