Trennen von zwei Populationen von der Probe

Ich versuche, zwei Wertegruppen von einem einzigen Datensatz zu trennen. Ich kann davon ausgehen, dass eine der Populationen normal verteilt ist und mindestens halb so groß wie die Stichprobe ist. Die Werte der zweiten sind beide niedriger oder höher als die Werte der ersten (Verteilung ist unbekannt). Was ich versuche, ist, die oberen und unteren Grenzen zu finden, die die normalverteilte Bevölkerung von der anderen einschließen würden.

Meine Vermutung gibt mir den Ausgangspunkt:

• Alle Punkte im Interquartilbereich der Stichprobe stammen aus der normalverteilten Grundgesamtheit.

Ich versuche, auf Ausreißer zu testen, indem ich sie aus dem Rest der Stichprobe nehme, bis sie nicht mehr in die 3. Stufe der normalverteilten Population passen. Was nicht ideal ist, aber vernünftig genug zu sein scheint.

Ist meine Annahme statistisch gesichert? Was wäre ein besserer Weg, dies zu tun?

ps bitte die tags jemand reparieren.

Wenn ich das richtig verstehe, können Sie den Daten einfach eine Mischung aus zwei Normalen zuordnen. Es gibt viele R-Pakete, die dazu zur Verfügung stehen. In diesem Beispiel wird das mixtools- Paket verwendet:

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

Das gibt:

Mischung aus zwei Normalen http://img294.imageshack.us/img294/4213/kernal.jpg

Das Paket enthält auch komplexere Methoden - überprüfen Sie die Dokumentation.

1. Für Daten im IQR-Bereich sollten Sie die abgeschnittene Normalverteilung (z. B. R-Paket gamlss.tr) verwenden, um die Parameter dieser Verteilung zu schätzen.
2. Ein anderer Ansatz ist die Verwendung von Mischungsmodellen mit 2 oder 3 Komponenten (Verteilungen). Sie können solche Modelle mit dem Paket gamlss.mx anpassen (Distributionen aus dem Paket gamlss.dist können für jede Komponente der Mischung angegeben werden).

Dies setzt voraus, dass Sie nicht einmal wissen, ob die zweite Verteilung normal ist oder nicht. Grundsätzlich gehe ich mit dieser Unsicherheit um, indem ich mich nur auf die Normalverteilung konzentriere. Dies kann der beste Ansatz sein oder auch nicht.

Wenn Sie davon ausgehen können, dass die beiden Populationen vollständig getrennt sind (dh, alle Werte aus der Verteilung A sind kleiner als alle Werte aus der Verteilung B), besteht ein Ansatz darin, die Optimize () - Funktion in R zu verwenden, um nach dem entsprechenden Haltepunkt zu suchen liefert Schätzungen des Mittelwerts und des SD der Normalverteilung, die die Daten am wahrscheinlichsten machen:

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

Wenn Sie nicht von einer vollständigen Trennung ausgehen können, müssen Sie meiner Meinung nach für die zweite Verteilung eine gewisse Verteilung annehmen und dann die Mischungsmodellierung verwenden. Beachten Sie, dass die Gemischmodellierung die einzelnen Datenpunkte nicht wirklich kennzeichnet, sondern Ihnen den Mischungsanteil und die Schätzungen der Parameter jeder Verteilung (z. B. Mittelwert, SD usw.) liefert.

Ich bin überrascht, dass niemand die offensichtliche Lösung vorgeschlagen hat:

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

Nun zur Erklärung: Die ltsRegFunktion im Paket robustbase, wenn sie mit der Option aufgerufen wird

nsamp="best"

ergibt die univariaten (exakten) MCD-Gewichte. (Dies sind 0-1 Gewichte eines n-Vektors, die in der$raw.weights Objekt . Der Algorithmus zur Identifizierung ist der MCD-Schätzer (1)).

Kurz gesagt, diese Gewichte sind 1 für die Mitglieder der Teilmenge von $h=\lceil(n+2)/2\rceil$ konzentriertesten Beobachtungen.

In der ersten Dimension werden zunächst alle Beobachtungen sortiert und dann das Maß aller zusammenhängenden Teilmengen von berechnet $h$ Beobachtungen: bezeichnen $x_{(i)}$ das $i^{th}$Eingabe des Vektors sortierter Beobachtungen, errechnet es das Maß von
(z$(x_{(1)},...,x_{(h+1)})$ dann $(x_{(2)},...,x_{(h+2)})$ und so weiter ...) behält dann den mit kleinerem Maß.

Bei diesem Algorithmus wird davon ausgegangen, dass Ihre Interessengruppe den größten Teil der ursprünglichen Stichprobe ausmacht und eine symmetrische Verteilung aufweist (es gibt jedoch keine Hypothese zur Verteilung der verbleibenden $n-h$ Überwachung).

(1) PJ Rousseeuw (1984). Geringster Median der Quadratregression, Journal of the American Statistical Association.