library(ggplot2)
library(readxl)
library(patchwork)

Histogramas para Distribuciones Discretas

set.seed(123) # reproducibilidad de resultados

# Preparamos los datos empíricos
valores_binomiales <- rbinom(n = 10000, size = 10, prob = 0.3)

datos_discretos <- data.frame(Datos = valores_binomiales)

# quiero concentrarlos en una tabla de frecuencias

tabla_frecuencias <- table(datos_discretos$Datos)
tabla_frecuencias <- as.data.frame(tabla_frecuencias)

tabla_frecuencias
##    Var1 Freq
## 1     0  300
## 2     1 1180
## 3     2 2326
## 4     3 2734
## 5     4 2010
## 6     5 1001
## 7     6  341
## 8     7   91
## 9     8   16
## 10    9    1
# cambiemos el nombre de las columnas
colnames(tabla_frecuencias) <- c("Valor", "Frecuencia")

tabla_frecuencias
##    Valor Frecuencia
## 1      0        300
## 2      1       1180
## 3      2       2326
## 4      3       2734
## 5      4       2010
## 6      5       1001
## 7      6        341
## 8      7         91
## 9      8         16
## 10     9          1
# factor, variable categórica
# necesitamos convertir la variable valor en numérica

tabla_frecuencias$Valor <- as.numeric(as.character(tabla_frecuencias$Valor))

tabla_frecuencias
##    Valor Frecuencia
## 1      0        300
## 2      1       1180
## 3      2       2326
## 4      3       2734
## 5      4       2010
## 6      5       1001
## 7      6        341
## 8      7         91
## 9      8         16
## 10     9          1
# Preparamos los datos teóricos
x_values <- 0:9
y_values <- dbinom(x_values, size = 10, prob = 0.3) * 10000
tabla_teorica <- data.frame (Valor = x_values, Frecuencia = y_values)

# Graficamos
ggplot() +
geom_segment (data = tabla_frecuencias, aes (x = Valor, y = Frecuencia, xend = Valor, yend = 0),
color = "steelblue", size = 3) +
geom_segment (data = tabla_teorica, aes (x = Valor, y = Frecuencia, xend = Valor, yend = 0), color =
"plum1", size = 1.5)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.


Al simular datos binomiales y construir su tabla de frecuencias, es posible visualizar la forma observada del fenómeno y contrastarla con los valores esperados según el modelo teórico.

La cercanía entre ambas representaciones sugiere un buen ajuste, lo que ilustra cómo este tipo de análisis permite validar si una distribución conocida describe adecuadamente los datos y sirve como base para modelado e inferencia estadística.

Problema de los autos (Tiempos de Interarribo)

Los tiempos de interarribo corresponden a los intervalos entre llegadas consecutivas de “clientes” a un sistema, donde un cliente puede representar personas, llamadas, paquetes o eventos, y el sistema puede ser un banco, servidor, hospital o aplicación.

Estos tiempos son fundamentales en el análisis de sistemas de espera, ya que describen la dinámica de llegada de la demanda. En el contexto de la teoría de colas, siempre que existen llegadas aleatorias y recursos limitados, surgen fenómenos de espera, por lo que entender estos intervalos permite modelar y analizar situaciones comunes en la vida real.

Este tipo de datos es clave porque define el comportamiento de la demanda: llegadas muy frecuentes pueden generar congestión, llegadas espaciadas implican capacidad ociosa, y una alta variabilidad puede producir picos de saturación. En este sentido, no basta con analizar promedios, ya que la variabilidad en los tiempos de interarribo es lo que realmente determina qué tan “estresado” está el sistema y cómo se comportan métricas como los tiempos de espera o la longitud de la cola.

Además, los tiempos de interarribo constituyen la base para modelar procesos reales y construir modelos probabilísticos, permitiendo ajustar distribuciones (como exponencial o gamma), validar supuestos como el proceso de Poisson y simular escenarios futuros. Su análisis tiene implicaciones directas en la toma de decisiones operativas, ya que permite evaluar el desempeño del sistema y responder preguntas clave como cuántos servidores son necesarios, cuándo podría saturarse el sistema o dónde se encuentran los cuellos de botella, facilitando así la optimización de recursos.

Trabajaremos con datos recopilados sobre el tiempo de arribo (\(X_i\) en minutos) entre autos que llegaron a una instalación bancaria con servicio de autos. Durante un intervalo fijo de 90 minutos, llegaron 220 autos y se registraron los tiempos de interarribo consecutivos entre los autos \(i\) y \(i+1\), para \(i=1,2,…,219\). Los datos fueron ordenados de menor a mayor y resumidos en el archivo tiempos.xlsx.

