Jak sprawdzić, czy moja dystrybucja jest multimodalna?

21

Kiedy wykreślam histogram moich danych, ma on dwa szczyty:

histogram

Czy to oznacza potencjalny rozkład multimodalny? Uruchomiłem dip.testw R ( library(diptest)), a dane wyjściowe to:

D = 0.0275, p-value = 0.7913

Czy mogę stwierdzić, że moje dane mają rozkład multimodalny?

DANE

10346 13698 13894 19854 28066 26620 27066 16658  9221 13578 11483 10390 11126 13487 
15851 16116 24102 30892 25081 14067 10433 15591  8639 10345 10639 15796 14507 21289 
25444 26149 23612 19671 12447 13535 10667 11255  8442 11546 15958 21058 28088 23827 
30707 19653 12791 13463 11465 12326 12277 12769 18341 19140 24590 28277 22694 15489 
11070 11002 11579  9834  9364 15128 15147 18499 25134 32116 24475 21952 10272 15404 
13079 10633 10761 13714 16073 23335 29822 26800 31489 19780 12238 15318  9646 11786 
10906 13056 17599 22524 25057 28809 27880 19912 12319 18240 11934 10290 11304 16092 
15911 24671 31081 27716 25388 22665 10603 14409 10736  9651 12533 17546 16863 23598 
25867 31774 24216 20448 12548 15129 11687 11581
użytkownik1260391
źródło
3
Użyj więcej pojemników na swoim histogramie. Sugeruję około dwa razy więcej
Glen_b
1
W tej odpowiedzi wspomniano o dziewięciu różnych testach , z których niektóre mogą być odpowiednie dla twojej sytuacji.
Glen_b
1
Ten artykuł może ci się przydać, jeśli jeszcze go nie widziałeś (także ten ciąg dalszy )
Eoin

Odpowiedzi:

15

@NickCox przedstawił interesującą strategię (+1). Mogę jednak uznać to za bardziej odkrywcze z uwagi na obawy, które wskazuje @whuber .

Pozwól, że zasugeruję inną strategię: możesz dopasować model Gaussa do mieszanki skończonej. Zauważ, że powoduje to bardzo silne założenie, że twoje dane pochodzą z co najmniej jednej prawdziwej normy. Jak zarówno @whuber, jak i @NickCox wskazują w komentarzach, bez merytorycznej interpretacji tych danych - popartej dobrze ugruntowaną teorią - na poparcie tego założenia, strategię tę należy również uznać za eksploracyjną.

Najpierw zastosujmy się do sugestii @ Glen_b i spójrz na twoje dane, używając dwa razy więcej pojemników:

wprowadź opis zdjęcia tutaj

Nadal widzimy dwa tryby; jeśli już, to przejawiają się tu wyraźniej. (Należy również zauważyć, że linia gęstości jądra powinna być identyczna, ale wydaje się bardziej rozłożona z powodu większej liczby pojemników).

Teraz dopasujmy model mieszanki skończonej Gaussa. W Rtym celu możesz użyć Mclustpakietu:

library(mclust)
x.gmm = Mclust(x)
summary(x.gmm)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
#   
# Mclust V (univariate, unequal variance) model with 2 components:
#   
#   log.likelihood   n df       BIC       ICL
#        -1200.874 120  5 -2425.686 -2442.719
# 
# Clustering table:
#  1  2 
# 68 52 

Dwa normalne komponenty optymalizują BIC. Dla porównania możemy wymusić dopasowanie jednego elementu i przeprowadzić test współczynnika wiarygodności:

