La simulación le pide otra cosa al ajuste

Una prueba de bondad de ajuste pregunta si un modelo describe la muestra disponible. Un modelo Monte Carlo le pide que genere valores que nunca se han visto. No son el mismo requisito, y una distribución puede satisfacer el primero y fallar el segundo.

La razón está en dónde pone la atención cada uno. Los estadísticos de bondad de ajuste están dominados por la región que contiene la mayor parte de los datos: el centro. La salida de una simulación suele estar dominada por la región opuesta: el peor percentil de un retraso, la probabilidad de superar una capacidad, el tamaño del colchón que absorbe los días malos. Al ejecutar un millón de réplicas, la respuesta que importa la deciden los pocos miles de extracciones de la cola lejana.

Esa cola es precisamente la parte de la distribución sobre la que la muestra menos sabe.

Qué tan poco dice la muestra sobre la cola

En un ajuste real de 500 observaciones (la ejecución completa), los dos modelos mejor clasificados fueron una gaussiana inversa y una lognormal. En todo el rango observado son casi idénticos:

Cuantil Empírico Gaussiana inversa Lognormal
0.50 3.3260 3.2595 3.2841
0.90 6.4504 6.6341 6.5475
0.95 8.0259 8.0596 7.9621
0.99 10.7939 11.4245 11.4915
0.999 16.4482 16.3888 17.3377

Ambos son excelentes hasta el percentil 95. En el percentil 99.9 difieren un 6%, y ahí la columna empírica carece de sentido de todos modos: con n = 500 está interpolando a partir de una única observación.

Más allá de los datos la brecha se amplía:

Umbral Empírico Gaussiana inversa Lognormal
P(X > 12) 0.0080 0.0076 0.0080
P(X > 15) 0.0040 0.0019 0.0024
P(X > 20) 0.0000 0.0002 0.0004

En 20 —por encima del máximo muestral de 16.49— un modelo afirma que ese evento es el doble de probable que el otro. Ambos pasaron todas las pruebas de bondad de ajuste. Ninguno queda contradicho por los datos, porque los datos guardan silencio ahí.

Esta es la posición honesta: con 500 observaciones hay aproximadamente 5 por encima del percentil 99, y ninguna información más allá del máximo. Cualquier afirmación que haga la simulación sobre esa región proviene del supuesto de forma de la familia ajustada, no de la muestra.

Qué hacer al respecto

Elegir la cola de forma deliberada. Puesto que los datos no pueden zanjarlo, la elección debe apoyarse en fundamentos explicables. ¿Tiene el proceso un máximo físico? ¿Existe un mecanismo que produzca valores extremos ocasionales — reintentos, escalamientos, efectos compuestos? Una lognormal, una Weibull con forma menor que 1 y una gamma implican respuestas muy distintas sobre eventos raros, y esa diferencia es un supuesto de modelado, no un resultado del ajuste.

Revisar el soporte antes que nada. Un modelo sin cota superior terminará generando un valor que el proceso no puede producir. En una ejecución larga, «terminará» llega. Si una duración no puede exceder la jornada, hay que usar una familia acotada o truncar explícitamente y documentarlo.

Ejecutar las mejores candidatas, no solo la ganadora. Es el hábito más útil de todos. Tomar cada modelo dentro de unos 2 puntos de AIC respecto al líder, correr la simulación con cada uno y comparar las salidas que realmente importan:

for name in ["inverse_gaussian", "lognormal", "inverse_gamma_3p"]:
    params = phi.sorted_distributions[name]["parameters"]
    # instanciar la distribución, extraer las réplicas
    # y registrar la métrica de la que depende la decisión

Si las respuestas coinciden, la elección entre ellas nunca fue importante y se puede dejar de lado. Si divergen, esa dispersión es el resultado. Reportar un número único a partir de una única distribución ajustada presenta una decisión de modelado como si fuera una medición.

Preferir Anderson–Darling entre las pruebas. Pondera las colas, que es la región que consume la simulación. Un modelo seleccionado solo por Kolmogorov–Smirnov se ha seleccionado por la parte de la distribución que menos importa aquí. La comparación de las tres pruebas explica por qué.

Considerar no ajustar en absoluto. Con suficientes observaciones, y si solo hace falta reproducir el rango ya observado, remuestrear la distribución empírica evita todos los problemas de extrapolación anteriores. Su limitación es exactamente la imagen especular: nunca podrá generar un valor mayor que el máximo. Conviene ajustar un modelo paramétrico cuando se necesita la cola más allá de los datos — y entonces dejar explícito que esa cola es un supuesto.

Ajustar marginales no es ajustar el proceso

Un fallo distinto, y más caro. Ajustar cada entrada por separado y muestrearlas de forma independiente descarta toda relación entre ellas. Si el tiempo de servicio y la tasa de llegadas se mueven juntos, si el retraso de hoy predice el de mañana, si dos componentes fallan por una causa común, las extracciones marginales independientes producirán una simulación demasiado bien comportada. Su varianza agregada quedará subestimada, y los escenarios extremos que justifican el modelo serán demasiado raros.

Ningún cuidado al ajustar las distribuciones individuales corrige esto. Hay que buscar correlación y autocorrelación en los datos originales antes de dar por bueno un modelo de marginales independientes.

Un protocolo de trabajo

  1. Ajustar, y leer el grupo de cabeza en lugar de la primera fila.
  2. Descartar todo lo que tenga un soporte incorrecto para el proceso.
  3. Comparar el comportamiento de las candidatas en los cuantiles que consume el modelo, no en la mediana.
  4. Ejecutar la simulación con cada candidata superviviente.
  5. Reportar la dispersión entre ellas como parte del resultado.
  6. Declarar por separado cuál es el supuesto sobre la cola y qué evidencia lo respalda, si es que hay alguna.

La distribución es una entrada del modelo, y una entrada incierta. Dar el ajuste por cerrado porque una prueba no lo rechazó desplaza esa incertidumbre fuera de la vista, no fuera de la respuesta.

Fuentes primarias