Introducción

Cuando se dispone de una variable aleatoria continua (tiempos, pesos, rendimientos, etc.), un paso habitual es identificar qué modelo de probabilidad describe mejor su comportamiento. Esto permite calcular probabilidades, cuantiles y simular el proceso.

El procedimiento sigue cuatro etapas:

Etapa Pregunta Herramienta
1 ¿Cómo se ven los datos? Estadística descriptiva, gráfico de Cullen y Frey
2 ¿Qué modelos son candidatos? fitdist() + gráficos comparativos
3 ¿Cuál ajusta mejor? Criterios de información (AIC, BIC)
4 ¿El ajuste es aceptable? Pruebas de bondad de ajuste (K-S, normalidad)

Distribuciones disponibles en fitdist() (se usa el nombre corto de R):

  • Continuas: weibull, gamma, lnorm, norm, logis, cauchy, exp, beta, unif, chisq, f, t
  • Discretas: binom, geom, hyper, nbinom, pois

El documento desarrolla el proceso completo con dos ejemplos simulados (uno con datos normales y otro con datos Weibull), para luego aplicarlo a datos propios.

Librerías

library(Analitica)     # descripYG(): estadística descriptiva
library(tidyverse)     # tratamiento y visualizacion de datos
library(readxl)        # lectura de archivos Excel
library(psych)         # estadística descriptiva
library(EnvStats)      # distribuciones adicionales
library(fitdistrplus)  # ajuste de distribuciones
library(nortest)       # pruebas de normalidad

Ejemplo 1: datos con distribución Normal

Los datos se simulan desde una \(N(\mu = 220,\ \sigma = 25)\), por lo que sabemos que el modelo correcto es la Normal. El objetivo es verificar que el proceso lo detecta.

set.seed(2020)
ej1 <- data.frame(id = 1:200,
                  y  = round(rnorm(200, mean = 220, sd = 25), 1))

Etapa 1: análisis descriptivo

descripYG(ej1, y)

##     n     Mean Median       SD Kurtosis    Skewness        CV   Min Max   P25
## 1 200 219.9185  221.7 28.23187 3.065778 -0.03745778 0.1283743 143.6 300 201.7
##       P75    IQR Fence_Low Fence_High
## 1 238.175 36.475  146.9875   292.8875

El gráfico de Cullen y Frey ubica los datos según su asimetría y curtosis. Para datos normales, la observación debe quedar cerca del punto marcado como normal (asimetría 0, curtosis 3).

descdist(ej1$y, boot = 500)

## summary statistics
## ------
## min:  143.6   max:  300 
## median:  221.7 
## mean:  219.9185 
## estimated sd:  28.23187 
## estimated skewness:  -0.03774144 
## estimated kurtosis:  3.098063

Etapa 2: ajuste de modelos candidatos

fitdist() estima los parámetros de cada distribución por máxima verosimilitud.

Weibull, Gamma y Log-normal solo admiten valores estrictamente positivos (\(y > 0\)).

f1 <- fitdist(ej1$y, "weibull")
f2 <- fitdist(ej1$y, "gamma")
f3 <- fitdist(ej1$y, "norm")
f4 <- fitdist(ej1$y, "lnorm")

leyenda <- c("Weibull", "Gamma", "Normal", "Lognormal")

# Parámetros estimados
lapply(list(Weibull = f1, Gamma = f2, Normal = f3, Lognormal = f4),
       function(f) round(f$estimate, 4))
## $Weibull
##    shape    scale 
##   8.4172 232.2612 
## 
## $Gamma
##   shape    rate 
## 59.3693  0.2700 
## 
## $Normal
##     mean       sd 
## 219.9185  28.1612 
## 
## $Lognormal
## meanlog   sdlog 
##  5.3848  0.1312

La Normal estimada (\(\hat\mu\) = 219.92, \(\hat\sigma\) = 28.16) se acerca a los valores reales usados en la simulación (\(\mu = 220\), \(\sigma = 25\)).

Comparación gráfica

