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

133

Tengo un conjunto de datos y me gustaría averiguar qué distribución se ajusta mejor a mis datos.

Utilicé la fitdistr()función para estimar los parámetros necesarios para describir la distribución supuesta (es decir, Weibull, Cauchy, Normal). Usando esos parámetros, puedo realizar una prueba de Kolmogorov-Smirnov para estimar si los datos de mi muestra son de la misma distribución que mi distribución supuesta.

Si el valor p es> 0.05, puedo suponer que los datos de la muestra se extraen de la misma distribución. Pero el valor p no proporciona ninguna información sobre la bondad del ajuste, ¿no es así?

Entonces, en caso de que el valor p de mis datos de muestra sea> 0.05 para una distribución normal y una distribución weibull, ¿cómo puedo saber qué distribución se ajusta mejor a mis datos?

Esto es básicamente lo que he hecho:

> 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

Los valores p son 0.8669 para la distribución de Weibull y 0.5522 para la distribución normal. Por lo tanto, puedo suponer que mis datos siguen un Weibull y una distribución normal. Pero, ¿qué función de distribución describe mejor mis datos?


Refiriéndome a elevendollar encontré el siguiente código, pero no sé cómo interpretar los resultados:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
tobibo
fuente
55
¿Por qué le gustaría averiguar qué distribución se ajusta mejor a sus datos?
Roland
66
Porque quiero generar números pseudoaleatorios siguiendo la distribución dada.
tobibo
66
No puede usar KS para verificar si una distribución con parámetros encontrados desde el conjunto de datos coincide con el conjunto de datos. Vea el # 2 en esta página, por ejemplo, más alternativas (y otras formas en que la prueba KS puede ser engañosa).
tpg2114
Otra discusión aquí con ejemplos de código sobre cómo aplicar la prueba KS cuando los parámetros se estiman a partir de la muestra.
Aksakal
1
I used the fitdistr() function ..... ¿Cuál es la fitdistrfunción? Algo de Excel? ¿O algo que escribiste en C?
Wolfies

Respuestas:

162

Primero, aquí hay algunos comentarios rápidos:

  • Los valores de una prueba de Kolmovorov-Smirnov (prueba KS) con parámetros estimados serán bastante erróneos. Desafortunadamente, no puede ajustar una distribución y luego usar los parámetros estimados en una prueba de Kolmogorov-Smirnov para analizar su muestra.p
  • Su muestra nunca seguirá exactamente una distribución específica. Entonces, incluso si sus valores de la prueba KS fueran válidos y , 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 nop>0.05
  • 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 hagamos un poco de exploración. Usaré el excelente fitdistrpluspaquete que ofrece algunas funciones agradables para el ajuste de distribución. Utilizaremos la función descdistpara obtener algunas ideas sobre posibles distribuciones de candidatos.

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)

Descdist

La curtosis y el sesgo al cuadrado de su muestra se representa 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 lo normal:

plot(fit.norm)

Ajuste normal

Y para el ajuste Weibull:

plot(fit.weibull)

Weibull fit

Ambos se ven bien, pero a juzgar por el QQ-Plot, el Weibull puede verse un poco mejor, especialmente en las colas. En consecuencia, el AIC del ajuste Weibull es más bajo 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 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 es el siguiente:

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

Estadísticas KS simuladas

Finalmente, nuestro valor usando la distribución nula simulada de las estadísticas KS es:p

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 explica aquí , podemos usar bootstrapping para agregar intervalos de confianza puntuales al PDF o CDF Weibull estimado:

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)

CI_Density

#-----------------------------------------------------------------------------
# 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")

CI_CDF


Distribución automática con GAMLSS

El gamlsspaquete Rofrece la posibilidad de probar muchas distribuciones diferentes y seleccionar la "mejor" según el 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, la configuración type = "realline"probará todas las distribuciones implementadas definidas en toda la línea real, mientras type = "realsplus"que solo intentará las distribuciones definidas en la línea positiva real. Otra opción importante es el parámetro , que es la penalización para el GAIC. En el ejemplo a continuación, configuro el parámetro que significa que la "mejor" distribución se selecciona de acuerdo con el AIC clásico. Puede configurar para cualquier cosa que desee, comokk=2klog(n) para el 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 la 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 WEI2se detalla en este documento en la página 279. Inspeccionemos el ajuste observando los residuos en un diagrama de gusanos (básicamente un diagrama QQ de tendencia):

WormPlot

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

COOLSerdash
fuente
1
+1 Buen análisis. Una pregunta, sin embargo. ¿La conclusión positiva sobre la compatibilidad con una distribución principal particular (Weibull, en este caso) permite descartar la posibilidad de la presencia de una distribución de mezcla? ¿O debemos realizar un análisis de mezcla adecuado y verificar GoF para descartar esa opción?
Aleksandr Blekh
18
@AleksandrBlekh Es imposible tener el poder suficiente para descartar una mezcla: cuando la mezcla tiene 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), se ajustan las distribuciones paramétricas para lograr descripciones aproximadas parsimoniosas de datos. Las mezclas no son ninguna de ellas: requieren demasiados parámetros y son demasiado flexibles para el propósito.
whuber
44
@whuber: +1 ¡Aprecio tu excelente explicación!
Aleksandr Blekh
1
@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 utilizando la función gofstaty el AIC. No hay consenso sobre cuál es la mejor manera de determinar la "mejor" distribución. Me gustan los métodos gráficos y el AIC.
COOLSerdash
1
@Lourenco ¿Te refieres a 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, elegí no mostrarlo para mantener la publicación corta. El lognormal muestra un peor ajuste en comparación con la distribución Weibull y Normal. El AIC es 537.59 y los gráficos tampoco se ven muy bien.
COOLSerdash
15

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

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

ingrese la descripción de la imagen aquí

La línea negra es la función empírica de distribución acumulativa y las líneas coloreadas son archivos PDF de diferentes distribuciones que utilizan parámetros que obtuve utilizando el método de Máxima Verosimilitud. Se puede ver fácilmente que la distribución exponencial y normal no se ajusta 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. Lamentablemente, las otras distribuciones son bastante cercanas. 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 el supuesto.

Si solo tiene dos distribuciones en competencia (por ejemplo, elegir las que parecen encajar mejor en la gráfica), podría usar una Prueba de relación de probabilidad para probar qué distribuciones se ajustan mejor.

elevendollar
fuente
20
¡Bienvenido a CrossValidated! Su respuesta podría ser más útil si pudiera editarla para incluir (a) el código que usó para producir el gráfico y (b) cómo se leería el gráfico.
Stephan Kolassa
2
¿Qué se está tramando allí? ¿Es eso una especie de trama de exponencialidad?
Glen_b
1
Pero, ¿cómo decide qué distribución se ajusta mejor a sus datos? Solo de acuerdo con el gráfico, no podría decirle si logNormal o weibull se ajusta mejor a sus datos.
tobibo
44
Si desea crear un generador de números pseudoaleatorios, ¿por qué no utilizar el cdf empírico? ¿Desea dibujar números que van más allá de su distribución observada?
elevendollar
66
Tomando su gráfico al pie de la letra, parece que ninguna de sus distribuciones candidatas se ajusta bien a los datos. Además, su ecdf parece tener una asíntota horizontal a menos de 0.03, lo que no tiene sentido, por lo que no estoy seguro de que realmente sea un ecdf en primer lugar.
Hong Ooi