x.gmm.1 = Mclust(x, G=1)
logLik(x.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(x.gmm)-logLik(x.gmm.1)
# 'log Lik.' 25.36657 (df=5)
1-pchisq(25.36657, df=3)  # [1] 1.294187e-05

Sugeruje to, że jest bardzo mało prawdopodobne, abyś znalazł dane tak dalekie od unimodalnego jak twoje, gdyby pochodziły z jednego prawdziwego rozkładu normalnego.

Niektórzy ludzie nie czują się komfortowo, używając testu parametrycznego (chociaż jeśli założenia się utrzymują, nie znam żadnego problemu). Jedną z bardzo szeroko stosowanych technik jest użycie metody parametrycznego dopasowania krzyżowego (opisuję tutaj algorytm ). Możemy spróbować zastosować to do tych danych:

x.gmm$parameters
# $mean
# 12346.98 23322.06 
# $variance$sigmasq
# [1]  4514863 24582180
x.gmm.1$parameters
# $mean
# [1] 17520.91
# $variance$sigmasq
# [1] 43989870

set.seed(7809)
B = 10000;    x2.d = vector(length=B);    x1.d = vector(length=B)
for(i in 1:B){
  x2      = c(rnorm(68, mean=12346.98, sd=sqrt( 4514863)), 
              rnorm(52, mean=23322.06, sd=sqrt(24582180)) )
  x1      = rnorm( 120, mean=17520.91, sd=sqrt(43989870))
  x2.d[i] = Mclust(x2, G=2)$loglik - Mclust(x2, G=1)$loglik
  x1.d[i] = Mclust(x1, G=2)$loglik - Mclust(x1, G=1)$loglik
}
x2.d = sort(x2.d);  x1.d = sort(x1.d)
summary(x1.d)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.29070 -0.02124  0.41460  0.88760  1.36700 14.01000 
summary(x2.d)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  9.006  23.770  27.500  27.760  31.350  53.500 

wprowadź opis zdjęcia tutaj

Statystyki podsumowujące i wykresy gęstości jądra dla rozkładów próbkowania pokazują kilka interesujących cech. Prawdopodobieństwo dziennika dla modelu jednoskładnikowego rzadko jest większe niż w przypadku dopasowania dwuskładnikowego, nawet jeśli proces generowania prawdziwych danych ma tylko jeden składnik, a gdy jest większy, ilość jest trywialna. Pomysł porównania modeli różniących się zdolnością do dopasowania danych jest jedną z motywów stojących za PBCM. Te dwa rozkłady próbkowania prawie się nie pokrywają; tylko 0,35% z nich x2.djest mniej niż maksimumx1.dwartość. Jeśli wybrałeś model dwuskładnikowy, jeśli różnica w prawdopodobieństwie logarytmu wynosiła> 9,7, niepoprawnie wybrałbyś model jednoskładnikowy .01% i model dwuskładnikowy w 0,02% przypadków. Są to wysoce dyskryminujące. Jeśli natomiast zdecydujesz się zastosować model jednoskładnikowy jako hipotezę zerową, zaobserwowany wynik jest wystarczająco mały, aby nie pojawił się w empirycznym rozkładzie próbkowania w 10 000 iteracjach. Możemy użyć reguły 3 (patrz tutaj ), aby ustalić górną granicę wartości p, a mianowicie szacujemy, że twoja wartość p jest mniejsza niż 0,0003. To jest bardzo istotne.

Rodzi to pytanie, dlaczego wyniki tak bardzo odbiegają od testu zanurzeniowego. (Aby odpowiedzieć na twoje wyraźne pytanie, twój test zanurzeniowy nie dostarcza dowodów na istnienie dwóch rzeczywistych trybów.) Szczerze mówiąc, nie znam testu zanurzeniowego, więc trudno powiedzieć; to może być słabe. Myślę jednak, że prawdopodobną odpowiedzią jest to, że takie podejście zakłada, że ​​twoje dane są generowane przez prawdziwe wartości normalne. Test Shapiro-Wilka dla danych jest bardzo istotny ( ), a także bardzo istotny dla optymalnej transformacji danych Box-Coxa (odwrotny pierwiastek kwadratowy; ). Jednak dane nigdy nie są tak naprawdę normalne (por. Ten słynny cytatp<.000001p<.001), a komponenty leżące u ich podstaw, jeśli takie istnieją, nie są również gwarantowane jako całkowicie normalne. Jeśli uznasz, że uzasadnione jest, że twoje dane mogą pochodzić z dodatnio wypaczonego rozkładu, a nie normalnego, ten poziom bimodalności może również mieścić się w typowym zakresie zmienności, co, jak podejrzewam, mówi test zanurzeniowy.

gung - Przywróć Monikę
źródło
1
Problem z tym podejściem polega na tym, że alternatywa, z którą porównujesz mieszaninę Gaussa, nie jest zbyt rozsądna. Bardziej rozsądne byłoby to, że dystrybucja jest rodzajem przesuniętym w prawo, takim jak gamma. Jest prawie pewne, że mieszanina zmieści się w prawie każdym wypaczonym zbiorze danych „znacznie” lepiej niż pojedynczy Gaussian .
whuber
Masz rację, @whuber. Starałem się to wyraźnie podkreślić. Nie jestem pewien, jak zrobić Gamma FMM, ale byłoby lepiej.
Gung - Przywróć Monikę
1
Ponieważ jest to eksploracyjne, jedną z myśli byłoby próba przekształcenia pierwotnego rozkładu w symetrię (być może z przesuniętą transformacją Box-Coxa, solidnie oszacowaną na podstawie kilku kwantyli danych) i ponowna próba podejścia. Oczywiście nie mówilibyście o „znaczeniu” per se, ale analiza prawdopodobieństwa wciąż może być odkrywcza.
whuber
@ Whuber, zrobiłem to, ale wspomniałem o tym tylko mimochodem. (Optymalną transformacją Box-Coxa jest odwrotny pierwiastek kwadratowy.) Otrzymujesz ten sam wynik, ale wartości p są (nadal bardzo wysokie, ale) mniej znaczące.
Gung - Przywróć Monikę
3
Bardzo podoba mi się pomysł, że powinieneś modelować to, co Twoim zdaniem jest procesem generowania. Mój problem polega na tym, że nawet gdy mieszanki Gaussa działają dobrze, uważam, że powinna istnieć merytoryczna interpretacja. Gdyby OP powiedział nam więcej nawet o tym, jakie dane są lepsze, możliwe są domysły.
Nick Cox,
10

Kontynuując pomysły w odpowiedzi i komentarzach @ Nicka, możesz zobaczyć, jak szerokie musi być pasmo, aby po prostu spłaszczyć tryb dodatkowy:

wprowadź opis zdjęcia tutaj

Weź to oszacowanie gęstości jądra jako proksymalną wartość zerową - rozkład najbliższy danym, ale wciąż zgodny z hipotezą zerową, że jest to próbka z populacji nieimodalnej - i symuluj na jej podstawie. W symulowanych próbkach tryb wtórny często nie wygląda tak wyraźnie i nie trzeba poszerzać szerokości pasma, aby go spłaszczyć.

<code> wprowadź opis obrazu tutaj </code>

Sformalizowanie tego podejścia prowadzi do testu przeprowadzonego przez Silvermana (1981), „Używanie oszacowań gęstości jądra do badania modalności”, JRSS B , 43 , 1. silvermantestPakiet Schwaiger i Holzmann realizuje ten test, a także procedurę kalibracji opisaną przez Hall & York ( 2001), „O kalibracji testu Silvermana dla multimodalności”, Statistica Sinica , 11 , str. 515, która dostosowuje się do asymptotycznego konserwatyzmu. Wykonanie testu na danych z zerową hipotezą nieimodalności daje wartości p 0,08 bez kalibracji i 0,02 z kalibracją. Nie jestem wystarczająco zaznajomiony z testem zanurzeniowym, aby odgadnąć, dlaczego może się różnić.

Kod R:

  # kernel density estimate for x using Sheather-Jones method to estimate b/w:
density(x, kernel="gaussian", bw="SJ") -> dens.SJ
  # tweak b/w until mode just disappears:
density(x, kernel="gaussian", bw=3160) -> prox.null
  # fill matrix with simulated samples from the proximal null:
x.sim <- matrix(NA, nrow=length(x), ncol=10)
for (i in 1:10){
  x.sim[ ,i] <- rnorm(length(x), sample(x, size=length(x), replace=T), prox.null$bw)
}
  # perform Silverman test without Hall-York calibration:
require(silvermantest)
silverman.test(x, k=1, M=10000, adjust=F)
  # perform Silverman test with Hall-York calibration:
silverman.test(x, k=1, M=10000, adjust=T)
Scortchi - Przywróć Monikę
źródło
+1. Ciekawy! Jakie jądro jest tutaj używane? Jak słabo pamiętam, istnieją subtelne powody, dla których jądra gaussowskie powinny być używane w formalnych wariantach tego podejścia.
Nick Cox,
@Nick: Jądro Gaussa, ale nie pamiętam, czy jest ku temu ważny powód. Każda symulowana próbka jest przeskalowywana i istnieje korekta konserwatywnego błędu systematycznego oryginalnego testu - opracowanego przez kogoś zwanego Storey.
Scortchi - Przywróć Monikę
@NickCox: Przepraszam, wcale nie Storey.
Scortchi - Przywróć Monikę
@Scortchi, trochę poprawiłem twój tekst i kod. Mam nadzieję, że nie masz nic przeciwko. +1. Używasz również przerażającego operatora przypisywania prawej strzałki ?! Och, ludzkości ...
gung - Przywróć Monikę
2
Tak naprawdę nie jest ani lepsza, ani gorsza, ale konwencja w programowaniu polega na podawaniu zmiennych po lewej stronie i przypisywaniu im tych po prawej. Wielu ludzi jest przerażonych ->; Jestem po prostu oszołomiony.
Gung - Przywróć Monikę
7

Rzeczy, o które należy się martwić to:

  1. Rozmiar zestawu danych. Nie jest mały, nie jest duży.

  2. Zależność tego, co widzisz, od początku histogramu i szerokości pojemnika. Mając tylko jeden oczywisty wybór, ty (i my) nie mamy pojęcia o wrażliwości.

  3. Zależność tego, co widzisz, od rodzaju i szerokości jądra oraz od wszelkich innych wyborów dokonanych dla ciebie w szacowaniu gęstości. Mając tylko jeden oczywisty wybór, ty (i my) nie mamy pojęcia o wrażliwości.

W innym miejscu wstępnie zasugerowałem, że wiarygodność trybów jest wspierana (ale nie ustalona) przez merytoryczną interpretację i przez zdolność rozróżniania tej samej modalności w innych zestawach danych o tym samym rozmiarze. (Im większy, tym lepiej ....)

Nie możemy komentować żadnego z tych tutaj. Jednym małym uchwytem dotyczącym powtarzalności jest porównanie tego, co otrzymujesz z próbkami bootstrap tego samego rozmiaru. Oto wyniki eksperymentu z użyciem tokena przy użyciu Staty, ale to, co widzisz, jest arbitralnie ograniczone do domyślnych ustawień Staty, które same są udokumentowane jako wyrwane z powietrza . Otrzymałem oszacowania gęstości dla oryginalnych danych i dla 24 próbek bootstrap z tego samego.

Wskazanie (nie więcej, nie mniej) jest tym, co myślę, że doświadczeni analitycy zgadną w jakikolwiek sposób z twojego wykresu. Tryb lewej ręki jest wysoce powtarzalny, a prawa ręka jest wyraźnie bardziej delikatna.

Zauważ, że jest to nieuniknione: ponieważ jest mniej danych w pobliżu trybu po prawej stronie, nie zawsze pojawią się one ponownie w próbce ładowania początkowego. Ale to także kluczowy punkt.

wprowadź opis zdjęcia tutaj

Należy zauważyć, że punkt 3. powyżej pozostaje nietknięty. Ale wyniki są gdzieś pomiędzy unimodal i bimodal.

Dla zainteresowanych jest to kod:

clear 
set scheme s1color 
set seed 2803 

mat data = (10346, 13698, 13894, 19854, 28066, 26620, 27066, 16658, 9221, 13578, 11483, 10390, 11126, 13487, 15851, 16116, 24102, 30892, 25081, 14067, 10433, 15591, 8639, 10345, 10639, 15796, 14507, 21289, 25444, 26149, 23612, 19671, 12447, 13535, 10667, 11255, 8442, 11546, 15958, 21058, 28088, 23827, 30707, 19653, 12791, 13463, 11465, 12326, 12277, 12769, 18341, 19140, 24590, 28277, 22694, 15489, 11070, 11002, 11579, 9834, 9364, 15128, 15147, 18499, 25134, 32116, 24475, 21952, 10272, 15404, 13079, 10633, 10761, 13714, 16073, 23335, 29822, 26800, 31489, 19780, 12238, 15318, 9646, 11786, 10906, 13056, 17599, 22524, 25057, 28809, 27880, 19912, 12319, 18240, 11934, 10290, 11304, 16092, 15911, 24671, 31081, 27716, 25388, 22665, 10603, 14409, 10736, 9651, 12533, 17546, 16863, 23598, 25867, 31774, 24216, 20448, 12548, 15129, 11687, 11581)
set obs `=colsof(data)' 
gen data = data[1,_n] 

gen index = . 

quietly forval j = 1/24 { 
    replace index = ceil(120 * runiform()) 
    gen data`j' = data[index]
    kdensity data`j' , nograph at(data) gen(xx`j' d`j') 
} 

