¿Cómo comprobar si mi distribución es multimodal?

Question

Cuando trazo un histograma de mis datos, tiene dos picos:

¿Significa eso una posible distribución multimodal? He comprobado el dip.test en R ( library(diptest) ), y la salida es:

D = 0.0275, p-value = 0.7913

¿Puedo concluir que mis datos tienen una distribución multimodal?

DATOS

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

Answer 1

@NickCox ha presentado una estrategia interesante (+1). Sin embargo, podría considerarla más bien de carácter exploratorio, debido a la preocupación que @whuber señala .

Permítanme sugerir otra estrategia: Podría ajustar un modelo de mezcla finita gaussiana. Tenga en cuenta que esto hace la suposición muy fuerte que sus datos se extraen de uno o más verdaderos normales. Como señalan tanto @whuber como @NickCox en los comentarios, sin una interpretación sustantiva de estos datos -apoyada por una teoría bien establecida- que respalde esta suposición, esta estrategia debería considerarse también exploratoria.

En primer lugar, sigamos la sugerencia de @Glen_b y analicemos los datos utilizando el doble de bins:

Seguimos viendo dos modos; en todo caso, aquí se ven más claramente. (Obsérvese también que la línea de densidad del núcleo debería ser idéntica, pero aparece más extendida debido al mayor número de intervalos).

Ahora vamos a ajustar un modelo de mezcla finita gaussiana. En R puede utilizar el Mclust paquete para hacerlo:

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 

Dos componentes normales optimizan el BIC. Para comparar, podemos forzar un ajuste de un componente y realizar una prueba de razón de verosimilitud:

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

Esto sugiere que es extremadamente improbable que se encuentren datos tan alejados de la unimodalidad como los suyos si provienen de una única distribución normal verdadera.

Algunas personas no se sienten cómodas utilizando una prueba paramétrica en este caso (aunque si los supuestos se mantienen, no conozco ningún problema). Una técnica muy aplicable es utilizar el método de ajuste cruzado de Bootstrap paramétrico (describo el algoritmo ici ). Podemos intentar aplicarlo a estos datos:

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 

Los estadísticos de resumen y los gráficos de densidad del núcleo para las distribuciones de muestreo muestran varias características interesantes. La probabilidad logarítmica para el modelo de un solo componente rara vez es mayor que la del ajuste de dos componentes, incluso cuando el verdadero proceso de generación de datos tiene un solo componente, y cuando es mayor, la cantidad es trivial. La idea de comparar modelos que difieren en su capacidad de ajuste de los datos es una de las motivaciones del PBCM. Las dos distribuciones de muestreo apenas se solapan; sólo el 0,35% de x2.d son menores que el máximo x1.d valor. Si selecciona un modelo de dos componentes si la diferencia en la probabilidad logarítmica es >9,7, seleccionará incorrectamente el modelo de un componente en un 0,01% y el modelo de dos componentes en un 0,02% de las veces. Estos son altamente discriminables. Si, por otro lado, usted eligió utilizar el modelo de un componente como hipótesis nula, su resultado observado es lo suficientemente pequeño como para no aparecer en la distribución de muestreo empírica en 10.000 iteraciones. Podemos utilizar la regla de 3 (véase ici ) para poner un límite superior al valor p, es decir, estimamos que su valor p es inferior a 0,0003. Es decir, es altamente significativo.

Esto plantea la cuestión de por qué estos resultados divergen tanto de su prueba de inmersión. (Para responder a su pregunta explícita, su prueba de inmersión no proporciona ninguna prueba de que haya dos modos reales). Sinceramente, no conozco la prueba de inmersión, así que es difícil de decir; puede que no tenga la potencia suficiente. Sin embargo, creo que la respuesta más probable es que este enfoque asume que sus datos son generados por verdaderas normales. Una prueba de Shapiro-Wilk para tus datos es altamente significativa ( $p < .000001$ ), y también es muy significativo para la transformación óptima de Box-Cox de sus datos (la raíz cuadrada inversa; $p < .001$ ). Sin embargo, los datos nunca son realmente normales (cf, esta famosa cita ), y los componentes subyacentes, en caso de que existan, tampoco se garantiza que sean perfectamente normales. Si te parece razonable que tus datos puedan provenir de una distribución sesgada positivamente, en lugar de una normal, este nivel de bimodalidad bien puede estar dentro del rango típico de variación, que es lo que sospecho que dice la prueba de inmersión.