Gráfico Qué compara Ajuste bueno si…
denscomp Histograma vs. densidad teórica la curva sigue la forma del histograma
qqcomp Cuantiles empíricos vs. teóricos los puntos caen sobre la diagonal (sensible a las colas)
cdfcomp Distribución acumulada la curva sigue la escalera empírica
ppcomp Probabilidades empíricas vs. teóricas los puntos caen sobre la diagonal (sensible al centro) 
par(mfrow = c(2, 2))
denscomp(list(f1, f2, f3, f4), legendtext = leyenda)
qqcomp(list(f1, f2, f3, f4),   legendtext = leyenda)
cdfcomp(list(f1, f2, f3, f4),  legendtext = leyenda)
ppcomp(list(f1, f2, f3, f4),   legendtext = leyenda)

par(mfrow = c(1, 1))

Etapa 3: selección por criterios de información

Se elige el modelo con el menor AIC (criterio de Akaike) o BIC. Ambos premian el buen ajuste y penalizan la cantidad de parámetros.

gof <- gofstat(list(f1, f2, f3, f4), fitnames = leyenda)
gof
## Goodness-of-fit statistics
##                                 Weibull      Gamma     Normal  Lognormal
## Kolmogorov-Smirnov statistic 0.06716626 0.05807421 0.04170254 0.06666453
## Cramer-von Mises statistic   0.16957378 0.11258201 0.05468902 0.16087707
## Anderson-Darling statistic   1.19715845 0.60107048 0.31677079 0.86555390
## 
## Goodness-of-fit criteria
##                                 Weibull    Gamma   Normal Lognormal
## Akaike's Information Criterion 1918.106 1909.826 1906.754  1913.145
## Bayesian Information Criterion 1924.702 1916.422 1913.350  1919.741

Regla práctica: si dos modelos difieren en menos de 2 unidades de AIC (\(\Delta AIC < 2\)), su ajuste es prácticamente equivalente. En ese caso se prefiere el más simple o el que tenga sentido según el contexto de la variable.

Etapa 4: pruebas de bondad de ajuste

Kolmogorov–Smirnov

\[ H_0: \text{los datos siguen la distribución propuesta} \qquad H_1: \text{no la siguen} \] Nota: Consideración importante, el test de Kolmogorov_smirnov presenta desventajas que son convenientes tener en cuenta a la hora de la definición. En este caso pierde potencia frente a tamaños de muestra grandes (n>100)

Se rechaza \(H_0\) si el valor-p \(< 0{,}05\).

ks1 <- ks.test(ej1$y, "pweibull", shape = f1$estimate["shape"], scale = f1$estimate["scale"])
ks2 <- ks.test(ej1$y, "pgamma",   shape = f2$estimate["shape"], rate  = f2$estimate["rate"])
ks3 <- ks.test(ej1$y, "pnorm",    mean  = f3$estimate["mean"],  sd    = f3$estimate["sd"])
ks4 <- ks.test(ej1$y, "plnorm",   meanlog = f4$estimate["meanlog"], sdlog = f4$estimate["sdlog"])

ks_ej1 <- data.frame(modelo  = leyenda,
                     D       = sapply(list(ks1, ks2, ks3, ks4), function(k) unname(k$statistic)),
                     valor_p = sapply(list(ks1, ks2, ks3, ks4), function(k) k$p.value))
ks_ej1$decision <- ifelse(ks_ej1$valor_p < 0.05, "Se rechaza H0", "No se rechaza H0")
knitr::kable(ks_ej1, digits = 4)
modelo D valor_p decision
Weibull 0.0672 0.3277 No se rechaza H0
Gamma 0.0581 0.5099 No se rechaza H0
Normal 0.0417 0.8775 No se rechaza H0
Lognormal 0.0667 0.3364 No se rechaza H0

️ Como los parámetros se estimaron con los mismos datos, el valor-p de K-S es conservador (tiende a no rechazar). Para la normal, la versión corregida es la prueba de Lilliefors.

Pruebas de normalidad

Todas contrastan \(H_0\): los datos provienen de una distribución normal.

