How to explain how I divided a bimodal distribution based on kernel density estimation

Question

I have a dataset of bimodal population. It contains a smaller peak, which is considered to be "bad", and a bigger peak. I try to separate the bad part of data from the rest of data. What I did was: first I did a kernel density estimation, then found the local maximum of this small peak, and the local minimum of the pit between two peaks, then I took the midpoint (arithmetic mean of x-coordinates) of them, and define it as a cutoff. Everything below this cutoff is considered "bad". The reason I took midpoint instead of pit is because I tried to be more conservative.

Now I would like to ask: Is what I did reasonable? If yes, how can I explain my action in a statisticians-favored manner? If not, how can I change? (Any other methods are welcome, especially those implemented in R.) Thank you!

Here's the figure.

Answer 1

You could fit a two-component mixture model using http://cran.r-project.org/web/packages/mixtools/index.html. Try using normalmixEM. You could then follow Erich Schubert's suggestions and find the region where Pr[data point was generated from the component with the smaller mean] >= 0.50.

Edit: example R code:

library(mixtools)

simulate <- function(lambda=0.3, mu=c(0, 4), sd=c(1, 1), n.obs=10^5) {
    x1 <- rnorm(n.obs, mu[1], sd[1])
    x2 <- rnorm(n.obs, mu[2], sd[2])    
    return(ifelse(runif(n.obs) < lambda, x1, x2))
}

x <- simulate()

model <- normalmixEM(x=x, k=2)
index.lower <- which.min(model$mu)  # Index of component with lower mean

find.cutoff <- function(proba=0.5, i=index.lower) {
    ## Cutoff such that Pr[drawn from bad component] == proba
    f <- function(x) {
        proba - (model$lambda[i]*dnorm(x, model$mu[i], model$sigma[i]) /
                     (model$lambda[1]*dnorm(x, model$mu[1], model$sigma[1]) + model$lambda[2]*dnorm(x, model$mu[2], model$sigma[2])))
        }
        return(uniroot(f=f, lower=-10, upper=10)$root)  # Careful with division by zero if changing lower and upper
}

cutoffs <- c(find.cutoff(proba=0.5), find.cutoff(proba=0.75))  # Around c(1.8, 1.5)

hist(x)
abline(v=cutoffs, col=c("red", "blue"), lty=2)

Answer 2

It would probably make more sense if you also estimated the "height" (actually, more of a weight) of both, and then set the threshold to the tipping point.

I.e. model the data as $$p_1 \cdot pdf(x, \mu_1, \sigma_1) + p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

and set the threshold to $x$ where $$p_1 \cdot pdf(x, \mu_1, \sigma_1) = p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

i.e. the object has the same chance of belonging to both classes.

You can still add a parameter to tune how conservative your method is, e.g. using $$p_1 \cdot pdf(x, \mu_1, \sigma_1) = c\cdot p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

where $c=2$ would put double weight on the second distribution.

Answer 3

I am using this example and I sometimes get this error

Error in uniroot(f = f, lower = -10, upper = 10) : f() values at end points not of opposite sign

and so i changed the lower to -1 and for some of the datasets it fixed it, but still errors out on others. Not sure if that can be dynamically set based on the input vector (i.e., x) ?