Importamos la base de datos

df_tiempos <- read_excel("tiempos.xlsx")
head(df_tiempos)
## # A tibble: 6 × 1
##   tiempos
##     <dbl>
## 1    0.01
## 2    0.01
## 3    0.01
## 4    0.01
## 5    0.01
## 6    0.01


Realizamos un análisis preliminar de la base de datos.

summary(df_tiempos) # Podenos observar que la Mediana < Media, esto implica que la distribución es asinétrica a la derecha.
##     tiempos      
##  Min.   :0.0100  
##  1st Qu.:0.1000  
##  Median :0.2700  
##  Mean   :0.3988  
##  3rd Qu.:0.5450  
##  Max.   :1.9600
hist(df_tiempos$tiempos)


Creamos tres histogramas con diferentes anchos de intervalos ∆b por ejemplo 0.05,0.075,0.1 y tres histogramas con bins aplicando las reglas. Cada histograma con su curva de densidad empírica.

g1 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               binwidth = 0.05) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)

g2 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               binwidth = 0.08) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)

g3 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               binwidth = 0.1) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)

k_sturges <- nclass.Sturges(df_tiempos$tiempos)
g4 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               bins = k_sturges) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)

k_scott <- nclass.scott(df_tiempos$tiempos)
g5 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               bins = k_scott) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)

k_FD <- nclass.FD(df_tiempos$tiempos)
g6 <- ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               bins = k_FD) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)


( g1 | g2 | g3 ) / ( g4 | g5 | g6 )


La elección del binning en un histograma no es meramente estética, ya que influye directamente en la interpretación de la distribución subyacente. Al construirlo, se está realizando una estimación de dicha distribución, por lo que el binwidth juega un papel crucial: valores demasiado grandes pueden ocultar patrones y hacer que los datos parezcan más “normales” de lo que realmente son, mientras que valores muy pequeños pueden exagerar el ruido y generar una apariencia caótica. Por ello, explorar distintas reglas y configuraciones permite evaluar la estabilidad de los patrones observados e identificar características relevantes como asimetría, multimodalidad, colas pesadas y posibles outliers, evitando así conclusiones sesgadas.

Asimismo, la incorporación de una curva de densidad no cumple solo una función visual, sino que proporciona una estimación suavizada de la distribución (por ejemplo, mediante kernel density), lo que facilita observar la forma general de los datos con menor dependencia del binning, comparar con distribuciones teóricas y fundamentar la selección de modelos probabilísticos adecuados para realizar inferencias válidas.

Concluimos que los tiempos de interarribo pudieran seguir una distribución:

  • Exponencial
  • Weibull
  • Gamma
  • Lognormal

Problema de artículos demandados

La siguiente tabla muestra valores y las frecuencias para n = 156 observaciones sobre la cantidad de artículos demandados en una semana de un inventario durante un periodo de 3 años, organizados en orden ascendente.

Cantidad 0 1 2 3 4 5 6 7 9 11
Frecuencia 59 26 24 18 12 5 4 3 3 2 
cantidades_demandadas <- c(0,1,2,3,4,5,6,7,8,9,10,11)

frecuencias <- c(59, 26, 24, 18, 12, 5, 4, 3, 0, 3, 0, 2)

df_demanda <- data.frame(Cantidad = cantidades_demandadas, Frecuencia = frecuencias)

head(df_demanda)
##   Cantidad Frecuencia
## 1        0         59
## 2        1         26
## 3        2         24
## 4        3         18
## 5        4         12
## 6        5          5
ggplot(data = df_demanda, aes(x = Cantidad, y = Frecuencia)) +
geom_segment(aes(xend = Cantidad, yend = 0), color = "skyblue", size = 3)


¿Qué distribución de probabilidad pudiera ajustar a estos datos?

  • Binomial donde el parámetro pudiera estar alrededor de 0.1
  • Geométrica donde el parámetro pudiera estar alrededor de 0.5
  • Poisson donde el paránetro pudiera estar alrededor de 0.5

Antes de pasar a la comparación de distribuciones retomemos el último tipo de estructura de datos de la lista.

Objetos de datos

En R, toda la información se almacena en objetos. Dependiendo de cómo estén organizados los datos, estos objetos pueden adoptar distintas estructuras. Existen cuatro estructuras principales en R:

I. Vectores

