El estudio Five easy equations for patient flow through an emergency department, de Madsen y Kofoed-Enevoldsen (2011), examina la atención de pacientes en un servicio de urgencias del Hospital Frederiksberg, en Dinamarca. Los autores trabajaron con 27142 llegadas registradas entre marzo de 2009 y febrero de 2010. Entre sus resultados se encuentran dos distribuciones útiles para analizar la variabilidad de las observaciones: una distribución de Poisson para las admisiones por hora en un período de alta ocupación y un comportamiento exponencial para los tiempos de tratamiento que superan los 15 minutos.

Se realiza la simulación por separado de dos variables, utilizando los parámetros que constan en el artículo. Se estudia cómo cambia el promedio de las observaciones al aumentar el tamaño de muestra a B = 10, 50, 100, 1000 y 10000. En cada distribución se obtiene primero su función teórica de probabilidad y después se comparan las muestras aleatorias con el valor esperado.

Variable aleatoria discreta: admisiones por hora

Se define X como el número de pacientes admitidos durante una hora en un período de alta ocupación del servicio de urgencias. Como se trata de un recuento, X es una variable aleatoria discreta y solo puede tomar valores enteros no negativos.

El artículo señala que la distribución de admisiones por hora durante un intervalo de alta utilización se ajustó bien a una distribución de Poisson con media de 4.1 admisiones por hora. Este dato corresponde a las admisiones observadas cuando el servicio estaba muy ocupado. Distribución: X ~ Poisson(4.1).

Parámetro: lambda = 4.1 admisiones por hora, valor reportado en el artículo.

Espacio muestral de la variable (soporte): X = {0, 1, 2, 3, …}. Aunque la media sea 4.1, cada observación representa una cantidad entera de admisiones.

Función de masa de probabilidad (FMP):

\[ P(X=x)=\frac{e^{-4.1}\,4.1^x}{x!}, \qquad x=0,1,2,\ldots \]

Función de distribución acumulada (FDA):

\[ F(x)=P(X\leq x)=\sum_{k=0}^{x}\frac{e^{-4.1}\,4.1^k}{k!}, \qquad x=0,1,2,\ldots \]

Valor esperado: E[X] = 4.1 admisiones por hora. Varianza: Var(X) = 4.1.

En R dpois() calcula las probabilidades individuales de la FMP, ppois() calcula la FDA y rpois() genera las admisiones simuladas. La semilla definida con set.seed() permite obtener los mismos resultados al ejecutar nuevamente el código. Los tamaños de muestra se establecen con B; son cortes de una misma secuencia de 10000 observaciones.

# Parametro y tamaños de muestra
lambda_pois <- 4.1
media_teorica_pois <- lambda_pois
B <- c(10, 50, 100, 1000, 10000)

# Una misma secuencia para observar como crece la muestra
set.seed(1234)
admisiones <- rpois(10000, lambda = lambda_pois)

En la gráfica de la FMP, las mayores probabilidades se presentan alrededor de cuatro admisiones por hora y disminuyen para cantidades mucho menores o mayores. Se observan barras porque la variable solo toma valores enteros. La FDA tiene forma escalonada: en cada número entero acumula la probabilidad correspondiente a ese valor. Se representa el intervalo de 0 a 16 admisiones para facilitar la lectura del gráfico, aunque el soporte teórico no tiene un límite superior.

x_pois <- 0:16
px_pois <- dpois(x_pois, lambda = lambda_pois)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(x_pois, px_pois, type = "h", lwd = 8, col = "#2475B8",
     main = "FMP Poisson: admisiones", xlab = "Admisiones por hora", ylab = "p(x)")
points(x_pois, px_pois, pch = 16, col = "#148E86")
plot(stepfun(x_pois, c(0, ppois(x_pois, lambda = lambda_pois))),
     main = "FDA Poisson: admisiones", xlab = "Admisiones por hora",
     ylab = "F(x)", col = "#148E86", lwd = 2, pch = 16,
     xlim = c(0, 16), ylim = c(0, 1))

par(op)

Simulaciones de la variable discreta