Prueba Función Observación
Shapiro–Wilk shapiro.test() La más potente hasta 25 datos; máximo 5.000 datos
Anderson–Darling ad.test() Da más peso a las colas
Cramér–von Mises cvm.test() Basada en la distribución acumulada
Lilliefors lillie.test() K-S corregido por estimar \(\mu\) y \(\sigma\) n > 70
Chi-cuadrado de Pearson pearson.test() Agrupa en clases; menos potente
Shapiro–Francia sf.test() Variante de Shapiro–Wilk; máximo 5.000 datos 
shapiro.test(ej1$y)
## 
##  Shapiro-Wilk normality test
## 
## data:  ej1$y
## W = 0.99546, p-value = 0.8147
ad.test(ej1$y)
## 
##  Anderson-Darling normality test
## 
## data:  ej1$y
## A = 0.31988, p-value = 0.5305
cvm.test(ej1$y)
## 
##  Cramer-von Mises normality test
## 
## data:  ej1$y
## W = 0.055678, p-value = 0.4296
lillie.test(ej1$y)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  ej1$y
## D = 0.042012, p-value = 0.5274
pearson.test(ej1$y)
## 
##  Pearson chi-square normality test
## 
## data:  ej1$y
## P = 10.12, p-value = 0.7534
sf.test(ej1$y)
## 
##  Shapiro-Francia normality test
## 
## data:  ej1$y
## W = 0.99573, p-value = 0.7738

Conclusión del ejemplo 1

  • K-S no rechaza \(H_0\) para: Weibull, Gamma, Normal, Lognormal.
  • Las pruebas de normalidad no rechazan \(H_0\) (valor-p de Shapiro–Wilk = 0.8147).

Lección: con datos simétricos, la Gamma y la Log-normal quedan relativamente cerca de la Normal en el AIC y también “pasan” la prueba K-S, porque cuando su forma es de campana se parecen mucho a ella. Por eso la decisión no se toma con una sola prueba: se combina el AIC, los gráficos Q-Q/P-P y las pruebas de normalidad. Aquí la evidencia apunta a la Normal, que es el modelo con el que se generaron los datos y el más simple de interpretar.

Ejemplo 2: datos con distribución Weibull

Los datos se simulan desde una \(\text{Weibull}(forma = 1{,}5;\ escala = 10)\): una distribución asimétrica a la derecha, típica de tiempos de falla y reparación.

set.seed(2026)
ej2 <- data.frame(id = 1:200,
                  y  = round(rweibull(200, shape = 1.5, scale = 10), 2))
head(ej2)
##   id     y
## 1  1  5.05
## 2  2  7.00
## 3  3 15.69
## 4  4 11.62
## 5  5  7.02
## 6  6 23.85

Funciones auxiliares

Para no repetir el código del ejemplo 1, el proceso se agrupa en funciones que reciben el vector de datos y. Se pueden reutilizar con cualquier conjunto de datos.

# Ajusta los cuatro modelos candidatos
ajustar_modelos <- function(y) {
  m <- mean(y); v <- var(y)
  list(Weibull   = fitdist(y, "weibull"),
       Gamma     = fitdist(y, "gamma", start = list(shape = m^2 / v, rate = m / v)),
       Normal    = fitdist(y, "norm"),
       Lognormal = fitdist(y, "lnorm"))
}

# Cuatro gráficos comparativos en un panel
graficar_ajustes <- function(modelos) {
  par(mfrow = c(2, 2))
  denscomp(modelos, legendtext = names(modelos))
  qqcomp(modelos,   legendtext = names(modelos))
  cdfcomp(modelos,  legendtext = names(modelos))
  ppcomp(modelos,   legendtext = names(modelos))
  par(mfrow = c(1, 1))
}

# Tabla AIC/BIC ordenada de menor a mayor AIC
tabla_aic <- function(modelos) {
  g <- gofstat(modelos, fitnames = names(modelos))
  r <- data.frame(modelo = names(modelos), AIC = round(g$aic, 2), BIC = round(g$bic, 2))
  r <- r[order(r$AIC), ]
  r$delta_AIC <- r$AIC - min(r$AIC)
  r
}

