Wie teste ich, ob meine Distribution multimodal ist?

21

Wenn ich ein Histogramm meiner Daten zeichne, hat es zwei Peaks:

Histogramm

Bedeutet das eine mögliche multimodale Verteilung? Ich habe das dip.testin R ( library(diptest)) ausgeführt und die Ausgabe ist:

D = 0.0275, p-value = 0.7913

Kann ich daraus schließen, dass meine Daten eine multimodale Verteilung haben?

DATEN

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
user1260391
quelle
3
Verwenden Sie mehr Behälter in Ihrem Histogramm. Ich schlage etwa doppelt so viele vor
Glen_b -Reinstate Monica
1
In dieser Antwort werden neun verschiedene Tests erwähnt , von denen einige für Ihre Situation relevant sein können.
Glen_b -Reinstate Monica
1
Dieses Papier ist wahrscheinlich nützlich für Sie, wenn Sie es noch nicht gesehen haben (auch dieses Follow-up )
Eoin

Antworten:

15

@NickCox hat eine interessante Strategie vorgestellt (+1). Ich könnte es für mehr explorativen in der Natur jedoch aufgrund der Sorge , dass @whuber weist darauf hin .

Lassen Sie mich eine andere Strategie vorschlagen: Sie könnten ein Gaußsches Modell für endliche Gemische verwenden. Beachten Sie, dass dies zu der starken Annahme führt, dass Ihre Daten aus einer oder mehreren echten Normalen stammen. Wie sowohl @whuber als auch @NickCox in den Kommentaren darlegen, sollte diese Strategie auch als explorativ betrachtet werden, ohne dass eine inhaltliche Interpretation dieser Daten - unterstützt durch eine gut etablierte Theorie - diese Annahme stützt.

Folgen wir zunächst dem Vorschlag von @ Glen_b und sehen uns Ihre Daten mit doppelt so vielen Behältern an:

Bildbeschreibung hier eingeben

Wir sehen immer noch zwei Modi; wenn überhaupt, kommen sie hier deutlicher durch. (Beachten Sie auch, dass die Kerneldichtelinie identisch sein sollte, jedoch aufgrund der größeren Anzahl von Fächern breiter erscheint.)

Passen wir nun ein Gaußsches Modell mit endlicher Mischung an. In Rkönnen Sie das MclustPaket verwenden, um dies zu tun:

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 

Zwei normale Komponenten optimieren den BIC. Zum Vergleich können wir eine Einkomponentenanpassung erzwingen und einen Likelihood-Ratio-Test durchführen:

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

Dies legt den Schluss nahe, dass es äußerst unwahrscheinlich ist, dass Sie Daten finden, die so unimodal sind wie Ihre, wenn sie aus einer einzigen echten Normalverteilung stammen.

Einige Leute fühlen sich mit einem parametrischen Test hier nicht wohl (obwohl ich, wenn die Annahmen stimmen, kein Problem kenne). Eine sehr allgemein anwendbare Technik ist die Verwendung der parametrischen Bootstrap-Cross-Fitting-Methode (den Algorithmus beschreibe ich hier ). Wir können versuchen, es auf diese Daten anzuwenden:

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 

Bildbeschreibung hier eingeben

Die zusammenfassende Statistik und die Kernel-Dichtediagramme für die Stichprobenverteilungen zeigen einige interessante Merkmale. Die Protokollwahrscheinlichkeit für das Einzelkomponentenmodell ist selten größer als die der Zweikomponentenanpassung, selbst wenn der Prozess zur Erzeugung echter Daten nur eine einzelne Komponente enthält und wenn sie größer ist, ist die Menge trivial. Die Idee, Modelle zu vergleichen, die sich in ihrer Fähigkeit unterscheiden, Daten anzupassen, ist eine der Beweggründe für die PBCM. Die beiden Stichprobenverteilungen überlappen sich kaum; nur .35% von x2.dsind kleiner als das Maximumx1.dWert. Wenn Sie ein Zweikomponentenmodell ausgewählt haben, bei dem der Unterschied in der Protokollwahrscheinlichkeit> 9,7 war, haben Sie in den meisten Fällen fälschlicherweise das Einkomponentenmodell 0,01% und das Zweikomponentenmodell 0,02% ausgewählt. Diese sind sehr diskriminierbar. Wenn Sie sich dagegen für die Verwendung des Einkomponentenmodells als Nullhypothese entschieden haben, ist das beobachtete Ergebnis so gering, dass es in der empirischen Stichprobenverteilung in 10.000 Iterationen nicht auftaucht. Wir können die Regel 3 (siehe hier ) verwenden, um eine Obergrenze für den p-Wert festzulegen. Wir schätzen nämlich, dass Ihr p-Wert kleiner als 0,0003 ist. Das heißt, das ist sehr wichtig.