En los cinco cortes se toman las primeras B observaciones de la misma secuencia simulada. En cada figura se presenta el histograma a la izquierda y el promedio acumulado a la derecha, calculado con cumsum() y dividido para el número de observaciones. La línea roja punteada señala E[X] = 4,1 y los puntos verdes permiten comparar las frecuencias simuladas con la FMP de Poisson. Así se observa tanto la distribución de los recuentos individuales como la fluctuación del promedio al aumentar la muestra.

Simulación con B = 10

# Muestra de 10 observaciones
n <- 10
muestra_pois <- admisiones[1:n]
promedio_acumulado <- cumsum(muestra_pois) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_pois,
     breaks = seq(-0.5, max(muestra_pois) + 0.5, by = 1),
     prob = TRUE, col = "#B5D8F3", border = "white",
     main = "Histograma: B = 10", xlab = "Admisiones por hora",
     ylab = "Frecuencia relativa",
     xlim = c(-0.5, max(16, max(muestra_pois)) + 0.5))
points(x_pois, px_pois,
       pch = 16, col = "#148E86")
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 10",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de admisiones")
abline(h = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 50

# Muestra de 50 observaciones
n <- 50
muestra_pois <- admisiones[1:n]
promedio_acumulado <- cumsum(muestra_pois) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_pois,
     breaks = seq(-0.5, max(muestra_pois) + 0.5, by = 1),
     prob = TRUE, col = "#B5D8F3", border = "white",
     main = "Histograma: B = 50", xlab = "Admisiones por hora",
     ylab = "Frecuencia relativa",
     xlim = c(-0.5, max(16, max(muestra_pois)) + 0.5))
points(x_pois, px_pois,
       pch = 16, col = "#148E86")
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 50",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de admisiones")
abline(h = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 100

# Muestra de 100 observaciones
n <- 100
muestra_pois <- admisiones[1:n]
promedio_acumulado <- cumsum(muestra_pois) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_pois,
     breaks = seq(-0.5, max(muestra_pois) + 0.5, by = 1),
     prob = TRUE, col = "#B5D8F3", border = "white",
     main = "Histograma: B = 100", xlab = "Admisiones por hora",
     ylab = "Frecuencia relativa",
     xlim = c(-0.5, max(16, max(muestra_pois)) + 0.5))
points(x_pois, px_pois,
       pch = 16, col = "#148E86")
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 100",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de admisiones")
abline(h = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 1000

# Muestra de 1000 observaciones
n <- 1000
muestra_pois <- admisiones[1:n]
promedio_acumulado <- cumsum(muestra_pois) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_pois,
     breaks = seq(-0.5, max(muestra_pois) + 0.5, by = 1),
     prob = TRUE, col = "#B5D8F3", border = "white",
     main = "Histograma: B = 1000", xlab = "Admisiones por hora",
     ylab = "Frecuencia relativa",
     xlim = c(-0.5, max(16, max(muestra_pois)) + 0.5))
points(x_pois, px_pois,
       pch = 16, col = "#148E86")
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 1000",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de admisiones")
abline(h = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 10000

# Muestra de 10000 observaciones
n <- 10000
muestra_pois <- admisiones[1:n]
promedio_acumulado <- cumsum(muestra_pois) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_pois,
     breaks = seq(-0.5, max(muestra_pois) + 0.5, by = 1),
     prob = TRUE, col = "#B5D8F3", border = "white",
     main = "Histograma: B = 10000", xlab = "Admisiones por hora",
     ylab = "Frecuencia relativa",
     xlim = c(-0.5, max(16, max(muestra_pois)) + 0.5))
points(x_pois, px_pois,
       pch = 16, col = "#148E86")
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 10000",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de admisiones")
abline(h = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)

Los resultados de los cinco tamaños muestrales se comparan con el valor esperado teórico. El error absoluto expresa la diferencia en admisiones por hora y el error relativo muestra esa diferencia en porcentaje.

promedios_pois <- c(mean(admisiones[1:10]),
                    mean(admisiones[1:50]),
                    mean(admisiones[1:100]),
                    mean(admisiones[1:1000]),
                    mean(admisiones[1:10000]))
