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

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

¿Eso significa una posible distribución multimodal? Ejecuté dip.testen 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

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

Permítanme sugerir otra estrategia: podría ajustarse a un modelo de mezcla finita gaussiana. Tenga en cuenta que esto supone que se supone que sus datos provienen de una o más normales verdaderas. Como tanto @whuber como @NickCox señalan en los comentarios, sin una interpretación sustantiva de estos datos, respaldada por una teoría bien establecida, para apoyar esta suposición, esta estrategia también debe considerarse exploratoria.

Primero, sigamos la sugerencia de @ Glen_b y miremos sus datos usando el doble de contenedores:

Todavía vemos dos modos; en todo caso, aparecen más claramente aquí. (Tenga en cuenta también que la línea de densidad del núcleo debe ser idéntica, pero parece más extendida debido a la mayor cantidad de contenedores).

Ahora ajustemos un modelo de mezcla finita gaussiana. En R, puede usar el Mclustpaquete para hacer esto:

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. A modo de comparación, podemos forzar un ajuste de un componente y realizar una prueba de razón de probabilidad:

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 encuentre datos tan distantes de ser unimodales como los suyos si provienen de una única distribución normal verdadera.

Algunas personas no se sienten cómodas usando una prueba paramétrica aquí (aunque si los supuestos se mantienen, no conozco ningún problema). Una técnica muy ampliamente aplicable es utilizar el método de ajuste cruzado de Bootstrap paramétrico (describo el algoritmo aquí ). 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 

Las estadísticas de resumen y las gráficas de densidad del núcleo para las distribuciones de muestreo muestran varias características interesantes. La probabilidad de registro 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 para ajustar datos es una de las motivaciones detrás de la PBCM. Las dos distribuciones de muestreo apenas se superponen en absoluto; solo el .35% x2.dson inferiores al máximox1.dvalor. Si seleccionó un modelo de dos componentes si la diferencia en la probabilidad de registro fuera> 9.7, seleccionaría incorrectamente el modelo de un componente .01% y el modelo de dos componentes .02% del tiempo. Estos son altamente discriminables. Si, por otro lado, elige usar el modelo de un componente como una hipótesis nula, el resultado observado es lo suficientemente pequeño como para no aparecer en la distribución de muestreo empírico en 10,000 iteraciones. Podemos usar la regla de 3 (ver aquí ) para colocar un límite superior en el valor p, es decir, estimamos que su valor p es menor que .0003. Es decir, esto es muy significativo.

$p < .000001$$p < .001$), y los componentes subyacentes, si existen, tampoco se garantiza que sean perfectamente normales. Si considera razonable que sus datos puedan provenir de una distribución sesgada positivamente, en lugar de una normal, este nivel de bimodalidad puede estar dentro del rango típico de variación, que es lo que sospecho que dice la prueba de inmersión.

El seguimiento de las ideas en @ respuesta y los comentarios de Nick, se puede ver cómo las necesidades de ancho de banda amplia que sea a solo aplanar el modo secundario:

Tome esta estimación de la densidad del núcleo como nula proximal (la distribución más cercana a los datos, pero aún así es 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 a menudo no se ve tan distinto, y no es necesario ampliar tanto el ancho de banda para aplanarlo.

La formalización de este enfoque lleva a la prueba dada en Silverman (1981), "Uso de estimaciones de densidad del núcleo para investigar la modalidad", JRSS B , 43 , 1. El silvermantestpaquete de Schwaiger & Holzmann implementa esta prueba, y también el procedimiento de calibración descrito por Hall & York ( 2001), "Sobre la calibración de la prueba de Silverman para multimodalidad", Statistica Sinica , 11 , p 515, que se ajusta al conservadurismo asintótico. Realizar la prueba en sus datos con una hipótesis nula de resultados de unimodalidad en 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é puede ser diferente.

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

Las cosas de las que preocuparse incluyen:

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

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

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

En otras partes he sugerido tentativamente que la credibilidad de los modos está respaldada (pero no establecida) por una interpretación sustantiva y por 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 sobre ninguno de los que están aquí. Un pequeño control sobre la repetibilidad es comparar lo que obtienes con muestras de bootstrap del mismo tamaño. Estos son los resultados de un experimento simbólico con Stata, pero lo que ves está limitado arbitrariamente a los valores predeterminados de Stata, que están documentados como extraídos del aire . Obtuve estimaciones de densidad para los datos originales y para 24 muestras de bootstrap de la misma.

La indicación (ni más ni menos) es lo que creo que los analistas experimentados podrían adivinar de cualquier manera de su gráfico. El modo izquierdo es altamente repetible y el derecho es claramente más frágil.

Tenga en cuenta que esto es inevitable: como hay menos datos más cerca del modo de la derecha, no siempre volverá a aparecer en una muestra de arranque. Pero este es también el punto clave.

Tenga en cuenta que el punto 3. anterior permanece intacto. Pero los resultados están en algún lugar 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' 

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

Modos MaxEnt [triángulos de color rojo en la trama]: 12783.36 y 24654.28.

Modos L2 [triángulos de color verde en la trama]: 13054.70 y 24111.61.

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

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