Oddzielenie dwóch populacji od próbki

13

Próbuję oddzielić dwie grupy wartości od jednego zestawu danych. Mogę założyć, że jedna z populacji jest normalnie rozmieszczona i ma co najmniej połowę wielkości próbki. Wartości drugiego są zarówno niższe, jak i wyższe niż wartości pierwszego (rozkład jest nieznany). Staram się znaleźć górne i dolne granice, które obejmowałyby normalnie rozłożoną populację od drugiej.

Moje założenie zapewnia mi punkt wyjścia:

  • wszystkie punkty w zakresie międzykwartylowym próbki pochodzą z populacji normalnie rozmieszczonej.

Próbuję przetestować pod kątem wartości odstających, pobierając je z reszty próbki, dopóki nie zmieszczą się w 3 st.dev normalnie rozłożonej populacji. Co nie jest idealne, ale wydaje się, że daje wystarczająco rozsądny wynik.

Czy moje założenie jest uzasadnione statystycznie? Jaki byłby lepszy sposób na to?

ps proszę naprawić tagi kogoś.

SilentGhost
źródło
Czy możesz założyć, że pozostałe dwie grupy pochodzą z różnych rozkładów normalnych?
csgillespie
@cgillespie: to ta sama grupa, chyba tylko z dwoma trybami i dlatego prawdopodobnie nie mogę tego założyć.
SilentGhost
1
Czy wiesz, że członkowie drugiej grupy nie są uwzględnieni w pierwszej grupie, czy po prostu chcesz omyłkowo oznaczyć tych członków jako należących do pierwszej grupy?
Christian

Odpowiedzi:

10

Jeśli dobrze rozumiem, możesz po prostu dopasować do danych mieszaninę dwóch normalnych. Dostępnych jest wiele pakietów R. W tym przykładzie użyto pakietu 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)

To daje:

Mieszanka dwóch normalnych http://img294.imageshack.us/img294/4213/kernal.jpg

Pakiet zawiera również bardziej wyrafinowane metody - sprawdź dokumentację.

csgillespie
źródło
1
Załączony obraz wygasł.
naktinis,
3
  1. W przypadku danych w zakresie IQR należy użyć skróconego rozkładu normalnego (na przykład pakiet R gamlss.tr), aby oszacować parametry tego rozkładu.
  2. Innym podejściem jest stosowanie modeli mieszanin z 2 lub 3 komponentami (rozkładami). Możesz dopasować takie modele za pomocą pakietu gamlss.mx (dla każdego składnika mieszanki można określić dystrybucje z pakietu gamlss.dist).
Wojtek
źródło
2

Zakłada się, że nawet nie wiesz, czy druga dystrybucja jest normalna, czy nie; Zasadniczo radzę sobie z tą niepewnością, koncentrując się tylko na rozkładzie normalnym. To może być najlepsze podejście.

Jeśli możesz założyć, że dwie populacje są całkowicie oddzielone (tj. Wszystkie wartości z rozkładu A są mniejsze niż wszystkie wartości z rozkładu B), wówczas jednym z podejść jest użycie funkcji optimize () w R do wyszukania punktu przerwania, który daje oszacowania średniej i sd rozkładu normalnego, które sprawiają, że dane są najbardziej prawdopodobne:

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

Jeśli nie możesz założyć całkowitego rozdzielenia, myślę, że będziesz musiał założyć rozkład dla drugiego rozkładu, a następnie użyć modelowania mieszanki. Zauważ, że modelowanie mieszaniny nie będzie właściwie oznaczać poszczególnych punktów danych, ale da ci proporcję mieszaniny i oszacowanie parametrów każdego rozkładu (np. Średnia, sd, itp.).

Mike Lawrence
źródło
optimizejak rozumiem, wymaga dwóch dystrybucji obok siebie. W moim przypadku jedno jest w drugim, tzn. Wartości z drugiej populacji znajdują się po obu stronach limitów.
SilentGhost,
1

Dziwię się, że nikt nie zasugerował oczywistego rozwiązania:

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

Teraz wyjaśnienie: ltsRegfunkcja w pakiecie robustbase, gdy zostanie wywołana z opcją

nsamp="best"

daje jednowymiarowe (dokładne) wagi MCD. (są to wagi n-wektor 0-1 przechowywane w $raw.weightsobiekcie. Algorytmem do ich identyfikacji jest estymator MCD (1)).

W skrócie, wagi te wynoszą 1 dla członków podzbioru najbardziej skoncentrowanych obserwacji.h=(n+2)/2

W jednym wymiarze, zaczyna się przez sortowanie wszystkie obserwacje następnie oblicza miarę stycznymi podzbiorów obserwacji: oznaczający wprowadzenie wektora sortowanych obserwacji, to oblicza miarę ( np. następnie i tak dalej ... ) zachowuje ten o mniejszej miary.x ( i ) i T h ( x ( 1 ) , . . . , x ( h + 1 ) ) ( x ( 2 ) , . . . , x ( H + 2 ) )hx(i)ith
(x(1),...,x(h+1))(x(2),...,x(h+2))

Algorytm ten zakłada, że ​​twoja grupa zainteresowań stanowi zdecydowaną większość oryginalnej próbki i że ma symetryczny rozkład (ale nie ma hipotezy o rozkładzie pozostałych obserwacji ).nh

(1) PJ Rousseeuw (1984). Najmniejsza mediana regresji kwadratów, Journal of American Statistics Association.

użytkownik603
źródło