Comment tester si ma distribution est multimodale?

21

Lorsque je trace un histogramme de mes données, il a deux pics:

histogramme

Cela signifie-t-il une distribution multimodale potentielle? J'ai exécuté le dip.testdans R ( library(diptest)), et la sortie est:

D = 0.0275, p-value = 0.7913

Je peux conclure que mes données ont une distribution multimodale?

LES DONNÉES

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
la source
3
Utilisez plus de bacs dans votre histogramme. Je suggère environ deux fois plus
Glen_b -Reinstate Monica
1
Il y a mention de neuf tests différents dans cette réponse , dont certains peuvent être pertinents pour votre situation.
Glen_b -Reinstate Monica
1
Ce document est susceptible de vous être utile, si vous ne l'avez pas déjà vu (également ce suivi )
Eoin

Réponses:

15

@NickCox a présenté une stratégie intéressante (+1). Je pourrais cependant le considérer de nature plus exploratoire, en raison de la préoccupation que souligne @whuber .

Permettez-moi de suggérer une autre stratégie: vous pourriez adapter un modèle de mélange fini gaussien. Notez que cela fait l'hypothèse très forte que vos données sont tirées d'une ou plusieurs vraies normales. Comme le soulignent @whuber et @NickCox dans les commentaires, sans une interprétation substantielle de ces données - étayée par une théorie bien établie - pour étayer cette hypothèse, cette stratégie doit également être considérée comme exploratoire.

Tout d'abord, suivons la suggestion de @ Glen_b et examinons vos données en utilisant deux fois plus de casiers:

entrez la description de l'image ici

Nous voyons toujours deux modes; au contraire, ils ressortent plus clairement ici. (Notez également que la ligne de densité du noyau doit être identique, mais semble plus étalée en raison du plus grand nombre de cases.)

Permet maintenant d'ajuster un modèle de mélange fini gaussien. Dans R, vous pouvez utiliser le Mclustpackage pour ce faire:

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 

Deux composants normaux optimisent le BIC. À titre de comparaison, nous pouvons forcer un ajustement à un composant et effectuer un test de rapport de vraisemblance:

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

Cela suggère qu'il est extrêmement improbable que vous trouviez des données aussi loin d'unimodales que les vôtres si elles provenaient d'une seule distribution normale normale.

Certaines personnes ne se sentent pas à l'aise d'utiliser un test paramétrique ici (bien que si les hypothèses se vérifient, je ne connais aucun problème). Une technique très largement applicable consiste à utiliser la méthode de cross-fit paramétrique Bootstrap (je décris ici l'algorithme ). Nous pouvons essayer de l'appliquer à ces données:

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 

entrez la description de l'image ici

Les statistiques récapitulatives et les graphiques de densité du noyau pour les distributions d'échantillonnage présentent plusieurs caractéristiques intéressantes. La probabilité logarithmique pour le modèle à un seul composant est rarement supérieure à celle de l'ajustement à deux composants, même lorsque le véritable processus de génération de données n'a qu'un seul composant, et lorsqu'il est supérieur, le montant est trivial. L'idée de comparer des modèles qui diffèrent par leur capacité à ajuster les données est l'une des motivations du PBCM. Les deux distributions d'échantillonnage se chevauchent à peine; seulement .35% x2.dsont inférieurs au maximumx1.dvaleur. Si vous avez sélectionné un modèle à deux composants si la différence de probabilité logarithmique était> 9,7, vous sélectionneriez incorrectement le modèle à un composant 0,01% et le modèle à deux composants 0,02% du temps. Ce sont très discriminables. Si, en revanche, vous avez choisi d'utiliser le modèle à une composante comme hypothèse nulle, votre résultat observé est suffisamment petit pour ne pas apparaître dans la distribution d'échantillonnage empirique à 10 000 itérations. Nous pouvons utiliser la règle de 3 (voir ici ) pour placer une limite supérieure sur la valeur de p, à savoir, nous estimons que votre valeur de p est inférieure à .0003. Autrement dit, c'est très important.

p<.000001p<.001), et les composants sous-jacents, s'ils existent, ne sont pas non plus garantis comme étant parfaitement normaux. Si vous trouvez qu'il est raisonnable que vos données proviennent d'une distribution asymétrique positive, plutôt que d'une normale, ce niveau de bimodalité peut très bien se situer dans la plage de variation typique, ce que je soupçonne de dire le test d'immersion.

