¿Cómo determinar qué distribución se adapta mejor a mis datos?

> mydata [1] 37.50 46.79 48.30 46.04 43.40 39.25 38.49 49.51 40.38 36.98 40.00 [12] 38.49 37.74 47.92 44.53 44.91 44.91 40.00 41.51 47.92 36.98 43.40 [23] 42.26 41.89 38.87 43.02 39.25 40.38 42.64 36.98 44.15 44.91 43.40 [34] 49.81 38.87 40.00 52.45 53.13 47.92 52.45 44.91 29.54 27.13 35.60 [45] 45.34 43.37 54.15 42.77 42.88 44.26 27.14 39.31 24.80 16.62 30.30 [56] 36.39 28.60 28.53 35.84 31.10 34.55 52.65 48.81 43.42 52.49 38.00 [67] 38.65 34.54 37.70 38.11 43.05 29.95 32.48 24.63 35.33 41.34 # estimate shape and scale to perform KS-test for weibull distribution > fitdistr(mydata, "weibull") shape scale 6.4632971 43.2474500 ( 0.5800149) ( 0.8073102) # KS-test for weibull distribution > ks.test(mydata, "pweibull", scale=43.2474500, shape=6.4632971) One-sample Kolmogorov-Smirnov test data: mydata D = 0.0686, p-value = 0.8669 alternative hypothesis: two-sided # KS-test for normal distribution > ks.test(mydata, "pnorm", mean=mean(mydata), sd=sd(mydata)) One-sample Kolmogorov-Smirnov test data: mydata D = 0.0912, p-value = 0.5522 alternative hypothesis: two-sided 

fits <- list(no = fitdistr(mydata, "normal"), we = fitdistr(mydata, "weibull")) sapply(fits, function(i) i$loglik) no we -259.6540 -257.9268 

Primero, aquí hay algunos comentarios rápidos:

• Los $ p $ -valores de un Kolmovorov -Smirnov-Test (KS-Test) con parámetros estimados será bastante incorrecto. Desafortunadamente, no puede simplemente ajustar una distribución y luego usar los parámetros estimados en una prueba de Kolmogorov-Smirnov para probar su muestra.
• Su muestra nunca seguirá una distribución exactamente. Por lo tanto, incluso si sus $ p $ -valores de la prueba KS fueran válidos y $ > 0.05 $ , solo significaría que no puede descartar que sus datos sigan esta distribución específica. Otra formulación sería que su muestra sea compatible con una determinada distribución. Pero la respuesta a la pregunta » ¿Mis datos siguen exactamente la distribución xy? » es siempre no.
• El objetivo aquí no puede ser determinar con certeza qué distribución sigue su muestra. El objetivo es lo que @whuber (en los comentarios) llama descripciones aproximadas parsimoniosas de los datos. Tener una distribución paramétrica específica puede ser útil como modelo de los datos.

---

Pero exploremos un poco. Usaré el excelente fitdistrplus paquete que ofrece algunas funciones interesantes para el ajuste de distribución. Usaremos la función descdist para obtener algunas ideas sobre posibles distribuciones candidatas.

library(fitdistrplus) library(logspline) x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00, 38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40, 42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40, 49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60, 45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30, 36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00, 38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34) 

Ahora usemos descdist:

descdist(x, discrete = FALSE) 

La curtosis y la asimetría al cuadrado de su muestra se trazan como un punto azul llamado » Observación «. Parece que las posibles distribuciones incluyen la distribución Weibull, Lognormal y posiblemente la distribución Gamma.

Ajustemos una Distribución de Weibull y una distribución normal:

fit.weibull <- fitdist(x, "weibull") fit.norm <- fitdist(x, "norm") 

Ahora inspeccione el ajuste para la normal:

plot(fit.norm) 

Y para el ajuste Weibull:

plot(fit.weibull) 

Ambos se ven bien, pero a juzgar por el QQ-Plot, el Weibull se ve un poco mejor, especialmente en las colas. En consecuencia, el AIC del ajuste de Weibull es menor en comparación con el ajuste normal:

fit.weibull$aic [1] 519.8537 fit.norm$aic [1] 523.3079 

---

Simulación de prueba de Kolmogorov-Smirnov

Usaré el procedimiento de @Aksakal explicado aquí para simular la estadística KS bajo el valor nulo.