II. Matrices

III. Data frames

IV. Listas

Continuamos con nuestro estudio de las listas.

Listas

Una lista en R es una estructura de datos flexible que permite almacenar objetos de distinto tipo y tamaño dentro de un mismo contenedor. A diferencia de vectores o matrices, los elementos de una lista pueden ser números, texto, vectores, data frames e incluso otras listas.

Para crear una lista se utiliza la función list() cuyos argumentos son:

  1. valores: Elementos a incluir en la lista (pueden ser de diferentes tipos)
  2. nombres (opcional): Etiquetas para identificar cada elemento
# Creación de una lista
mi_lista <- list(
  nombre = "María",
  edad = 25,
  notas = c(8, 9, 10)
)

mi_lista
## $nombre
## [1] "María"
## 
## $edad
## [1] 25
## 
## $notas
## [1]  8  9 10
mi_lista[1]
## $nombre
## [1] "María"
mi_lista[[1]]
## [1] "María"
# Acceso a elementos
mi_lista$nombre      # "María"
## [1] "María"
mi_lista[[2]]        # 25
## [1] 25
mi_lista["notas"]    # sublista con notas
## $notas
## [1]  8  9 10
# Modificación de elementos
mi_lista$edad <- 26              # cambiar valor
mi_lista$notas <- c(9, 10, 10)   # modificar vector
mi_lista$ciudad <- "CDMX"        # agregar nuevo elemento

clase <- list(e1 = mi_lista,
              e2 = list(nombre = "Luis",
                        edad = 20,
                        notas = c(8, 9, 10))
              )


# Lista clase de dos estudiantes
clase
## $e1
## $e1$nombre
## [1] "María"
## 
## $e1$edad
## [1] 26
## 
## $e1$notas
## [1]  9 10 10
## 
## $e1$ciudad
## [1] "CDMX"
## 
## 
## $e2
## $e2$nombre
## [1] "Luis"
## 
## $e2$edad
## [1] 20
## 
## $e2$notas
## [1]  8  9 10
# Añadimos ciudad para estudiante 2
clase$e2$ciudad <- "Xalapa"  
clase
## $e1
## $e1$nombre
## [1] "María"
## 
## $e1$edad
## [1] 26
## 
## $e1$notas
## [1]  9 10 10
## 
## $e1$ciudad
## [1] "CDMX"
## 
## 
## $e2
## $e2$nombre
## [1] "Luis"
## 
## $e2$edad
## [1] 20
## 
## $e2$notas
## [1]  8  9 10
## 
## $e2$ciudad
## [1] "Xalapa"
# Acceso a elementos
clase$nombre                         # "Ana"
## NULL
mi_lista$estudiante2$nombre             # "María"
## NULL
mi_lista$estudiante2$notas              # notas de María
## NULL
# Modificación
mi_lista$estudiante2$edad <- 25         # cambiar edad de María
mi_lista$estudiante2$notas[1] <- 10     # modificar una nota


Las listas se abordan al final porque muchas de las funciones más importantes en R, especialmente en el contexto de modelado estadístico, devuelven sus resultados en este formato. Por ejemplo, al entrenar un modelo, el objeto resultante suele ser una lista que contiene múltiples componentes, como coeficientes, residuos, valores ajustados y métricas de desempeño. Comprender cómo acceder y manipular estos elementos es fundamental para poder interpretar correctamente los resultados.

Después de determinar una o más distribuciones de probabilidad que podría ajustarse a nuestros datos observados, ahora debemos examinar estas distribuciones detenidamente para evaluar qué tan bien representan la verdadera distribución subyacente de nuestros datos. Si varias de estas distribuciones son “representativas”, también debemos determinar cuál de ellas proporciona el mejor ajuste.

Lo que realmente intentamos hacer es determinar una distribución que sea lo suficientemente precisa para los propósitos previstos del modelo.

La función fitdist del paquete fitdistrplus en R sirve para ajustar distribuciones de probabilidad a un conjunto de datos.

Su uso principal es estimar los parámetros de una distribución teórica que mejor se ajustan a la distribución observada de tus datos. Esto lo hace utilizando diferentes métodos, siendo el más común y potente el de Máxima Verosimilitud (MLE).

library(fitdistrplus)