p<.000001p<.001) und die zugrunde liegenden Komponenten, falls vorhanden, sind ebenfalls nicht unbedingt normal. Wenn Sie es für sinnvoll halten, dass Ihre Daten aus einer positiv verzerrten Verteilung stammen und nicht aus einer normalen, dann kann dieses Maß an Bimodalität durchaus innerhalb des typischen Variationsbereichs liegen, was ich im Dip-Test vermute.

gung - Wiedereinsetzung von Monica
quelle
1
Das Problem bei diesem Ansatz ist, dass die Alternative, mit der Sie eine Gaußsche Mischung vergleichen, nicht sehr vernünftig ist. Eine vernünftigere wäre, dass die Verteilung eine Art von rechtwinkliger Verteilung ist, wie z. B. eine Gamma-Verteilung. Es ist fast eine Selbstverständlichkeit, dass eine Mischung zu fast jedem verzerrten Datensatz "signifikant" besser passt als zu einem einzelnen Gaußschen Datensatz .
Whuber
Du hast recht, @whuber. Ich habe versucht, diesen Punkt explizit zu machen. Ich bin nicht sicher, wie man ein Gamma-FMM macht, aber das wäre besser.
gung - Reinstate Monica
1
Da dies explorativ ist, besteht ein Gedanke darin, zu versuchen, die ursprüngliche Verteilung in Symmetrie umzuwandeln (möglicherweise mit einer versetzten Box-Cox-Transformation, die robust aus einigen Quantilen der Daten geschätzt wird), und Ihren Ansatz erneut zu versuchen. Natürlich würden Sie nicht über "Signifikanz" per se sprechen, aber die Analyse der Wahrscheinlichkeit könnte immer noch aufschlussreich sein.
Whuber
@whuber, das habe ich gemacht, aber ich habe es nur beiläufig erwähnt. (Die optimale Box-Cox-Transformation ist die inverse Quadratwurzel.) Sie erhalten das gleiche Ergebnis, aber die p-Werte sind (immer noch hoch, aber) weniger signifikant.
gung - Reinstate Monica
3
Mir gefällt die Idee sehr, dass Sie das modellieren sollten, was Sie für den Erzeugungsprozess halten. Mein Problem ist, dass selbst wenn Gaußsche Mischungen gut funktionieren, ich der Meinung bin, dass es eine inhaltliche Interpretation geben sollte. Wenn das OP uns mehr darüber sagt, was die Daten sind, könnten einige bessere Vermutungen möglich sein.
Nick Cox
10

Im Anschluss an den Ideen in @ Nick Antwort und Kommentare können Sie sehen , wie breit der Bandbreitenbedarf sein , um nur den sekundären Modus verflachen:

Bildbeschreibung hier eingeben

Nehmen Sie diese Schätzung der Kerneldichte als proximale Null - die Verteilung, die den Daten am nächsten kommt, aber dennoch mit der Nullhypothese übereinstimmt, dass es sich um eine Stichprobe aus einer unimodalen Population handelt - und simulieren Sie daraus. In den simulierten Samples sieht der sekundäre Modus nicht so deutlich aus, und Sie müssen die Bandbreite nicht so stark vergrößern, um ihn zu reduzieren.

<code> Bildbeschreibung hier eingeben </ code>

Die Formalisierung dieses Ansatzes führt zu dem Test in Silverman (1981), "Verwendung von Kerndichteschätzungen zur Untersuchung der Modalität", JRSS B , 43 , 1. Das silvermantestPaket von Schwaiger & Holzmann implementiert diesen Test und das von Hall & York ( 2001), "Zur Kalibrierung von Silvermans Test auf Multimodalität", Statistica Sinica , 11 , S. 515, der sich auf asymptotischen Konservatismus einstellt. Wenn Sie den Test für Ihre Daten mit einer Nullhypothese der Unimodalität durchführen, erhalten Sie p-Werte von 0,08 ohne Kalibrierung und 0,02 mit Kalibrierung. Ich kenne den Dip-Test nicht gut genug, um zu erraten, warum er sich unterscheiden könnte.

