7
$\begingroup$

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.

enter image description here

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

3 Answers 3

8
$\begingroup$

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)
$\endgroup$
3
  • $\begingroup$ 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$? $\endgroup$
    – Eoin
    May 14, 2014 at 20:41
  • 1
    $\begingroup$ 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. $\endgroup$
    – Eoin
    May 14, 2014 at 21:14
  • 1
    $\begingroup$ 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. $\endgroup$
    – Adrian
    May 15, 2014 at 7:59
2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – 琼 朱
    May 3, 2013 at 19:41
  • $\begingroup$ It's not a technical term. Does the turning point work for you? $\endgroup$ May 6, 2013 at 7:10
0
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$
    – Momo
    Dec 3, 2013 at 9:52
  • $\begingroup$ 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. $\endgroup$ Dec 4, 2013 at 6:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.