Separar dos poblaciones de la muestra.

Estoy tratando de separar dos grupos de valores de un solo conjunto de datos. Puedo suponer que una de las poblaciones está normalmente distribuida y tiene al menos la mitad del tamaño de la muestra. Los valores del segundo son más bajos o más altos que los valores del primero (se desconoce la distribución). Lo que intento hacer es encontrar los límites superior e inferior que encerrarían a la población normalmente distribuida de la otra.

Mi suposición me proporciona un punto de partida:

• Todos los puntos dentro del rango intercuartil de la muestra son de la población normalmente distribuida.

Estoy tratando de detectar valores atípicos tomándolos del resto de la muestra hasta que no encajen en el 3 st.dev de la población normalmente distribuida. Lo cual no es ideal, pero parece producir un resultado lo suficientemente razonable.

¿Mi suposición es estadísticamente sólida? ¿Cuál sería una mejor manera de hacer esto?

PD: arregla las etiquetas a alguien.

Si entiendo correctamente, entonces solo puede ajustar una mezcla de dos normales a los datos. Hay muchos paquetes R disponibles para hacer esto. Este ejemplo usa el paquete mixtools :

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

Esto da:

Mezcla de dos normales http://img294.imageshack.us/img294/4213/kernal.jpg

El paquete también contiene métodos más sofisticados: consulte la documentación.

1. Para los datos en el rango IQR, debe usar la distribución normal truncada (por ejemplo, el paquete R gamlss.tr) para estimar los parámetros de esta distribución.
2. Otro enfoque es utilizar modelos de mezcla con 2 o 3 componentes (distribuciones). Puede ajustar dichos modelos usando el paquete gamlss.mx (las distribuciones del paquete gamlss.dist se pueden especificar para cada componente de la mezcla).

Esto supone que ni siquiera sabe si la segunda distribución es normal o no; Básicamente manejo esta incertidumbre enfocándome solo en la distribución normal. Este puede o no ser el mejor enfoque.

Si puede suponer que las dos poblaciones están completamente separadas (es decir, todos los valores de la distribución A son menores que todos los valores de la distribución B), entonces un enfoque es utilizar la función optimizar () en R para buscar el punto de ruptura que produce estimaciones de la media y la desviación estándar de la distribución normal que hacen que los datos sean más probables:

#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)[2],max(mix))
    , maximum = TRUE
)$maximum

#label the data
labelled_mix = data.frame(
    x = mix
    , source = ifelse(mix<break_point,'A','B')
)
print(labelled_mix)

Si no puede asumir una separación completa, entonces creo que tendrá que asumir alguna distribución para la segunda distribución y luego usar el modelado de mezclas. Tenga en cuenta que el modelado de mezclas en realidad no etiquetará los puntos de datos individuales, pero le dará la proporción de la mezcla y las estimaciones de los parámetros de cada distribución (por ejemplo, media, SD, etc.).

Me sorprende que nadie sugiera la solución obvia:

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

Ahora para la explicación: la ltsRegfunción en el paquete robustbase, cuando se llama con la opción

nsamp="best"

produce los pesos univariados (exactos) de MCD. (estos son pesos n-vector 0-1 almacenados en el $raw.weightsobjeto. El algoritmo para identificarlos es el estimador MCD (1)).

En pocas palabras, estos pesos son 1 para los miembros del subconjunto de $h=\lceil(n+2)/2\rceil$ observaciones más concentradas.

En la dimensión uno, comienza ordenando todas las observaciones y luego calcula la medida de todos los subconjuntos contiguos de $h$ observaciones: denotando $x_{(i)}$ el $i^{th}$entrada del vector de observaciones ordenadas, calcula la medida de
(p. ej.$(x_{(1)},...,x_{(h+1)})$ luego $(x_{(2)},...,x_{(h+2)})$ y así sucesivamente ...) luego retiene el que tiene una medida menor.

Este algoritmo supone que su grupo de interés cuenta con una mayoría estricta de la muestra original y que tiene una distribución simétrica (pero no hay hipótesis sobre la distribución del resto) $n-h$ observación).

(1) PJ Rousseeuw (1984). Mínima mediana de regresión de cuadrados, Journal of the American Statistical Association.