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

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.

• how did you do to found the local maximum of this small peak, and the local minimum of the pit between two peaks?
– Lea
Oct 5, 2013 at 16:40

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)  • This is excellent. However, could you explain why there are two cutoffs, rather than just a single value for which$P(Distribution_1) = P(Distribution_2) = .5$? – Eoin May 14, 2014 at 20:41 • Too late to Edit: I figured it out (I think) - the first value is the 50/50 probability, the second is 75/25, set by the values passed to find.cutoff. – Eoin May 14, 2014 at 21:14 • That's right. The function find.cutoff(prob) gives you a point at which there is a a given probability that the observation was drawn from the lower component, i.e. it returns the x such that Pr[X came from lower component | X=x] = prob. May 15, 2014 at 7:59 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.

• Thank you Erich! This is a good piece of advice! Can you elaborate a bit how I can calculate tipping point (inflection point?)? I use R. The turning points are computed using the "pastecs" package. However I Googled but couldn't find a package to compute tipping points of a density function.
– 琼 朱
May 3, 2013 at 19:41
• It's not a technical term. Does the turning point work for you? May 6, 2013 at 7:10

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) ?

• Welcome @StudentofScience. Not sure whether this constitutes an answer to the OP's question. Would you mind giving more explanation about your idea and how this adresses the question?
– Momo
Dec 3, 2013 at 9:52
• Hi, I did not have enough rep to add a comment to ask Adrian but it helps in making the answer better i believe because the answer hit an error. Dec 4, 2013 at 6:27