Séparer deux populations de l'échantillon

J'essaie de séparer deux groupes de valeurs d'un même ensemble de données. Je peux supposer que l'une des populations est normalement distribuée et représente au moins la moitié de la taille de l'échantillon. Les valeurs du second sont à la fois inférieures ou supérieures aux valeurs du premier (la distribution est inconnue). Ce que j'essaie de faire, c'est de trouver les limites supérieure et inférieure qui enfermeraient la population normalement répartie de l'autre.

Mon hypothèse me fournit un point de départ:

• tous les points dans l'intervalle interquartile de l'échantillon proviennent de la population normalement distribuée.

J'essaie de tester les valeurs aberrantes en les tirant du reste de l'échantillon jusqu'à ce qu'elles ne rentrent pas dans le 3 e rang de la population normalement distribuée. Ce qui n'est pas idéal, mais semble produire un résultat assez raisonnable.

Mon hypothèse est-elle statistiquement valable? Quelle serait la meilleure façon de procéder?

ps veuillez corriger les balises quelqu'un.

Si je comprends bien, vous pouvez simplement adapter un mélange de deux normales aux données. Il existe de nombreux packages R disponibles pour ce faire. Cet exemple utilise le package 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)

Cela donne:

Mélange de deux normales http://img294.imageshack.us/img294/4213/kernal.jpg

Le paquet contient également des méthodes plus sophistiquées - consultez la documentation.

1. Pour les données dans la plage IQR, vous devez utiliser la distribution normale tronquée (par exemple le package R gamlss.tr) pour estimer les paramètres de cette distribution.
2. Une autre approche consiste à utiliser des modèles de mélange à 2 ou 3 composants (distributions). Vous pouvez adapter ces modèles à l'aide du package gamlss.mx (les distributions du package gamlss.dist peuvent être spécifiées pour chaque composant du mélange).

Cela suppose que vous ne savez même pas si la deuxième distribution est normale ou non; Je gère essentiellement cette incertitude en me concentrant uniquement sur la distribution normale. Cela peut ou non être la meilleure approche.

Si vous pouvez supposer que les deux populations sont complètement séparées (c'est-à-dire que toutes les valeurs de la distribution A sont inférieures à toutes les valeurs de la distribution B), alors une approche consiste à utiliser la fonction optimise () dans R pour rechercher le point de rupture qui donne des estimations de la moyenne et du sd de la distribution normale qui rendent les données les plus 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 vous ne pouvez pas supposer une séparation complète, je pense que vous devrez assumer une certaine distribution pour la deuxième distribution, puis utiliser la modélisation du mélange. Notez que la modélisation du mélange ne marquera pas réellement les points de données individuels, mais vous donnera la proportion du mélange et les estimations des paramètres de chaque distribution (par exemple, moyenne, sd, etc.).

Je suis surpris que personne n'ait suggéré la solution évidente:

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

Maintenant pour l'explication: la ltsRegfonction dans le package robustbase, lorsqu'elle est appelée avec l'option

nsamp="best"

donne les poids MCD univariés (exacts). (il s'agit d'un poids n-vecteur 0-1 stocké dans l' $raw.weightsobjet. L'algorithme pour les identifier est l'estimateur MCD (1)).

En un mot, ces poids sont de 1 pour les membres du sous-ensemble de $h=\lceil(n+2)/2\rceil$ observations les plus concentrées.

Dans la dimension un, il commence par trier toutes les observations puis calcule la mesure de tous les sous-ensembles contigus de $h$ observations: dénoter $x_{(i)}$ le $i^{th}$entrée du vecteur d'observations triées, il calcule la mesure de
(par exemple$(x_{(1)},...,x_{(h+1)})$ ensuite $(x_{(2)},...,x_{(h+2)})$ et ainsi de suite ...) conserve alors celui avec une mesure plus petite.

Cet algorithme suppose que votre groupe d'intérêt numérote une majorité stricte de l'échantillon d'origine et qu'il a une distribution symétrique (mais il n'y a pas d'hypothèse sur la distribution du reste $n-h$ observation).

(1) PJ Rousseeuw (1984). La moindre médiane de régression des carrés, Journal of the American Statistical Association.