# Prueba K-S para cada modelo, con sus parámetros estimados
tabla_ks <- function(y, modelos) {
  r <- t(sapply(modelos, function(f) {
    k <- do.call(ks.test, c(list(y, paste0("p", f$distname)), as.list(f$estimate)))
    c(D = unname(k$statistic), valor_p = k$p.value)
  }))
  r <- data.frame(modelo = names(modelos), r)
  r$decision <- ifelse(r$valor_p < 0.05, "Se rechaza H0", "No se rechaza H0")
  r
}

# Seis pruebas de normalidad en una tabla
tabla_normalidad <- function(y) {
  p <- list("Shapiro-Wilk"     = shapiro.test(y),
            "Anderson-Darling" = ad.test(y),
            "Cramer-von Mises" = cvm.test(y),
            "Lilliefors"       = lillie.test(y),
            "Pearson"          = pearson.test(y),
            "Shapiro-Francia"  = sf.test(y))
  r <- data.frame(prueba      = names(p),
                  estadistico = sapply(p, function(t) unname(t$statistic)),
                  valor_p     = sapply(p, function(t) t$p.value))
  r$decision <- ifelse(r$valor_p < 0.05, "No normal", "Normal")
  r
}

Etapa 1: análisis descriptivo

descripYG(ej2, y)

##     n    Mean Median       SD Kurtosis  Skewness        CV Min   Max   P25
## 1 200 10.0845   8.93 6.637357 3.999857 0.9893702 0.6581742 0.4 35.64 5.065
##      P75  IQR Fence_Low Fence_High
## 1 14.145 9.08    -8.555     27.765

A diferencia del ejemplo 1, el histograma muestra una cola larga a la derecha y el diagrama de caja presenta valores alejados en el extremo superior.

descdist(ej2$y, boot = 500)

## summary statistics
## ------
## min:  0.4   max:  35.64 
## median:  8.93 
## mean:  10.0845 
## estimated sd:  6.637357 
## estimated skewness:  0.9968623 
## estimated kurtosis:  4.055922

En el gráfico de Cullen y Frey la observación se aleja del punto normal y queda en la zona de las distribuciones asimétricas (Gamma, Log-normal, Weibull).

Etapa 2: ajuste de modelos candidatos

modelos2 <- ajustar_modelos(ej2$y)
lapply(modelos2, function(f) round(f$estimate, 4))
## $Weibull
##   shape   scale 
##  1.5615 11.2278 
## 
## $Gamma
##  shape   rate 
## 2.0697 0.2052 
## 
## $Normal
##    mean      sd 
## 10.0845  6.6207 
## 
## $Lognormal
## meanlog   sdlog 
##  2.0504  0.8038

La Weibull estimada (forma = 1.56, escala = 11.23) se acerca a los valores reales de la simulación (forma = 1,5; escala = 10).

graficar_ajustes(modelos2)

En el gráfico Q-Q, la Normal se separa claramente de la diagonal en ambas colas.

Etapa 3: selección por criterios de información

ranking2 <- tabla_aic(modelos2)
knitr::kable(ranking2, row.names = FALSE)
modelo AIC BIC delta_AIC
Weibull 1275.81 1282.41 0.00
Gamma 1277.99 1284.59 2.18
Lognormal 1304.35 1310.95 28.54
Normal 1327.66 1334.26 51.85

Modelo con menor AIC: Weibull

Etapa 4: pruebas de bondad de ajuste

ks_ej2 <- tabla_ks(ej2$y, modelos2)
knitr::kable(ks_ej2, digits = 4, row.names = FALSE)
modelo D valor_p decision
Weibull 0.0235 0.9999 No se rechaza H0
Gamma 0.0346 0.9706 No se rechaza H0
Normal 0.0915 0.0704 No se rechaza H0
Lognormal 0.0826 0.1304 No se rechaza H0 
norm_ej2 <- tabla_normalidad(ej2$y)
knitr::kable(norm_ej2, digits = 4, row.names = FALSE)
prueba estadistico valor_p decision
Shapiro-Wilk 0.9357 0.0000 No normal
Anderson-Darling 2.9128 0.0000 No normal
Cramer-von Mises 0.4587 0.0000 No normal
Lilliefors 0.0915 0.0003 No normal
Pearson 33.7500 0.0022 No normal
Shapiro-Francia 0.9362 0.0000 No normal