error_pois <- abs(promedios_pois - media_teorica_pois)
error_rel_pois <- 100 * error_pois / media_teorica_pois

tabla_pois <- data.frame(B = B,
                         Promedio = round(promedios_pois, 3),
                         Error_absoluto = round(error_pois, 3),
                         Error_porcentaje = round(error_rel_pois, 2))
print(tabla_pois, row.names = FALSE)
##      B Promedio Error_absoluto Error_porcentaje
##     10    3.900          0.200             4.88
##     50    3.920          0.180             4.39
##    100    3.670          0.430            10.49
##   1000    4.171          0.071             1.73
##  10000    4.098          0.002             0.05

Los promedios simulados fueron 3.9, 3.92, 3.67, 4.17 y 4.1 admisiones por hora para B = 10, 50, 100, 1000 y 10000, respectivamente. Estas diferencias muestran que el promedio depende de los recuentos obtenidos al azar, sobre todo al comenzar la simulación. En los cortes de 1000 y 10000, los errores relativos fueron 1.73 % y 0.05 %. El acercamiento a 4,1 no tiene que ser continuo en cada corte: la incorporación de nuevas observaciones puede producir pequeños retrocesos. La cantidad atendida en una hora sigue siendo un número entero, aunque su valor esperado sea decimal.

Variable aleatoria continua: tiempo adicional de tratamiento

Se define Y como los minutos adicionales de tratamiento después de los primeros 15 minutos, considerando únicamente tratamientos que superan ese umbral. Y es una variable aleatoria continua, ya que el tiempo puede expresarse con valores reales y no se limita a cantidades enteras.

En el estudio se encontró que, para duraciones superiores a 15 minutos, las frecuencias mostraban un decaimiento exponencial con tasa de 0,0185 por minuto. Por ese motivo, el modelo se aplica al tiempo que excede los primeros 15 minutos. No corresponde a la espera en la fila ni describe por sí solo la totalidad de los tiempos de tratamiento registrados.

Distribución: Y ~ Exponencial(0.0185).

Parámetro: lambda = 0.0185 por minuto, valor reportado en el artículo.

Espacio muestral de la variable (soporte): Y pertenece al intervalo [0, +infinito), expresado en minutos adicionales.

Función de densidad de probabilidad (FDP):

\[ f(y)=0.0185\,e^{-0.0185y},\qquad y\geq0 \]

Función de distribución acumulada (FDA):

\[ F(y)=1-e^{-0.0185y},\qquad y\geq0 \]

Valor esperado: E[Y] = 1/0.0185 = 54.05 minutos adicionales. Varianza: Var(Y) = 1/(0.0185)^2 minutos cuadrados.

La media obtenida para esta variable no equivale al promedio general de tratamiento de 61 minutos citado en el artículo. En R, dexp() calcula la FDP, pexp() calcula la FDA y rexp() genera tiempos adicionales aleatorios. El gráfico de la densidad no representa la probabilidad de un valor exacto, sino la forma de la distribución; las probabilidades se obtienen para intervalos de tiempo.

# Parametro
lambda_exp <- 0.0185
media_teorica_exp <- 1 / lambda_exp
var_exp <- 1 / lambda_exp^2

set.seed(5678)
tiempos <- rexp(10000, rate = lambda_exp)

cat("Media teorica adicional:", round(media_teorica_exp, 2), "minutos\n")
## Media teorica adicional: 54.05 minutos

La FDP comienza con una densidad mayor cerca de cero y disminuye a medida que aumenta el tiempo adicional. Esto representa la mayor frecuencia de duraciones adicionales cortas frente a duraciones muy largas. La FDA parte de cero y aumenta hacia uno; indica la probabilidad acumulada de que el tiempo adicional no supere un valor determinado. Las dos curvas corresponden a la misma distribución exponencial.

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
curve(dexp(x, rate = lambda_exp), from = 0, to = 300,
      col = "#2475B8", lwd = 2, main = "FDP Exponencial",
      xlab = "Minutos adicionales", ylab = "f(y)")
curve(pexp(x, rate = lambda_exp), from = 0, to = 300,
      col = "#148E86", lwd = 2, main = "FDA Exponencial",
      xlab = "Minutos adicionales", ylab = "F(y)", ylim = c(0, 1))

