Separating two populations from the sample

I'm trying to separate two groups of values from a single data set. I can assume that one of the populations is normally distributed and is at least half the size of the sample. The values of the second one are both lower or higher than the values from the first one (distribution is unknown). What I'm trying to do is to find the upper and lower limits that would enclose the normally-distributed population from the other.

My assumption provide me with starting point:

• all points within the interquartile range of the sample are from the normally-distributed population.

I'm trying to test for outliers taking them from the rest of the sample until they don't fit into the 3 st.dev of the normally-distributed population. Which is not ideal, but seem to produce reasonable enough result.

p.s. please fix the tags someone.

• Can you assume that the other two groups are from different Normal distributions? Jul 28 '10 at 14:37
• @cgillespie: it is the same group, just with two modes, I guess, and therefore I probably cannot assume this. Jul 28 '10 at 16:03
• Do you know that members of the second group aren't included in the first group or are you just willing to mistakenly label those members as belonging to the first group? Jul 29 '10 at 20:54

If I understand correctly, then you can just fit a mixture of two Normals to the data. There are lots of R packages that are available to do this. This example uses the mixtools package:

#Taken from the documentation
library(mixtools)
data(faithful)
attach(faithful)

#Fit two Normals
wait1 = normalmixEM(waiting, lambda = 0.5)
plot(wait1, density=TRUE, loglik=FALSE)

This gives:

Mixture of two Normals http://img294.imageshack.us/img294/4213/kernal.jpg

The package also contains more sophisticated methods - check the documentation.

• The image you've attached has expired. Oct 10 '15 at 9:13
1. For data in IQR range you should use truncated normal distribution (for example R package gamlss.tr) to estimate parameters of this distribution.
2. Another approach is using mixture models with 2 or 3 components (distributions). You can fit such models using gamlss.mx package (distributions from package gamlss.dist can be specified for each component of mixture).

This assumes that you don't even know if the second distribution is normal or not; I basically handle this uncertainty by focusing only on the normal distribution. This may or may not be the best approach.

If you can assume that the two populations are completely separated (i.e., all values from distribution A are less than all values from distribution B), then one approach is to use the optimize() function in R to search for the break-point that yields estimates of the mean and sd of the normal distribution that make the data most likely:

#generate completely separated data
a = rnorm(100)
b = rnorm(100,10)
while(!all(a<b)){
a = rnorm(100)
b = rnorm(100,10)
}

#create a mix
mix = c(a,b)

#"forget" the original distributions
rm(a)
rm(b)

#try to find the break point between the distributions
break_point = optimize(
f = function(x){
data_from_a = mix[mix<x]
likelihood = dnorm(data_from_a,mean(data_from_a),sd(data_from_a))
SLL = sum(log(likelihood))
return(SLL)
}
, interval = c(sort(mix),max(mix))
, maximum = TRUE
)$maximum #label the data labelled_mix = data.frame( x = mix , source = ifelse(mix<break_point,'A','B') ) print(labelled_mix) If you can't assume complete separation, then I think you'll have to assume some distribution for the second distribution and then use mixture modelling. Note that mixture modelling won't actually label the individual data points, but will give you the mixture proportion and estimates of the parameters of each distribution (eg. mean, sd, etc.). • optimize requires two distributions to be side-by-side as I understand. In my case one is inside the other, i.e., the values from the second population are on both side of the limits. Jul 28 '10 at 16:09 I'm surprised nobody suggested the obvious solution: #generate completely separated data library(robustbase) set.seed(123) x<-rnorm(200) x[1:40]<-x[1:40]+10 x[41:80]<-x[41:80]-10 Rob<-ltsReg(x~1,nsamp="best") #all the good guys which(Rob$raw.weights==1)

Now for the explanation: the ltsReg function in package robustbase, when called with the option

nsamp="best"

yields the univariate (exact) MCD weights. (these are a n-vector 0-1 weights stored in the $raw.weights object. The algorithm to identify them is the MCD estimator (1)). In a nutshell, these weights are 1 for the members of the subset of$h=\lceil(n+2)/2\rceil$most concentrated observations. In dimension one, it starts by sorting all the observations then computes the measure of all contiguous subsets of$h$observations: denoting$x_{(i)}$the$i^{th}$entry of the vector of sorted observations, it computes the measure of (e.g.$(x_{(1)},...,x_{(h+1)})$then$(x_{(2)},...,x_{(h+2)})$and so forth...) then retains the one with smaller measure. This algorithm assumes that your group of interest numbers a strict majority of the original sample and that it has a symmetrical distribution (but there no no hypothesis on the distribution of the remaining$n-h\$ observation).

(1) P.J. Rousseeuw (1984). Least median of squares regression, Journal of the American Statistical Association.