ajuste_exponencial <- fitdist(df_tiempos$tiempos, "exp", method = "mle")
ajuste_exponencial$estimate
##     rate 
## 2.507442
ajuste_weibull <- fitdist(df_tiempos$tiempos, "weibull", method = "mle")
ajuste_weibull$estimate
##     shape     scale 
## 1.0287171 0.4034665
ajuste_gamma <- fitdist(df_tiempos$tiempos, "gamma", method = "mle")
ajuste_gamma$estimate
##    shape     rate 
## 1.048944 2.630106
# Añadamos las líneas de densidad teóricas
k_FD <- nclass.FD(df_tiempos$tiempos)
ggplot(data = df_tiempos, aes(x = tiempos)) +
  geom_histogram(aes(y = after_stat(density)), 
               fill = "plum1", 
               color = "black", 
               bins = k_FD) +
  stat_function(fun = dexp, args = list(ajuste_exponencial$estimate), color = "blue", linewidth = 1.5) +
  stat_function(fun = dweibull, args = list(shape = ajuste_weibull$estimate[1], scale = ajuste_weibull$estimate[2]), color = "purple", linewidth = 1.5) +
  stat_function(fun = dgamma, args = list(shape = ajuste_gamma$estimate[1], rate = ajuste_gamma$estimate[2]), color = "limegreen", linewidth = 1.5) +
  geom_density (color = "deeppink", 
                linewidth = 1.5)


Criterios de selección de modelos

El AIC (Akaike Information Criterion) es una métrica utilizada para comparar modelos estadísticos que busca un equilibrio entre la bondad de ajuste de un modelo a los datos (qué tan bien predice los valores observados) y la complejidad de ese modelo (cuántos parámetros utiliza).

El modelo con el valor de AIC más bajo es el preferido. Esto se debe a que un AIC más bajo indica un mejor equilibrio entre un buen ajuste a los datos y una menor complejidad.

Si un modelo tiene muchos parámetros, podría ajustarse “demasiado bien” al ruido de los datos (sobreajuste), lo que no es deseable. El AIC penaliza ese exceso de complejidad.

ajuste_exponencial$aic # 37.36273
## [1] 37.36273
ajuste_weibull$aic # 39.68275
## [1] 39.08275
ajuste_gamma$aic # 39.04695
## [1] 39.04695

Regremos ahora el ejemplo discreto para la demanda semanal de un producto.

cantidades_expandidas <- rep(cantidades_demandadas, times = frecuencias)
cantidades_expandidas
##   [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [51]  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
##  [76]  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
## [101]  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
## [126]  3  3  4  4  4  4  4  4  4  4  4  4  4  4  5  5  5  5  5  6  6  6  6  7  7
## [151]  7  9  9  9 11 11
ajuste_geometrico <- fitdist(cantidades_expandidas, "geom", method = "mle")

ajuste_geometrico$estimate
##     prob 
## 0.345898
ajuste_poisson <- fitdist(cantidades_expandidas, "pois"
, method = "mle")
ajuste_poisson$estimate
##   lambda 
## 1.891026
ajuste_binomial <- fitdist(cantidades_expandidas, "binom", method = "mle", fix.arg = list(size = 11), start = list(prob =0.5))

ajuste_binomial$estimate
##      prob 
## 0.1719129
# La binomial tiene dos paráemtros, uno de tamaño y una proba de cada ensayo Bernouillo
# start es el punto de partida para el proceso de iteración y
# encuentre el verdadero parámetro para el éxito
ajuste_geometrico$aic
## [1] 583.673
ajuste_poisson$aic
## [1] 688.6274
ajuste_binomial$aic
## [1] 767.6454


Al utilizar herramientas como el paquete fitdistrplus, primero se ajustan distintas distribuciones candidatas de manera individual. Cada ajuste genera un objeto (generalmente una lista) con información relevante del modelo.

Posteriormente, estos objetos se pueden pasar a funciones como gofstat() para comparar su desempeño mediante distintas métricas. De esta forma, el manejo de listas se vuelve esencial para trabajar de manera eficiente con múltiples modelos y sus resultados.

# Del problema de los tiempos de interarribo
resultados_ajustes <- list()
resultados_ajustes[["Exponencial"]] <- ajuste_exponencial
resultados_ajustes[["Weibull"]] <- ajuste_weibull
resultados_ajustes[["Gamma"]] <- ajuste_gamma