R-Code:

  # 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 - Wiedereinsetzung von Monica
quelle
+1. Interessant! Welcher Kernel wird hier verwendet? Soweit ich mich recht erinnere, gibt es subtile Gründe, warum Gauß'sche Kerne für formale Varianten dieses Ansatzes verwendet werden sollten.
Nick Cox
@Nick: Gauß'scher Kernel, aber ich kann mich nicht erinnern, ob es einen zwingenden Grund dafür gibt. Jedes simulierte Sample wird neu skaliert und es gibt eine Korrektur für eine konservative Tendenz, die der ursprüngliche Test hat - erarbeitet von jemandem namens Storey, glaube ich.
Scortchi - Wiedereinsetzung von Monica
@ NickCox: Sorry, überhaupt nicht Storey.
Scortchi - Wiedereinsetzung von Monica
@ Scortchi, ich habe deinen Text & Code leicht überarbeitet. Ich hoffe es macht dir nichts aus. +1. Außerdem verwenden Sie den gefürchteten Rechtspfeil-Zuweisungsoperator ?! Oh die Menschheit ...
gung - Wiedereinsetzung von Monica
2
Eigentlich ist es nicht besser oder schlechter, aber die Konvention bei der Programmierung besteht darin, Ihre Variablen links anzugeben und ihnen rechts das zuzuweisen, was ihnen zugewiesen ist. Viele Leute sind durch ausgeflippt ->; Ich bin nur amüsiert.
gung - Wiedereinsetzung von Monica
7

Die Dinge, über die man sich Sorgen machen muss, sind:

  1. Die Größe des Datensatzes. Es ist nicht klein, nicht groß.

  2. Die Abhängigkeit dessen, was Sie sehen, von Histogrammursprung und Behälterbreite. Mit nur einer Wahl haben Sie (und wir) keine Ahnung von Sensibilität.

  3. Die Abhängigkeit dessen, was Sie sehen, von Kerneltyp und -breite und was auch immer für andere Auswahlmöglichkeiten Sie bei der Dichteschätzung treffen. Mit nur einer Wahl haben Sie (und wir) keine Ahnung von Sensibilität.

An anderer Stelle habe ich vorläufig vorgeschlagen, dass die Glaubwürdigkeit von Verkehrsträgern durch eine inhaltliche Auslegung und durch die Fähigkeit, die gleiche Modalität in anderen Datensätzen derselben Größe zu erkennen, gestützt (aber nicht nachgewiesen) wird. (Größer ist auch besser ....)

Beides können wir hier nicht kommentieren. Ein kleiner Punkt bei der Wiederholbarkeit ist der Vergleich mit Bootstrap-Beispielen derselben Größe. Hier sind die Ergebnisse eines Token-Experiments mit Stata, aber was Sie sehen, ist willkürlich auf die Standardeinstellungen von Stata beschränkt, die selbst als aus der Luft gezupft dokumentiert sind . Ich habe Dichteschätzungen für die Originaldaten und für 24 Bootstrap-Beispiele von denselben erhalten.

Der Indikator (nicht mehr und nicht weniger) ist, wie erfahrene Analysten meiner Meinung nach nur vermuten, aus Ihrer Grafik heraus. Der linke Modus ist sehr wiederholbar und der rechte Modus ist deutlich anfälliger.

Beachten Sie, dass dies unvermeidlich ist: Da sich weniger Daten in der Nähe des rechten Modus befinden, werden sie in einem Bootstrap-Beispiel nicht immer wieder angezeigt. Das ist aber auch der entscheidende Punkt.

Bildbeschreibung hier eingeben

Beachten Sie, dass Punkt 3. oben unberührt bleibt. Die Ergebnisse liegen jedoch irgendwo zwischen unimodal und bimodal.