Acotación: K-S contradice a las pruebas de normalidad

En este ejemplo los resultados se contradicen: K-S no rechaza que los datos sean normales, mientras que las pruebas de normalidad la rechazan en 6 de 6 casos. Sabemos que los datos provienen de una Weibull, de modo que K-S se equivoca (error tipo II).

comparacion <- data.frame(
  prueba  = c("Kolmogorov-Smirnov", "Lilliefors"),
  D       = c(ks_ej2["Normal", "D"],
              unname(lillie.test(ej2$y)$statistic)),
  valor_p = c(ks_ej2["Normal", "valor_p"],
              lillie.test(ej2$y)$p.value))
comparacion$decision <- ifelse(comparacion$valor_p < 0.05, "No normal", "Normal")
knitr::kable(comparacion, digits = 4, row.names = FALSE)
prueba D valor_p decision
Kolmogorov-Smirnov 0.0915 0.0704 Normal
Lilliefors 0.0915 0.0003 No normal

Ambas pruebas usan prácticamente el mismo estadístico \(D\) (la máxima distancia entre la distribución acumulada empírica y la teórica). La diferencia está en cómo se obtiene el valor-p:

  • K-S supone que los parámetros de la Normal (\(\mu\), \(\sigma\)) son conocidos de antemano. Aquí se estimaron con la misma muestra, de modo que la curva teórica queda “hecha a la medida” de los datos y \(D\) resulta más pequeño de lo que sería con parámetros externos. Al comparar ese \(D\) con la tabla de K-S, el valor-p sale inflado y la prueba rechaza mucho menos de lo que debería: es conservadora y pierde potencia.
  • Lilliefors es la misma prueba K-S, pero con valores críticos corregidos para el caso en que \(\mu\) y \(\sigma\) se estiman desde los datos. Con el mismo \(D\) entrega un valor-p mucho menor.
  • Además, \(D\) mide la mayor diferencia en la zona central de la distribución, mientras que Shapiro–Wilk y Anderson–Darling son más sensibles a la asimetría y a las colas, que es justamente donde la Weibull se aparta de la Normal.

La siguiente simulación lo verifica: se generan 1.000 muestras verdaderamente normales, a cada una se le estiman \(\mu\) y \(\sigma\) y se calcula \(D\). La proporción de esos \(D\) que supera al observado es el valor-p correcto (bootstrap paramétrico):

set.seed(123)
n     <- length(ej2$y)
D_obs <- ks_ej2["Normal", "D"]

D_sim <- replicate(1000, {
  z <- rnorm(n, mean = mean(ej2$y), sd = sd(ej2$y))
  unname(ks.test(z, "pnorm", mean(z), sd(z))$statistic)
})

c(D_observado        = D_obs,
  valor_p_KS_tabla   = ks_ej2["Normal", "valor_p"],
  valor_p_bootstrap  = mean(D_sim >= D_obs))
##       D_observado  valor_p_KS_tabla valor_p_bootstrap 
##        0.09147628        0.07036494        0.00000000

El \(D\) observado queda en el extremo de la distribución simulada (un valor-p bootstrap de 0 significa que ninguna de las 1.000 muestras normales alcanzó ese \(D\), es decir, \(p < 0{,}001\)): con el valor-p correcto, la normalidad sí se rechaza, en concordancia con Lilliefors, Shapiro–Wilk y Anderson–Darling.

Recomendación: para evaluar normalidad, no usar ks.test() con parámetros estimados; usar lillie.test(), shapiro.test() o ad.test(). Para los demás modelos (Weibull, Gamma, Log-normal), el valor-p de K-S debe leerse con cautela y apoyarse siempre en el AIC y en los gráficos Q-Q y P-P.

Conclusión del ejemplo 2

  • El modelo con menor AIC es Weibull.
  • K-S no rechaza \(H_0\) para: Weibull, Gamma, Normal, Lognormal.
  • Las pruebas de normalidad rechazan \(H_0\) en 6 de 6 casos: los datos no son normales.