par(op)

Simulaciones de la variable continua

Se utilizan los mismos cinco tamaños de muestra: B = 10, 50, 100, 1000 y 10000. Los histogramas permiten observar la concentración de tiempos adicionales pequeños y la aparición ocasional de tiempos prolongados. La curva verde corresponde a la FDP exponencial y la línea roja punteada señala E[Y] = 54.05 minutos adicionales. A la derecha se presenta el promedio acumulado de la muestra, calculado con cumsum() y dividido para el número de observaciones.

Simulación con B = 10

# Muestra de 10 observaciones
n <- 10
muestra_exp <- tiempos[1:n]
promedio_acumulado <- cumsum(muestra_exp) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_exp, breaks = 12, prob = TRUE,
     col = "#B5D8F3", border = "white",
     main = "Histograma: B = 10", xlab = "Minutos adicionales",
     ylab = "Densidad", xlim = c(0, max(100, max(muestra_exp))))
curve(dexp(x, rate = lambda_exp), add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 10",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de minutos")
abline(h = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 50

# Muestra de 50 observaciones
n <- 50
muestra_exp <- tiempos[1:n]
promedio_acumulado <- cumsum(muestra_exp) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_exp, breaks = 12, prob = TRUE,
     col = "#B5D8F3", border = "white",
     main = "Histograma: B = 50", xlab = "Minutos adicionales",
     ylab = "Densidad", xlim = c(0, max(100, max(muestra_exp))))
curve(dexp(x, rate = lambda_exp), add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 50",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de minutos")
abline(h = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 100

# Muestra de 100 observaciones
n <- 100
muestra_exp <- tiempos[1:n]
promedio_acumulado <- cumsum(muestra_exp) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_exp, breaks = 18, prob = TRUE,
     col = "#B5D8F3", border = "white",
     main = "Histograma: B = 100", xlab = "Minutos adicionales",
     ylab = "Densidad", xlim = c(0, max(100, max(muestra_exp))))
curve(dexp(x, rate = lambda_exp), add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 100",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de minutos")
abline(h = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 1000

# Muestra de 1000 observaciones
n <- 1000
muestra_exp <- tiempos[1:n]
promedio_acumulado <- cumsum(muestra_exp) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_exp, breaks = 30, prob = TRUE,
     col = "#B5D8F3", border = "white",
     main = "Histograma: B = 1000", xlab = "Minutos adicionales",
     ylab = "Densidad", xlim = c(0, max(100, max(muestra_exp))))
curve(dexp(x, rate = lambda_exp), add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 1000",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de minutos")
abline(h = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)

Simulación con B = 10000

# Muestra de 10000 observaciones
n <- 10000
muestra_exp <- tiempos[1:n]
promedio_acumulado <- cumsum(muestra_exp) / (1:n)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(muestra_exp, breaks = 30, prob = TRUE,
     col = "#B5D8F3", border = "white",
     main = "Histograma: B = 10000", xlab = "Minutos adicionales",
     ylab = "Densidad", xlim = c(0, max(100, max(muestra_exp))))