# resultados_ajustes
gofstat(resultados_ajustes, fitnames = c("Exponencial", "Weibull", "Gamma"))
## Goodness-of-fit statistics
##                              Exponencial   Weibull      Gamma
## Kolmogorov-Smirnov statistic  0.04725306 0.0574567 0.05731973
## Cramer-von Mises statistic    0.08438533 0.1033474 0.10650415
## Anderson-Darling statistic    0.55790821 0.6460445 0.65947691
## 
## Goodness-of-fit criteria
##                                Exponencial  Weibull    Gamma
## Akaike's Information Criterion    37.36273 39.08275 39.04695
## Bayesian Information Criterion    40.75180 45.86089 45.82509


Los valores de Kolmogorov-Smirnov, Cramer-von Mises y Anderson-Darling indican qué tan bien se ajusta cada distribución a los datos. Valores más pequeños indican mejor ajuste.

En este caso, la distribución exponencial tiene valores menores en todas las pruebas en comparación con la Weibull, lo que sugiere que se ajusta mejor.

AIC y BIC menores indican mejor ajuste con menor complejidad.

Con base en estos resultados, la distribución exponencial se ajusta mejor a los tiempos de interarribo que la Weibull, ya que tiene valores menores en todas las métricas de ajuste.

# Adicional

ajuste_uniforme <- fitdist(df_tiempos$tiempos, "unif", method = "mle")
ajuste_uniforme$estimate
##  min  max 
## 0.01 1.96
ajuste_lognormal <- fitdist(df_tiempos$tiempos, "lnorm", method = "mle")
ajuste_lognormal$estimate
##   meanlog     sdlog 
## -1.466505  1.194720
resultados_ajustes[["Uniforme"]] <- ajuste_uniforme
resultados_ajustes[["Lognormal"]] <- ajuste_lognormal

# resultados_ajustes
gofstat(resultados_ajustes, fitnames = c("Exponencial", "Weibull", "Gamma", "Uniforme", "Lognormal"))
## Goodness-of-fit statistics
##                              Exponencial   Weibull      Gamma   Uniforme
## Kolmogorov-Smirnov statistic  0.04725306 0.0574567 0.05731973  0.4845803
## Cramer-von Mises statistic    0.08438533 0.1033474 0.10650415 24.1546564
## Anderson-Darling statistic    0.55790821 0.6460445 0.65947691        Inf
##                               Lognormal
## Kolmogorov-Smirnov statistic 0.09755645
## Cramer-von Mises statistic   0.38423217
## Anderson-Darling statistic   2.37777308
## 
## Goodness-of-fit criteria
##                                Exponencial  Weibull    Gamma Uniforme Lognormal
## Akaike's Information Criterion    37.36273 39.08275 39.04695 296.5093  61.09104
## Bayesian Information Criterion    40.75180 45.86089 45.82509 303.2874  67.86918
# Del problema de demanda de un producto
resultados_ajustes_discretos <- list()
resultados_ajustes_discretos[["Binomial"]] <- ajuste_binomial

resultados_ajustes_discretos[["Geométrica"]]<-ajuste_geometrico 

resultados_ajustes_discretos[["Poisson"]] <- ajuste_poisson

gofstat(resultados_ajustes_discretos, fitnames = c("Binomial", "Geométrico", "Poisson"))
## Chi-squared statistic:  139.1333 4.11829 87.14346 
## Degree of freedom of the Chi-squared distribution:  4 4 4 
## Chi-squared p-value:  4.327007e-29 0.3902337 5.322291e-18 
##    the p-value may be wrong with some theoretical counts < 5  
## Chi-squared table:
##      obscounts theo Binomial theo Geométrico theo Poisson
## <= 0        59      19.58667       53.960089     23.54304
## <= 1        26      44.72865       35.295402     44.52050
## <= 2        24      46.42890       23.086793     42.09470
## <= 3        18      28.91627       15.101117     26.53405
## <= 4        12      12.00618        9.877671     12.54414
## > 4         17       4.33333       18.678929      6.76356
## 
## Goodness-of-fit criteria
##                                Binomial Geométrico  Poisson
## Akaike's Information Criterion 767.6454   583.6730 688.6274
## Bayesian Information Criterion 770.6952   586.7229 691.6772

Proyecto pt.2 incluirá…

  • Seleccionar la base de datos HIV
  • Seleccionar una variable numérica para poder ajustar a una distribución.
  • Realiza un histograma escalado a densidades de esa variable usando ggplot2
  • Agrega la curva de densidad empirica.
  • Hipotetizar al menos 2 distirbuciones de probabilidad.
  • Realiza los ajustes de estas distribuciones hipotetizadas para estinar sus parámetros por “mle”
  • Superponer las curvas de densidad teórica en el histograma.
  • Aplicar el AIC para seleccionar la mejor distribución.