Für Interessierte ist dies der Code:

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
quelle
+1 Ich mag Ihren Bootstrap-Ansatz: Die Anordnung der Diagramme hilft jedem, die Daten besser zu verstehen. Ich frage mich, ob diese Diagramme empfindlich dafür sind, wie Stata die Bandbreite schätzt. Ich vermute, dass dies zu einem Test mit zu geringer Leistung führen könnte, da seine Schätzung wahrscheinlich auf einer unimodalen Annahme beruht, die zu einer relativ großen Bandbreite führt. Selbst eine etwas engere Bandbreitenschätzung könnte den zweiten Modus in allen Bootstrap-Beispielen hervorheben.
whuber
2
@whuber Danke! Wie üblich konzentrieren Sie sich zielsicher auf die Schwächen, über die wir uns Sorgen machen müssen, und ich stimme zu. Wenn die Kernel-Bandbreiten zunehmen, besteht die Tendenz, dass Unimodalität unvermeidlich ist. Umgekehrt weisen kleine Bandbreiten oft nur auf unwiederholbare und / oder triviale Störungsmoden hin. Der Kompromiss ist wirklich heikel. Ich denke, der Hauptnutzen dieses Ansatzes ist die Rhetorik von "Ist das reproduzierbar, wenn wir wackeln?" Ich bin oft besorgt über die Bereitschaft von Softwarebenutzern, Standardergebnisse ohne Reflexion zu kopieren.
Nick Cox
2
Es gibt systematische Ansätze für dieses Problem, die darauf beruhen, die Bandbreite schrittweise zu ändern und das Auftreten und Verschwinden von Modi zu verfolgen, wenn sich die Bandbreite ändert. Im Wesentlichen bleibt ein glaubwürdiger Modus bestehen und ein weniger glaubwürdiger Modus nicht. Es ist eine nette Herangehensweise, aber manchmal wird ein Tunnelkonstruktor ausgelöst, wenn ein Spaten genügt. Wenn Sie beispielsweise die Auswahl des Histogramms ändern und der sekundäre Modus allzu schnell verschwindet (oder sich bewegt), glauben Sie es nicht.
Nick Cox
2

LP Nichtparametrische Modusidentifikation

LP Nonparametric Mode Identification (Name des LPMode- Algorithmus, siehe unten)

MaxEnt- Modi [Rote Farbdreiecke im Diagramm]: 12783,36 und 24654,28.

L2- Modi [Grüne Farbdreiecke im Diagramm]: 13054.70 und 24111.61.

Interessant sind die modalen Formen, insbesondere die zweite, die eine beträchtliche Schiefe aufweist (traditionelles Modell der Gaußschen Mischung, das hier wahrscheinlich versagt).

Mukhopadhyay, S. (2016) Identifizierung im großen Maßstab und datengetriebene Wissenschaften. https://arxiv.org/abs/1509.06428

Deep Mukherjee
quelle
1
Können Sie einen Kontext für die Einführung und Erklärung dieser Methoden erläutern? Es ist schön, einen Link zum Artikel zu haben, aber wir bevorzugen, dass unsere Antworten hier in sich geschlossen sind, insbesondere, wenn der Link nicht mehr funktioniert.
gung - Wiedereinsetzung von Monica
Der Kontext ist die ursprüngliche Frage: Gibt es Multimodalität? Wenn ja, Positionen. Die Relevanz einer neuen Methode ergibt sich aus der Tatsache, dass die nichtparametrische Suche nach Unebenheiten ein schwieriges Modellierungsproblem darstellt.
Deep Mukherjee
@gung bittet Sie, Ihre Antwort zu erweitern. Das Ergebnis stammt beispielsweise aus einer Methode, die in einem Artikel ohne öffentliche Version erläutert wurde.
Nick Cox
2
Nein, ich meine, was ist "LP Nonparametric Mode Identification"? Was ist "MaxEnt"? Wie funktioniert das in ein paar Sätzen? Warum / wann ist es anderen Methoden vorzuziehen? Ich bin mir bewusst, dass Sie einen Link zu dem Artikel haben, der sie erklärt, aber es wäre schön, ein paar Sätze zu haben, um sie hier vorzustellen, insbesondere, wenn der Link nicht mehr funktioniert, aber auch, um zukünftigen Lesern ein Gefühl dafür zu geben, ob sie es tun möchte diese Methode weiterverfolgen.
gung - Wiedereinsetzung von Monica
2
@DeepMukherjee, du musst auf keinen Fall das gesamte Papier in deinem Beitrag umschreiben. Fügen Sie einfach ein paar Sätze hinzu, in denen steht, was es ist und wie es funktioniert.
gung - Wiedereinsetzung von Monica