Lección: un K-S favorable no basta (ver la acotación anterior); para la Normal manda Lilliefors / Shapiro–Wilk. Cuando los datos son asimétricos, las pruebas de normalidad lo detectan con claridad, y el AIC permite elegir entre los modelos asimétricos que compiten (Weibull, Gamma, Log-normal).

Uso del modelo elegido

Con el modelo ajustado se pueden responder preguntas de probabilidad. Por ejemplo, con la Weibull estimada:

  1. ¿Probabilidad de que una reparación tome más de 15 horas?
  2. ¿Bajo qué tiempo se completa el 90 % de las reparaciones?
forma  <- modelos2$Weibull$estimate["shape"]
escala <- modelos2$Weibull$estimate["scale"]

# a) P(Y > 15)
pweibull(15, shape = forma, scale = escala, lower.tail = FALSE)
## [1] 0.2076462
# b) Percentil 90
qweibull(0.90, shape = forma, scale = escala)
## [1] 19.15413
# Intervalos de confianza bootstrap de los parámetros
summary(bootdist(modelos2$Weibull, niter = 500))
## Parametric bootstrap medians and 95% percentile CI 
##          Median      2.5%     97.5%
## shape  1.570362  1.426149  1.744788
## scale 11.218400 10.145621 12.295983

Aplicación a datos propios

Esta sección repite el proceso con un archivo propio. Los chunks se ejecutan solo si el archivo existe en la carpeta del proyecto (opción eval); basta con cambiar el nombre en archivo.

archivo <- "NOMBREDELARCCHIVO.xlsx"
existe  <- file.exists(archivo)
existe
## [1] FALSE

Lectura y preparación

Se recomienda que el data frame tenga un nombre corto (d0, d1), aunque el archivo original tenga un nombre largo. Si los nombres de las variables no son “cómodos”, se pueden cambiar; el índice [2] corresponde a la columna que se quiere renombrar.

d0 <- as.data.frame(read_excel(archivo))

# Si los datos ya están cargados en el entorno:
# d0 <- AG_P_tiempos_de_reparacion

names(d0)[2] <- "y"
head(d0)

Con más de 5.000 datos conviene trabajar con una submuestra: las pruebas de bondad de ajuste se vuelven extremadamente sensibles (rechazan por diferencias mínimas) y shapiro.test() no acepta más de 5.000 observaciones.

set.seed(123)

if (nrow(d0) > 5000) {
  # p = porcentaje de elementos que se conservan (0,20 = 20 %)
  indice <- createDataPartition(d0$y, p = 0.20)$Resample1
  d1 <- d0[indice, , drop = FALSE]
} else {
  d1 <- d0
}

# Se eliminan datos faltantes de la variable de interés
d1 <- d1[!is.na(d1$y), , drop = FALSE]

c(n_original = nrow(d0), n_analisis = nrow(d1))

Análisis

descripYG(d1, d1$y, NULL)
descdist(d1$y, boot = 500)
modelos <- ajustar_modelos(d1$y)
lapply(modelos, function(f) round(f$estimate, 4))
graficar_ajustes(modelos)
knitr::kable(tabla_aic(modelos), row.names = FALSE)
knitr::kable(tabla_ks(d1$y, modelos), digits = 4, row.names = FALSE)
knitr::kable(tabla_normalidad(d1$y), digits = 4, row.names = FALSE)

De acuerdo con lo anterior, se observa un comportamiento … (completar la interpretación)

Resumen

Señal Datos normales (ej. 1) Datos Weibull (ej. 2)
Histograma Simétrico, forma de campana Cola larga a la derecha
Cullen y Frey Cerca del punto normal Zona Gamma / Weibull / Log-normal
Q-Q de la Normal Puntos sobre la diagonal Se curva en las colas
Menor AIC
Pruebas de normalidad No rechazan \(H_0\) Rechazan \(H_0\)

¿Cómo decidir?

  1. Mirar los datos (histograma, caja, Cullen y Frey) para proponer candidatos.
  2. Ajustar los candidatos y compararlos gráficamente (sobre todo Q-Q y P-P).
  3. Elegir el de menor AIC.
  4. Confirmar que ese modelo no es rechazado por las pruebas de bondad de ajuste.