gung - Réintégrer Monica
la source
1
Le problème avec cette approche est que l'alternative à laquelle vous comparez un mélange gaussien n'est pas très raisonnable. Une solution plus raisonnable serait que la distribution est une sorte de distorsion à droite, comme un Gamma. Il est presque acquis qu'un mélange va s'adapter à presque n'importe quel ensemble de données asymétriques «de manière significative» mieux qu'un seul gaussien ne lui convient.
whuber
Tu as raison, @whuber. J'ai essayé de le faire explicitement. Je ne sais pas comment faire un Gamma FMM, mais ce serait mieux.
gung - Rétablir Monica
1
Étant donné qu'il s'agit d'un élément exploratoire, une idée serait de tenter de transformer la distribution d'origine en symétrie (peut-être avec une transformation Box-Cox décalée, estimée de manière robuste à partir de quelques quantiles des données) et d'essayer à nouveau votre approche. Bien sûr, vous ne parleriez pas de «signification» en soi, mais l'analyse de la probabilité pourrait encore être révélatrice.
whuber
@whuber, je l'ai fait, mais je ne l'ai mentionné qu'en passant. (La transformation Box-Cox optimale est la racine carrée inverse.) Vous obtenez le même résultat, mais les valeurs de p sont (toujours très élevées, mais) moins significatives.
gung - Rétablir Monica
3
J'aime beaucoup l'idée que vous devriez modéliser ce que vous pensez être le processus de génération. Mon problème est que même lorsque les mélanges gaussiens fonctionnent bien, je pense qu'il devrait y avoir une interprétation substantielle. Si le PO nous en disait plus, même sur ce que sont les données, de meilleures suppositions pourraient être possibles.
Nick Cox
10

Donnant suite aux idées @ réponse de Nick et les commentaires, vous pouvez voir la largeur des besoins en bande passante pour être à peu aplatir le mode secondaire:

entrez la description de l'image ici

Prenez cette estimation de la densité du noyau comme le zéro proximal - la distribution la plus proche des données mais toujours cohérente avec l'hypothèse nulle qu'il s'agit d'un échantillon d'une population unimodale - et simulez-la. Dans les échantillons simulés, le mode secondaire n'a pas souvent l'air si distinct, et vous n'avez pas besoin d'élargir autant la bande passante pour l'aplatir.

<code> entrez la description de l'image ici </code>

La formalisation de cette approche conduit au test donné dans Silverman (1981), "Using kernel density estimations to investigation modality ", JRSS B , 43 , 1. Le silvermantestpackage de Schwaiger & Holzmann met en œuvre ce test, ainsi que la procédure d'étalonnage décrite par Hall & York ( 2001), "Sur l'étalonnage du test de Silverman pour la multimodalité", Statistica Sinica , 11 , p 515, qui ajuste pour le conservatisme asymptotique. La réalisation du test sur vos données avec une hypothèse nulle d'unimodalité donne des valeurs de p de 0,08 sans étalonnage et de 0,02 avec étalonnage. Je ne connais pas assez bien le test d'immersion pour deviner pourquoi il pourrait différer.

Code 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 - Réintégrer Monica
la source
+1. Intéressant! Quel noyau est utilisé ici? Comme je me souviens vaguement, il existe des raisons subtiles pour lesquelles les noyaux gaussiens devraient être utilisés pour des variantes formelles de cette approche.
Nick Cox
@ Nick: noyau gaussien, mais je ne me souviens pas s'il y a une raison impérieuse à cela. Chaque échantillon simulé est redimensionné, et il y a une correction pour un biais conservateur que le test original a - élaboré par quelqu'un appelé Storey je pense.
Scortchi - Réintégrer Monica
@NickCox: Désolé, pas du tout Storey.
Scortchi - Réintégrer Monica
@Scortchi, j'ai légèrement modifié votre texte et votre code. J'espère que ça ne vous dérange pas. +1. En outre, vous utilisez l'opérateur d'affectation de flèche droite redoutée?! Oh l'humanité ...
gung - Réinstalle Monica
2
Ce n'est ni mieux ni pire, vraiment, mais la convention en programmation est de dire vos variables à gauche et d'avoir ce qui leur est assigné à droite. Beaucoup de gens sont flippés par ->; Je suis juste perplexe.
gung - Rétablir Monica
7

Les choses à craindre incluent:

  1. La taille de l'ensemble de données. Ce n'est ni minuscule, ni grand.

  2. La dépendance de ce que vous voyez sur l'origine de l'histogramme et la largeur du bac. Avec un seul choix évident, vous (et nous) n'avons aucune idée de la sensibilité.

  3. La dépendance de ce que vous voyez sur le type et la largeur du noyau et tout autre choix est fait pour vous dans l'estimation de la densité. Avec un seul choix évident, vous (et nous) n'avons aucune idée de la sensibilité.

Ailleurs, j'ai suggéré à titre provisoire que la crédibilité des modes est soutenue (mais non établie) par une interprétation substantielle et par la capacité de discerner la même modalité dans d'autres ensembles de données de la même taille. (Plus c'est gros, mieux c'est aussi ....)