kdensity data, nograph at(data) gen(xx d) 

local xstuff xtitle(data/1000) xla(10000 "10" 20000 "20" 30000 "30") sort 
local ystuff ysc(r(0 .0001)) yla(none) `ystuff'   

local i = 1 
local colour "orange" 
foreach v of var d d? d?? { 
    line `v' data, lc(`colour') `xstuff'  `ystuff' name(g`i', replace) 
    local colour "gs8" 
    local G `G' g`i' 
    local ++i 
} 

graph combine `G' 
Nick Cox
źródło
+1 Podoba mi się twoje podejście do ładowania: tablica wykresów pomaga każdemu lepiej zrozumieć dane. Zastanawiam się, czy te wykresy mogą być wrażliwe na to, jak Stata szacuje przepustowość. Podejrzewam, że może to doprowadzić do niedostatecznego zasilania, ponieważ jego szacunki prawdopodobnie oparte są na niejednoznacznym założeniu, co prowadzi do stosunkowo szerokiego pasma. Nawet nieco węższe oszacowanie przepustowości może sprawić, że drugi tryb będzie bardziej widoczny we wszystkich próbkach bootstrap.
whuber
2
@whuber Thanks! Jak zwykle skupiasz się bezbłędnie na słabościach, o które musimy się martwić, i zgadzam się. Wraz ze wzrostem przepustowości jądra pojawienie się niejednoznaczności ma tendencję do nieuchronności. I odwrotnie, małe szerokości pasma często po prostu wskazują fałszywe tryby, które są niepowtarzalne i / lub trywialne. Kompromis jest naprawdę delikatny. Myślę, że główną zaletą tego podejścia jest retoryka: „Czy da się to powtórzyć, jeśli będziemy drżeć?” Często martwię się o gotowość użytkowników oprogramowania do kopiowania domyślnych wyników bez refleksji.
Nick Cox,
2
Istnieją systematyczne podejścia do tego problemu polegające na stopniowym modyfikowaniu przepustowości oraz śledzeniu wyglądu i znikania trybów w zależności od przepustowości. Zasadniczo tryb wiarygodny utrzymuje się, a tryb mniej wiarygodny nie. To urocze podejście, ale czasami podpala się konstruktora tunelu, gdy zrobi to szpadel. Na przykład, jeśli przekręcisz opcje histogramu, a tryb dodatkowy znika (lub porusza się) zbyt łatwo, nie wierz w to.
Nick Cox,
2

Identyfikacja trybu nieparametrycznego LP

LP Nonparametric Mode Identification (nazwa algorytmu LPMode , odnośnik do artykułu podano poniżej)

Tryby MaxEnt [trójkąty w kolorze czerwonym na wykresie]: 12783.36 i 24654.28.

Tryby L2 [zielone trójkąty na wykresie]: 13054,70 i 24111,61.

Interesujące jest zwrócenie uwagi na kształty modalne, zwłaszcza ten drugi, który wykazuje znaczną skośność (tradycyjny model mieszanki gaussowskiej prawdopodobnie tutaj zawodzi).

Mukhopadhyay, S. (2016) Identyfikacja w trybie wielkoskalowym i nauki oparte na danych. https://arxiv.org/abs/1509.06428

Deep Mukherjee
źródło
1
Czy potrafisz opracować i podać kontekst do wprowadzenia i wyjaśnienia tych metod? Miło jest mieć link do gazety, ale wolimy, aby nasze odpowiedzi były niezależne, szczególnie jeśli link znika.
Gung - Przywróć Monikę
Kontekstem jest pierwotne pytanie: czy istnieje multimodalność? jeśli tak, to pozycje. istotność nowej metody wynika z faktu, że polowanie na nierówności w sposób nieparametryczny jest trudnym problemem modelowania.
Deep Mukherjee
@ gung prosi o rozwinięcie odpowiedzi. Na przykład wynik pochodzi z metody opisanej w artykule, który nie ma publicznej wersji.
Nick Cox
2
Nie, mam na myśli, co to jest „LP Nonparametric Mode Identification”? Co to jest „MaxEnt”? Itd. W kilku zdaniach, jak to działa? Dlaczego / kiedy może być lepsza od innych metod? Itd. Zdaję sobie sprawę, że link do artykułu, który je wyjaśnia, ale byłoby miło mieć kilka zdań, aby je tutaj przedstawić, szczególnie jeśli link znika, ale nawet jeśli nie, aby dać przyszłym czytelnikom poczucie, czy chcę zastosować tę metodę.
Gung - Przywróć Monikę
2
@DeepMukherjee, na pewno nie musisz przepisywać całego artykułu w swoim poście. Wystarczy dodać kilka zdań mówiących, co to jest i jak to działa.
Gung - Przywróć Monikę