Cutpoint Method using Mixture of Two Normals

cutpoint.r

normalmix_loglike <- function(x, y) {
  y <- as.factor(y)
  x1 <- x[y == levels(y)[1]]
  x2 <- x[y == levels(y)[2]]
  
  n1 <- length(x1)
  n2 <- length(x2)
  n <- n1 + n2
  
  w1 <- n1 / n
  w2 <- n2 / n
  
  xbar1 <- mean(x1)
  xbar2 <- mean(x2)
  
  var1 <- (n1 - 1) * var(x1) / n1
  var2 <- (n2 - 1) * var(x2) / n2
  
  f1 <- dnorm(x, mean = xbar1, sd = sqrt(var1))
  f2 <- dnorm(x, mean = xbar2, sd = sqrt(var2))
  
  sum(log(w1 * f1 + w2 * f2))
}


gate1d_loglike <- function(x, num_candidates = 100) {
  x <- as.vector(x)
  cutpoints <- seq(from = min(x), to = max(x), length = num_candidates)
  
  sapply(cutpoints, function(cutpoint) {
    y <- findInterval(x, cutpoint)
    normalmix_loglike(x, y)
  })
}

library(flowClust)
set.seed(42)

# Data with a shoulder
x <- SimulateMixture(N = 10000, w = c(0.6, 0.4), mu = matrix(c(0, 2.5), nrow = 1), sigma = array(c(1, 1), dim = c(2, 1, 1)), nu = 30)
gate_out <- gate1d_loglike(x, 500)

plot(density(x, adjust = 2))
abline(v = seq(from = min(x), to = max(x), length = 500)[which.max(gate_out)], col = "red")

plot(gate_out)