n.sims <- 5e4 stats <- replicate(n.sims, { r <- rweibull(n = length(x) , shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"] ) estfit.weibull <- fitdist(r, "weibull") # added to account for the estimated parameters as.numeric(ks.test(r , "pweibull" , shape= estfit.weibull$estimate["shape"] , scale = estfit.weibull$estimate["scale"])$statistic ) }) 

El ECDF de las estadísticas KS simuladas tiene el siguiente aspecto:

plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7) grid() 

Por último, nuestro $ p $ -valor utilizando la distribución nula simulada de las estadísticas de KS es:

fit <- logspline(stats) 1 - plogspline(ks.test(x , "pweibull" , shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"])$statistic , fit ) [1] 0.4889511 

Esto confirma nuestra conclusión gráfica de que la muestra es compatible con una distribución de Weibull.

Como se explicó aquí , podemos usar bootstrapping para agregar intervalos de confianza puntuales al PDF o CDF estimado de Weibull:

xs <- seq(10, 65, len=500) true.weibull <- rweibull(1e6, shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"]) boot.pdf <- sapply(1:1000, function(i) { xi <- sample(x, size=length(x), replace=TRUE) MLE.est <- suppressWarnings(fitdist(xi, distr="weibull")) dweibull(xs, shape=MLE.est$estimate["shape"], scale = MLE.est$estimate["scale"]) } ) boot.cdf <- sapply(1:1000, function(i) { xi <- sample(x, size=length(x), replace=TRUE) MLE.est <- suppressWarnings(fitdist(xi, distr="weibull")) pweibull(xs, shape= MLE.est$estimate["shape"], scale = MLE.est$estimate["scale"]) } ) #----------------------------------------------------------------------------- # Plot PDF #----------------------------------------------------------------------------- par(bg="white", las=1, cex=1.2) plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf), xlab="x", ylab="Probability density") for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1)) # Add pointwise confidence bands quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975)) min.point <- apply(boot.pdf, 1, min, na.rm=TRUE) max.point <- apply(boot.pdf, 1, max, na.rm=TRUE) lines(xs, quants[1, ], col="red", lwd=1.5, lty=2) lines(xs, quants[3, ], col="red", lwd=1.5, lty=2) lines(xs, quants[2, ], col="darkred", lwd=2) 

#----------------------------------------------------------------------------- # Plot CDF #----------------------------------------------------------------------------- par(bg="white", las=1, cex=1.2) plot(xs, boot.cdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.cdf), xlab="x", ylab="F(x)") for(i in 2:ncol(boot.cdf)) lines(xs, boot.cdf[, i], col=rgb(.6, .6, .6, .1)) # Add pointwise confidence bands quants <- apply(boot.cdf, 1, quantile, c(0.025, 0.5, 0.975)) min.point <- apply(boot.cdf, 1, min, na.rm=TRUE) max.point <- apply(boot.cdf, 1, max, na.rm=TRUE) lines(xs, quants[1, ], col="red", lwd=1.5, lty=2) lines(xs, quants[3, ], col="red", lwd=1.5, lty=2) lines(xs, quants[2, ], col="darkred", lwd=2) #lines(xs, min.point, col="purple") #lines(xs, max.point, col="purple") 

---

Ajuste automático de distribución con GAMLSS

El gamlss paquete para R ofrece la posibilidad de probar muchas distribuciones diferentes y seleccionar mejor » de acuerdo con GAIC (el criterio de información generalizado de Akaike). La función principal es fitDist. Una opción importante en esta función es el tipo de distribuciones que se prueban. Por ejemplo, al configurar type = "realline" se probarán todas las distribuciones implementadas definidas en toda la línea real, mientras que type = "realsplus" solo probará las distribuciones definidas en la línea positiva real . Otra opción importante es el parámetro $ k $ , que es la penalización por GAIC. En el siguiente ejemplo, configuré el parámetro $ k = 2 $ , lo que significa que el » mejor se selecciona de acuerdo con el AIC clásico. Puede configurar $ k $ como desee, como $ \ log (n) $ para BIC.

library(gamlss) library(gamlss.dist) library(gamlss.add) x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00, 38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40, 42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40, 49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60, 45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30, 36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00, 38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34) fit <- fitDist(x, k = 2, type = "realplus", trace = FALSE, try.gamlss = TRUE) summary(fit) ******************************************************************* Family: c("WEI2", "Weibull type 2") Call: gamlssML(formula = y, family = DIST[i], data = sys.parent()) Fitting method: "nlminb" Coefficient(s): Estimate Std. Error t value Pr(>|t|) eta.mu -24.3468041 2.2141197 -10.9962 < 2.22e-16 *** eta.sigma 1.8661380 0.0892799 20.9021 < 2.22e-16 *** 

Según el AIC, la distribución de Weibull (más específicamente WEI2, una parametrización especial de la misma ) se ajusta mejor a los datos. La parametrización exacta de la distribución WEI2 se detalla en este documento en la página 279. Inspeccionemos el ajuste por mirando los residuos en un gráfico de gusano (básicamente un gráfico QQ sin tendencia):

Esperamos que los residuos estén cerca de la línea horizontal media y el 95% de ellos se encuentren entre la línea superior y curvas de puntos inferiores, que actúan como intervalos de confianza puntuales del 95%. En este caso, el diagrama de gusano me parece bien, lo que indica que la distribución de Weibull es un ajuste adecuado.

Comentarios

• +1 Buen análisis. Sin embargo, una pregunta. La conclusión positiva sobre la compatibilidad con una distribución principal en particular (Weibull, en este caso) permite descartar la posibilidad de una distribución de mezcla ‘ s? O necesitamos realizar un análisis de mezcla adecuado y verificar GoF para ¿descartar esa opción?
• @AleksandrBlekh Es imposible tener suficiente poder para descartar una mezcla: cuando la mezcla es de dos distribuciones casi idénticas no se puede detectar y cuando todos los componentes menos uno tienen proporciones muy pequeñas tampoco se puede detectar. Típicamente (en ausencia de una teoría que pueda sugerir una forma de distribución), uno ajusta distribuciones paramétricas para lograr descripciones aproximadas parsimoniosas de datos. Las mezclas no son nada de eso: requieren demasiados parámetros y son demasiado flexibles para el propósito.
• @whuber: +1 ¡Aprecie su excelente explicación!
• @Lourenco Miré el gráfico de Cullen y Fey. El punto azul denota nuestra muestra. Verá que el punto está cerca de las líneas de Weibull, Lognormal y Gamma (que está entre Weibull y Gamma). Después de ajustar cada una de esas distribuciones, comparé las estadísticas de bondad de ajuste usando la función gofstat y el AIC. No hay ‘ un consenso sobre cuál es la mejor manera de determinar la » mejor » distribución es. Me gustan los métodos gráficos y el AIC.
• @Lourenco ¿Te refieres al lognormal? La distribución logística (el signo » + «) está bastante alejada de los datos observados. El lognormal también sería un candidato que ‘ normalmente miraría. Para este tutorial, ‘ he elegido no mostrarlo para que la publicación sea breve. El lognormal muestra un ajuste peor en comparación con la distribución de Weibull y Normal. El AIC es 537.59 y los gráficos tampoco ‘ se ven demasiado bien.

Los gráficos son principalmente una buena forma de tener una mejor idea de cómo se ven sus datos. En su caso, recomendaría graficar la función de distribución acumulativa empírica (ecdf) contra los CDF teóricos con los parámetros que obtuvo de fitdistr ().

Hice eso una vez para mis datos y también incluí los intervalos de confianza. Aquí está la imagen que obtuve usando ggplot2 ().

La línea negra es la función de distribución acumulativa empírica y las líneas de color son CDF de diferentes distribuciones usando parámetros que obtuve usando el método de máxima verosimilitud. Uno puede ver fácilmente que la distribución exponencial y normal no se ajustan bien a los datos, porque las líneas tienen una forma diferente a la ecdf y las líneas están bastante lejos de la ecdf. Desafortunadamente, las otras distribuciones están bastante cerca. Pero yo diría que la línea logNormal es la más cercana a la línea negra. Usando una medida de distancia (por ejemplo, MSE) se podría validar la suposición.

Si solo tiene dos distribuciones en competencia (por ejemplo, seleccionando las que parecen encajar mejor en la gráfica), podría usar un Likelihood-Ratio-Test para probar qué distribuciones encajan mejor.

Comentarios

¡Bienvenido a CrossValidated! Su respuesta podría ser más útil si pudiera editarla para incluir (a) el código que utilizó para producir el gráfico y (b) cómo se leería el gráfico.

¿Qué se está trazando allí? ¿Es una especie de diagrama de exponencialidad?

Pero, ¿cómo decide qué distribución se ajusta mejor a sus datos? Solo de acuerdo con el gráfico, no pude ‘ t decirle si logNormal o weibull se ajusta mejor a sus datos.

Si desea crear un generador de números pseudoaleatorios, ¿por qué ¿No utilizas el CDF empírico? ¿Quiere dibujar números que vayan más allá de su distribución observada?

Tomando su gráfico al pie de la letra, parecería que ninguna de sus distribuciones candidatas se ajusta bien a los datos. Además, su ecdf parece tener una asíntota horizontal de menos de 0,03, lo que no ‘ no tiene sentido, por lo que ‘ no estoy seguro de que en primer lugar, realmente es un ecdf.