curve(dexp(x, rate = lambda_exp), add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

plot(promedio_acumulado, type = "l", log = "x", lwd = 2,
     col = "#2475B8", main = "Promedio acumulado: B = 10000",
     xlab = "Observaciones acumuladas (escala log)",
     ylab = "Promedio de minutos")
abline(h = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)

Los promedios muestrales de los minutos adicionales se comparan con la media teórica de 54.05 minutos. También se calculan los errores absolutos y porcentuales.

promedios_exp <- c(mean(tiempos[1:10]),
                   mean(tiempos[1:50]),
                   mean(tiempos[1:100]),
                   mean(tiempos[1:1000]),
                   mean(tiempos[1:10000]))
error_exp <- abs(promedios_exp - media_teorica_exp)
error_rel_exp <- 100 * error_exp / media_teorica_exp

tabla_exp <- data.frame(B = B,
                        Promedio = round(promedios_exp, 3),
                        Error_absoluto = round(error_exp, 3),
                        Error_porcentaje = round(error_rel_exp, 2))
print(tabla_exp, row.names = FALSE)
##      B Promedio Error_absoluto Error_porcentaje
##     10   64.220         10.166            18.81
##     50   59.051          4.997             9.25
##    100   65.261         11.207            20.73
##   1000   55.314          1.260             2.33
##  10000   55.168          1.114             2.06

Los promedios obtenidos para B = 10, 50, 100, 1000 y 10000 fueron 64.22, 59.05, 65.26, 55.31 y 55.17 minutos adicionales, respectivamente. Se comparan con los 54.05 minutos del valor esperado. En B = 1000 y B = 10000, los errores relativos fueron 2.33 % y 2.06 %. En las muestras pequeñas, un tiempo adicional muy prolongado puede desplazar bastante el promedio. A medida que se incorporan observaciones, la influencia de un solo caso suele ser menor y el promedio acumulado tiende a acercarse a la media teórica.

Análisis de las simulaciones

¿Por qué se simulan los datos?

La simulación permite examinar cómo funcionan las distribuciones de probabilidad y cómo cambia el promedio muestral sin disponer de la base de datos individual del hospital. En este caso se conocen dos parámetros publicados: 4.1 admisiones por hora para una situación de alta ocupación y 0.0185 por minuto para el tiempo adicional de los tratamientos que exceden 15 minutos. A partir de esas distribuciones, R genera nuevas observaciones aleatorias que hacen posible comparar muestras pequeñas y grandes bajo las mismas condiciones del modelo.

El objetivo no es sustituir las observaciones reales, sino comprender qué parte de las diferencias entre los promedios simulados puede atribuirse al azar. En particular, se busca comprobar si una muestra de 10 o 50 observaciones proporciona una media tan estable como una de 1000 o 10000, y reconocer por qué los resultados varían aun cuando el parámetro poblacional permanece fijo.

¿Qué restricciones existen al tomar muestras aleatorias de distintos tamaños?

El primer requisito es que las observaciones simuladas sean independientes y se generen con una misma ley de probabilidad. Para la distribución de Poisson, mantener lambda = 4.1 significa representar condiciones comparables de alta ocupación; no sería adecuado utilizar esa tasa como si describiera todas las horas del hospital. En la exponencial, lambda = 0.0185 describe únicamente el tiempo adicional después de 15 minutos. Estos límites deben tenerse en cuenta al interpretar los resultados.

El tamaño de muestra también importa. Con B = 10 o B = 50, pocos valores altos o bajos pueden modificar notablemente el promedio. Con B = 100 se dispone de más información, aunque todavía pueden observarse oscilaciones. Al aumentar a B = 1000 o B = 10000, los promedios suelen ser menos sensibles a una observación aislada. Sin embargo, el error no disminuye obligatoriamente en cada corte. En una investigación real, conseguir más registros requiere tiempo, acceso a la información y revisión de datos faltantes; el artículo, por ejemplo, informa que para alrededor de 4000 pacientes no constaba la hora de admisión y se aplicó un procedimiento de imputación. En R se pueden generar muestras grandes con facilidad, pero su validez depende del modelo utilizado.

¿En qué tamaño de muestra se aprecia la estabilización del promedio?

Al comparar los cinco cortes, se toma B = 1000 como referencia inicial de estabilización visual en ambas variables: en ese tramo las oscilaciones acumuladas se vuelven más pequeñas en comparación con las primeras observaciones. El corte B = 10000 permite observar una trayectoria más sostenida alrededor de E[X] = 4,1 y E[Y] = 54.05, respectivamente. En la distribución de Poisson, el error relativo al terminar B = 1000 es 1.73 %; en la exponencial es 2.33 %. Por ello, este punto se interpreta como una referencia del ejercicio y no como una regla fija para cualquier muestra o distribución.

Como apoyo para interpretar la estabilización se considera la desviación estándar teórica del promedio muestral en relación con la media. En Poisson, esa dispersión relativa es aproximadamente 4,94 % con B = 100 y 1,56 % con B = 1000; en la exponencial pasa de 10,00 % a 3,16 % en los mismos cortes. Así, B = 1000 es el primero de los tamaños seleccionados en el que las dos distribuciones tienen una dispersión teórica relativa menor al 5 %. Este 5 % es un criterio de comparación adoptado para la actividad, no un resultado ni una exigencia del estudio hospitalario. Tampoco debe confundirse con el error observado en una simulación particular, que se presenta en las tablas anteriores.

La Ley de los Grandes Números explica este comportamiento. Cuando las observaciones son independientes y siguen una misma distribución con media finita, el promedio de una muestra grande tiende a acercarse a su valor esperado. Esta propiedad no significa que la aleatoriedad desaparezca en B = 1000 ni que los promedios mejoren en todos los pasos, sino que las fluctuaciones respecto de la media teórica tienden a reducirse conforme aumenta la cantidad de observaciones.

Teorema del Límite Central

El Teorema del Límite Central se analiza mediante las medias de muchas muestras independientes. Con replicate() se calculan 1000 promedios de cada distribución para n = 10, 100, 1000 y 10000, usando mean(), rpois() y rexp(). En esta parte, n corresponde a las observaciones de cada muestra y las 1000 repeticiones permiten construir un histograma de medias. Las 1000 repeticiones son una decisión del procedimiento de simulación, no un dato del artículo.

Los cuatro histogramas de cada distribución muestran cómo cambia la forma y la dispersión de esos promedios. La curva verde indica la aproximación normal teórica, calculada con dnorm(), y la línea roja marca el valor esperado. Se utiliza una escala horizontal propia para cada panel, por lo que es necesario observar también las desviaciones estándar de las tablas para comparar cuánto se estrecha la distribución. La aproximación normal corresponde a las medias muestrales, no a que las admisiones o los tiempos de tratamiento individuales cambien de distribución.

Distribución de los promedios muestrales de Poisson.

# Teorema del límite central: Poisson
repeticiones <- 1000
set.seed(2345)

# 1000 promedios para cada tamaño de muestra
medias_pois_10 <- replicate(repeticiones, mean(rpois(10, lambda_pois)))
medias_pois_100 <- replicate(repeticiones, mean(rpois(100, lambda_pois)))
medias_pois_1000 <- replicate(repeticiones, mean(rpois(1000, lambda_pois)))
medias_pois_10000 <- replicate(repeticiones, mean(rpois(10000, lambda_pois)))

op <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

hist(medias_pois_10, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Poisson: n = 10",
     xlab = "Promedio de admisiones", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_pois, sd = sqrt(lambda_pois / 10)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

hist(medias_pois_100, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Poisson: n = 100",
     xlab = "Promedio de admisiones", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_pois, sd = sqrt(lambda_pois / 100)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

hist(medias_pois_1000, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Poisson: n = 1000",
     xlab = "Promedio de admisiones", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_pois, sd = sqrt(lambda_pois / 1000)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

hist(medias_pois_10000, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Poisson: n = 10000",
     xlab = "Promedio de admisiones", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_pois, sd = sqrt(lambda_pois / 10000)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_pois, col = "red", lty = 2, lwd = 2)

par(op)
n_tlc <- c(10, 100, 1000, 10000)
sd_teorica_pois <- sqrt(lambda_pois / n_tlc)
sd_simulada_pois <- c(sd(medias_pois_10), sd(medias_pois_100),
                     sd(medias_pois_1000), sd(medias_pois_10000))

tabla_tlc_pois <- data.frame(n = n_tlc,
                             DE_teorica = round(sd_teorica_pois, 4),
                             DE_simulada = round(sd_simulada_pois, 4))
print(tabla_tlc_pois, row.names = FALSE)
##      n DE_teorica DE_simulada
##     10     0.6403      0.6087
##    100     0.2025      0.2017
##   1000     0.0640      0.0665
##  10000     0.0202      0.0199

En Poisson, con n = 10 los promedios pueden diferir notablemente entre una repetición y otra. Al aumentar n a 100, 1000 y 10000, se agrupan más cerca de 4,1 admisiones por hora y la distribución de esos promedios muestra una forma aproximadamente normal. La desviación estándar teórica pasa de 0.64 para n = 10 a 0.02 para n = 10000. Esto describe la variabilidad entre promedios de muestras independientes, no la variabilidad de las admisiones individuales.

Distribución de los promedios muestrales de la exponencial.

# Teorema del límite central: Exponencial
set.seed(3456)

# 1000 promedios para cada tamaño de muestra
medias_exp_10 <- replicate(repeticiones, mean(rexp(10, rate = lambda_exp)))
medias_exp_100 <- replicate(repeticiones, mean(rexp(100, rate = lambda_exp)))
medias_exp_1000 <- replicate(repeticiones, mean(rexp(1000, rate = lambda_exp)))
medias_exp_10000 <- replicate(repeticiones, mean(rexp(10000, rate = lambda_exp)))

op <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

hist(medias_exp_10, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Exponencial: n = 10",
     xlab = "Promedio de minutos", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_exp, sd = sqrt(var_exp / 10)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

hist(medias_exp_100, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Exponencial: n = 100",
     xlab = "Promedio de minutos", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_exp, sd = sqrt(var_exp / 100)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

hist(medias_exp_1000, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Exponencial: n = 1000",
     xlab = "Promedio de minutos", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_exp, sd = sqrt(var_exp / 1000)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

hist(medias_exp_10000, breaks = 20, prob = TRUE,
     col = "#B5D8F3", border = "white", main = "Exponencial: n = 10000",
     xlab = "Promedio de minutos", ylab = "Densidad")
curve(dnorm(x, mean = media_teorica_exp, sd = sqrt(var_exp / 10000)),
      add = TRUE, col = "#148E86", lwd = 2)
abline(v = media_teorica_exp, col = "red", lty = 2, lwd = 2)

par(op)
sd_teorica_exp <- sqrt(var_exp / n_tlc)
sd_simulada_exp <- c(sd(medias_exp_10), sd(medias_exp_100),
                    sd(medias_exp_1000), sd(medias_exp_10000))

tabla_tlc_exp <- data.frame(n = n_tlc,
                            DE_teorica = round(sd_teorica_exp, 4),
                            DE_simulada = round(sd_simulada_exp, 4))
print(tabla_tlc_exp, row.names = FALSE)
##      n DE_teorica DE_simulada
##     10    17.0934     16.8122
##    100     5.4054      5.3966
##   1000     1.7093      1.7001
##  10000     0.5405      0.5567

En la exponencial, la diferencia de forma entre los tamaños muestrales puede ser más evidente. Con n = 10, los promedios todavía presentan cierta asimetría, debido a la posibilidad de observar tiempos adicionales largos. Con n = 100, 1000 y 10000, su distribución se aproxima progresivamente a una normal y se concentra cerca de 54.05 minutos adicionales. La desviación estándar teórica se reduce de 17.09 minutos con n = 10 a 0.54 minutos con n = 10000.

Los histogramas del Teorema del Límite Central muestran la distribución de los promedios de muestras repetidas, mientras que los histogramas anteriores representan los valores individuales simulados. La Ley de los Grandes Números permite interpretar el acercamiento del promedio acumulado al parámetro poblacional y el Teorema del Límite Central explica por qué, al utilizar muestras grandes, la distribución de los promedios se aproxima a la normal.

Conclusión

En los dos ejemplos se observa que el tamaño de muestra influye en la estabilidad del promedio, aunque las variables se midan de manera diferente. Las admisiones por hora se representan con recuentos enteros y una distribución de Poisson; el tiempo adicional de tratamiento se mide en minutos y sigue el modelo exponencial adoptado. Los histogramas permiten reconocer sus formas de distribución y las gráficas acumuladas muestran el efecto de incrementar B. La comparación de los cinco cortes sugiere una estabilización visual más clara alrededor de 1000 observaciones, confirmada con mayor detalle al ampliar la simulación a 10000. Esta conclusión se limita a las muestras generadas y a los parámetros que se emplearon, sin constituir una nueva estimación del funcionamiento real del hospital.

Referencia bibliográfica: Madsen, T. L., y Kofoed-Enevoldsen, A. (2011). Five easy equations for patient flow through an emergency department. Danish Medical Bulletin, 58(10), A4318.