Nous ne pouvons commenter ni l'un ni l'autre ici. Une petite poignée sur la répétabilité consiste à comparer ce que vous obtenez avec des échantillons de bootstrap de la même taille. Voici les résultats d'une expérience symbolique utilisant Stata, mais ce que vous voyez est arbitrairement limité aux valeurs par défaut de Stata, qui sont elles-mêmes documentées comme arrachées . J'ai obtenu des estimations de densité pour les données originales et pour 24 échantillons de bootstrap de la même chose.

L'indication (ni plus, ni moins) est ce que je pense que les analystes expérimentés devineraient tout simplement à partir de votre graphique. Le mode gauche est hautement reproductible et la main droite est nettement plus fragile.

Notez qu'il y a une fatalité à ce sujet: comme il y a moins de données plus près du mode droit, elles ne réapparaîtront pas toujours dans un échantillon bootstrap. Mais c'est aussi le point clé.

entrez la description de l'image ici

Notez que le point 3. ci-dessus reste intact. Mais les résultats se situent quelque part entre unimodal et bimodal.

Pour les personnes intéressées, voici le 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
la source
+1 J'aime votre approche bootstrap: le tableau de tracés aide tout le monde à mieux comprendre les données. Je me demande si ces graphiques pourraient être sensibles à la façon dont Stata estime la bande passante. Je soupçonne que cela pourrait entraîner un test sous-alimenté car son estimation est probablement basée sur une hypothèse unimodale, conduisant à une bande passante relativement large. Même une estimation de bande passante légèrement plus étroite pourrait rendre le deuxième mode plus important dans tous les échantillons de bootstrap.
whuber
2
@whuber Merci! Comme d'habitude, vous vous concentrez infailliblement sur les faiblesses dont nous devons nous inquiéter, et je suis d'accord. À mesure que les bandes passantes du noyau augmentent, l'apparition d'unimodalité tend à devenir inévitable. Inversement, de petites bandes passantes indiquent souvent des modes parasites qui ne peuvent pas être répétés et / ou triviaux. Le compromis est vraiment délicat. Je pense que le principal mérite de cette approche est la rhétorique de "Est-ce reproductible si nous trémoussons?" Je suis souvent préoccupé par la volonté des utilisateurs de logiciels de copier les résultats par défaut sans réflexion.
Nick Cox
2
Il existe des approches systématiques à ce problème basées sur la modification progressive de la bande passante et le traçage de l'apparition et de la disparition des modes à mesure que la bande passante varie. Essentiellement, un mode crédible persiste et pas un mode moins crédible. C'est une approche mignonne, mais il faudrait parfois lancer un constructeur de tunnel quand une pelle fera l'affaire. Par exemple, si vous tournez les choix d'histogramme et que le mode secondaire disparaît (ou se déplace) trop facilement, ne le croyez pas.
Nick Cox
2

Identification du mode non paramétrique LP

Identification du mode non paramétrique LP (nom de l'algorithme LPMode , la référence du papier est donnée ci-dessous)

Modes MaxEnt [Triangles de couleur rouge dans l'intrigue]: 12783.36 et 24654.28.

Modes L2 [triangles de couleur verte dans l'intrigue]: 13054.70 et 24111.61.

Il est intéressant de noter les formes modales, en particulier la seconde qui présente une asymétrie considérable (modèle de mélange gaussien traditionnel susceptible d'échouer ici).

Mukhopadhyay, S. (2016) Identification des modes à grande échelle et sciences des données. https://arxiv.org/abs/1509.06428

Deep Mukherjee
la source
1
Pouvez-vous élaborer et fournir un contexte pour introduire et expliquer ces méthodes? C'est bien d'avoir un lien vers le document, mais nous préférons que nos réponses ici soient autonomes, surtout si le lien devient mort.
gung - Rétablir Monica
Le contexte est la question d'origine: existe-t-il une multimodalité? si tel est le cas. et la pertinence d'une nouvelle méthode vient du fait que la chasse aux bosses de manière non paramétrique est un problème de modélisation difficile.
Deep Mukherjee
@gung vous demande d'élargir votre réponse. Par exemple, le résultat provient d'une méthode expliquée dans un document qui n'a pas de version publique.
Nick Cox
2
Non, je veux dire ce qu'est "Identification du mode non paramétrique LP"? Qu'est-ce que "MaxEnt"? Etc. En quelques phrases, comment ça marche? Pourquoi / quand pourrait-il être préférable à d'autres méthodes? Etc. Je sais que vous avez un lien vers le document qui les explique, mais ce serait bien d'avoir quelques phrases pour les présenter ici, surtout si le lien disparaît, mais même si ce n'est pas pour donner aux futurs lecteurs une idée s'ils veulent poursuivre cette méthode.
gung - Rétablir Monica
2
@DeepMukherjee, vous n'avez certainement pas besoin de réécrire l'intégralité du document dans votre message. Ajoutez simplement quelques phrases expliquant ce que c'est et comment cela fonctionne.
gung - Rétablir Monica