Answer 2

Siguiendo con las ideas de la respuesta de @Nick y los comentarios, se puede ver la amplitud que debe tener el ancho de banda para sólo aplanar el modo secundario:

Tome esta estimación de la densidad del kernel como la nula proximal -la distribución más cercana a los datos, pero aún consistente con la hipótesis nula de que es una muestra de una población unimodal- y simule a partir de ella. En las muestras simuladas, el modo secundario no suele ser tan distinto, y no es necesario ampliar tanto el ancho de banda para aplanarlo.

La formalización de este enfoque conduce a la prueba dada en Silverman (1981), "Using kernel density estimates to investigate modality", JRSS B , 43 1. Schwaiger & Holzmann's silvermantest paquete implementa esta prueba, y también el procedimiento de calibración descrito por Hall & York (2001), "On the calibration of Silverman's test for multimodality", Statistica Sinica , 11 , p 515, que ajusta el conservadurismo asintótico. La realización de la prueba en sus datos con una hipótesis nula de unimodalidad da lugar a valores p de 0,08 sin calibración y 0,02 con calibración. No estoy lo suficientemente familiarizado con la prueba de inmersión para adivinar por qué podría diferir.

Código 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=TRUE), 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)

Answer 3

Las cosas de las que hay que preocuparse son:

1. El tamaño del conjunto de datos. No es pequeño, ni grande.

2. La dependencia de lo que se ve en el origen del histograma y el ancho de la bandeja. Con sólo una opción evidente, usted (y nosotros) no tiene idea de la sensibilidad.

3. La dependencia de lo que se ve en el tipo y la anchura del núcleo y cualquier otra elección que se haga para la estimación de la densidad. Con una sola elección evidente, usted (y nosotros) no tiene idea de la sensibilidad.

En otro lugar he sugerido provisionalmente que la credibilidad de las modalidades se apoya (pero no se establece) en una interpretación sustantiva y en la capacidad de discernir la misma modalidad en otros conjuntos de datos del mismo tamaño. (Más grande es mejor también....)

No podemos comentar ninguno de ellos aquí. Un pequeño asidero en la repetibilidad es comparar lo que se obtiene con muestras bootstrap del mismo tamaño. Aquí están los resultados de un experimento simbólico utilizando Stata, pero lo que se ve se limita arbitrariamente a los valores por defecto de Stata, que a su vez están documentados como sacados del aire . He obtenido estimaciones de densidad para los datos originales y para 24 muestras bootstrap de los mismos.

La indicación (ni más ni menos) es la que creo que los analistas experimentados adivinarían de cualquier manera a partir de su gráfico. El modo de la izquierda es muy repetible y el de la derecha es claramente más frágil.

Tenga en cuenta que esto es inevitable: como hay menos datos cerca de la moda de la derecha, no siempre reaparecerá en una muestra bootstrap. Pero este es también el punto clave.

Obsérvese que el punto 3. anterior permanece intacto. Pero los resultados están entre unimodal y bimodal.

Para los interesados, este es el código:

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'

Answer 4

Identificación del modo no paramétrico del LP (nombre del algoritmo LPMode (la referencia del documento se indica a continuación)

MaxEnt Modos [Triángulos de color rojo en el gráfico]: 12783,36 y 24654,28.

L2 Modos [Triángulos de color verde en el gráfico]: 13054,70 y 24111,61.

Es interesante observar las formas modales, especialmente la segunda, que muestra una considerable asimetría (Es probable que el modelo tradicional de mezcla gaussiana falle aquí).

Mukhopadhyay, S. (2016) Identificación de modos a gran escala y ciencias basadas en datos. https://arxiv.org/abs/1509.06428

Answer 5

Añadiré que puedes sustituir las barras de progreso, que no se iluminan mucho, por un recuento de iteraciones con esto:

    x2.d[i] = Mclust(x2, G=2, verbose=FALSE)$loglik - Mclust(x2, 
                     G=1, verbose=FALSE)$loglik
    x1.d[i] = Mclust(x1, G=2, verbose=FALSE)$loglik - Mclust(x1, 
                     G=1, verbose=FALSE)$loglik
    cat(sprintf("\rIteration %s of %s (%.1f%% complete)", i, B, 
              i/B*100))