APUNTE DOCENTE

MODELACIÓN DEL INPUT,

VERIFICACIÓN, CALIBRACIÓN
Y VALIDACIÓN DE MODELOS
DE SIMULACIÓN
PARTE 2

Favereau, M. (2022). Modelación del Input, verificación,

calibración y validación de modelos de simulación parte 2
Universidad Andrés Bello, Santiago, Chile.
Contenido

1 Pruebas de Bondad de Ajuste

2 Problemas Comunes en Modelamiento del Input

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 2 / 32
Contenido

1 Pruebas de Bondad de Ajuste

2 Problemas Comunes en Modelamiento del Input

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 3 / 32
Pruebas de Bondad de Ajuste

Hasta el momento, hemos adivinado una distribución de probabilidad razonable y

luego estimado sus parámetros asociados (con MLE y MOM). Ahora, realicemos
una prueba formal para ver cuán exitosos han sido nuestros esfuerzos.
Particularmente, nuestra prueba es
iid
H0 : X1 , X2 , . . . , Xn ∼ pmf/pdf f (x),

a un nivel de significancia

α ≡ P(rechazar H0 |H0 verdadera) = P(error tipo I)

(usualmente, α = 0,05 o 0,01)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 4 / 32
Pruebas de Bondad de Ajuste

Un tı́pico procedimiento de prueba de bondad de ajuste (GOF) es como sigue:

1 Dividir el dominio de f (x) en k conjuntos, digamos A1 , A2 , . . . , Ak (puntos
distintos si X es discreta o intervalos si X es continua)

2 Contar el número de observaciones que entran en cada conjunto, digamos,

Oi , i = 1, 2, . . . , k. Si pi ≡ P(X ∈ Ai ), entonces Oi ∼ Bin(n, pi )

3 Dterminar el número esperado de observaciones que habrı́an de entrar en cada

conjunto si H0 fuera verdadera, es decir, Ei = E[Oi ] = npi , i = 1, 2, . . . , k

4 Calcular un estadı́stico de prueba en base a las diferentes entre los Ei ’s y

Oi ’s. El estadı́stico de prueba GOF chi-cuadrado es
k
X (Oi − Ei )2
χ20 ≡
i=1
Ei

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 5 / 32
Pruebas de Bondad de Ajuste
5 Un valor alto de χ20 indica un mal ajuste. Ası́,
I Rechazamos H0 si χ20 > χ2α,k−1−s , donde
F s es el número de parámetros desconocidos de f (x) que tienen que ser
estimados. E.g., si X ∼ N (µ, σ 2 ), entonces s = 2
F χ2α,ν es el cuantil (1 − α) de la distribución χ2ν , i.e., P(χ2ν < χ2α,ν ) = 1 − α

I Si χ20 ≤ χ2α,k−1−s , fallamos en rechazar H0

Observaciones
Recomendación tı́pica: para que funcione la prueba GOF χ2 , elegir k, n tales
que Ei ≥ 5 y n ≥ 30
Si los grados de libertad ν = k − 1 − s son muy grandes, entonces
 r 
2 2 2
χα,ν ≈ ν 1 − + zα ,
9ν 9ν

donde zα es el cuantil normal estándar apropiado

Otras pruebas GOF: Kolmogorov-Smirnov, Anderson-Darling, Shapiro-Wilk,
etc.
Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 6 / 32
Pruebas de Bondad de Ajuste

Ejemplo ilustrativo
Probar H0 : Xi ’s son U (0, 1). Suponer que tenemos n = 1000 observaciones
divididas en k = 5 intervalos:

Intervalo [0,0; 0,2) [0,2; 0,4) [0,4; 0,6) [0,6; 0,8) [0,8; 1,0]
Ei 200 200 200 200 200
Oi 179 208 222 199 192

Pk (Oi −Ei )2
Resulta que χ20 ≡ i=1 Ei = 5,27.
Tomemos α = 0,05. No hay parámetros inciertos, ası́ que s = 0. Luego, de las
tablas χ2 , tenemos
χ2α,k−1−s = χ20,05,4 = 9,49.
Ya que χ20 < χ2α,k−1−s , fallamos en rechazar H0 , y asumiremos que las
observaciones son aproximadamente uniforme.
¿Recuerdan este ejemplo? (Unidad 3)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 7 / 32
Pruebas de Bondad de Ajuste

Ejemplo discreto
Suponemos que el número de tiros libres que se necesitan hasta que se anote un
gol sigue una distribución Geo(p). Se ha recolectado una muestra aleatoria de
n = 70 conjuntos de tiros libres intentados, y el número de tiros hasta anotar un
gol se observa para cada uno de los 70 conjuntos

Tiros hasta el gol (xi ) Frecuencia (Oi )

1 34
2 18
3 2
4 9
5 7
70

Probaremos H0 : las observaciones son Geom(p)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 8 / 32
Pruebas de Bondad de Ajuste

Ejemplo discreto (cont.)

Comenzaremos estimando p a través de MLE. La función likelihood es
n
Y n
Y Pn
L(p) = f (xi ) = (1 − p)xi −1 p = (1 − p) i=1 xi −n n
p
i=1 i=1
n
X 


ln L(p) = xi − n ln(1 − p) + n ln(p)
i=1

Pn
∂ ln L(p) − i=1 xi + n n
= + = 0.
∂p 1−p p
Resolviendo para p da el MLE
1 1 70
p̂ = = P5 ≡ = 0,476
X i=1 xi · (Oi /n)
1(34) + 2(18) + 3(2) + 4(9) + 5(7)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 9 / 32
Pruebas de Bondad de Ajuste

Ejemplo discreto (cont.)

Obtengamos el estadı́stico de prueba GOF, χ20 . Crearemos una pequeña tabla,
asumiendo que p̂ = 0,476 es correcto. Por la propiedad de invarianza de los MLEs,
el número esperado de conjuntos de tiros que tienen un cierto valor x es
Ex = nP(X = x) = n(1 − p̂)x−1 p̂. Notar que combinaremos las entradas de la
última fila (≥ 5) para que las probabilidades sumen uno

x P(X = x) Ex Ox
1 0,4762 33,33 34
2 0,2494 17,46 18
3 0,1307 9,15 2
4 0,0684 4,79 9
≥5 0,0752 5,27 7
1,00 70 70

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 10 / 32
Pruebas de Bondad de Ajuste

Ejemplo discreto (cont.)

Realmente, deberı́amos combinar las dos últimas celdas también, pues
E4 = 4,79 < 5. Hagámoslo para obtener la siguiente tabla mejorada

x P(X = x) Ex Ox
1 0,4762 33,33 34
2 0,2494 17,46 18
3 0,1307 9,15 2
≥4 0,1436 10,6 16
1,00 70 70

Por lo tanto, el estadı́stico de prueba es

4
X (Ox − Ex )2 (34 − 33,33)2
χ20 = = + · · · = 9,12
x=1
Ex 33,33

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 11 / 32
Pruebas de Bondad de Ajuste

Ejemplo discreto (cont.)

Sea k = 4 el número de celdas (con las que terminamos), y sea s = 1 el número
de parámetros que estimamos.
Suponer el nivel α = 0,05. Luego, comparamos χ20 contra
χ2α,k−1−s = χ20,05,2 = 5,99.

Ya que χ20 > χ2α,k−1−s , rechazamos H0 . Esto significa que el número de tiros
probablemente no es geométrico

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 12 / 32
Pruebas de Bondad de Ajuste

Pruebas GOF para distribuciones continuas

Para el caso continuo, denotemos los intervalos Ai ≡ (ai−1 , ai ], i = 1, 2, . . . , k.
Por conveniencia, elegimos los ai ’s para asegurar que tengan intervalos de igual
probabilidad, i.e.,
1
pi = P(X ∈ Ai ) = P(ai−1 < X ≤ ai ) = ∀i.
k
En este caso, inmediatamente tenemos Ei = n/k para todo i. Luego,
k
2
X Oi − (n/k)
χ20 = .
i=1
n/k

El problema es que los ai ’s pueden depender de los parámetros desconocidos

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 13 / 32
Pruebas de Bondad de Ajuste

Ejemplo continuo
Suponer que estamos interesados en ajustar una distribución a una serie de
tiempos entre llegadas. ¿Podrán ser exponenciales?
iid
H0 : X1 , X2 , . . . , Xn ∼ Exp(λ).

Hagamos una prueba GOF χ2 con intervalos equiprobables. Esto equivale a elegir
ai ’s tales que
i
F (ai ) = P(X ≤ ai ) = 1 − e−λai = , i = 1, 2, . . . , k.
k
Esto es,  
1 i
ai = − ln 1 − , i = 1, 2, . . . , k.
λ k
Bien, pero λ es desconocido (ası́ que, estimaremos s = 1 parámetros)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 14 / 32
Pruebas de Bondad de Ajuste

Ejemplo continuo (cont.)

Buenas noticias: sabemos que el MLE es λ̂ = 1/X. Por lo tanto, por la propiedad
de invarianza, los MLEs de los ai ’s son
   
1 i i
âi = − ln 1 − = −X ln 1 − , i = 1, 2, . . . , k.
λ̂ k k

Suponer que tomamos n = 100 observaciones y las dividimos en k = 5 intervalos.

Además, suponer que la media muestral basada en las 100 observaciones es
X = 9. Luego,
âi = −9 ln(1 − 0,2i), i = 1, 2, . . . , 5

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 15 / 32
Pruebas de Bondad de Ajuste

Ejemplo continuo (cont.)

Suponer que determinamos a qué intervalo pertenece cada una de las 100
observaciones y las sumamos para obtener los Oi ’s

Intervalo (âi−1 , âi ] Oi Ei = n/k

(0; 2,01] 25 20
(2,01; 4,60] 27 20
(4,60; 8,25] 23 20
(8, 25; 14,48] 13 20
(14,48; ∞] 12 20
100 100

Pk (Oi −Ei )2
Entonces, tenemos χ20 = i=1 Ei = 9,8 y χ2α,k−1−s = χ20,05,3 = 7,81.
Ası́ que, rechazamos H0 . Estas observaciones no son exponenciales

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 16 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Hay muchas pruebas GOF que se pueden usar en vez de la χ2 . La

Kolmogorov-Smirnov (K-S) es una que trabaja bien en situaciones de pocos
datos
Probaremos
iid
H0 : X1 , X2 , . . . , Xn ∼ alguna distribución con cdf F (x)

Recordar la definición de la cdf empı́rica de los datos X1 , X2 , . . . , Xn

(Capı́tulo 1),
número de Xi ’s ≤ x
F̂n (x) ≡
n

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 17 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov
Notar que F̂n (x) es una función escalonada con saltos de altura 1/n (cada
vez que ocurre una observación)
Por ejemplo, aquı́ está la cdf empı́rica de 10, 30, 50 y 500 observaciones
N (0, 1) que generé (junto con la cdf N (0, 1) teórica)
n = 10 observaciones n = 30 observaciones
1.0 cdf (0, 1) cdf (0, 1)
cdf empírica cdf empírica
0.8
0.6
F(x)

0.4
0.2
0.0
n = 50 observaciones n = 500 observaciones
1.0 cdf (0, 1) cdf (0, 1)
cdf empírica cdf empírica
0.8
0.6
F(x)

0.4
0.2
0.0
4 2 0 2 4 4 2 0 2 4
x x

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 18 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

El lema Glivenko-Cantelli dice que F̂n (x) → F (x) para todo x cuando
n → ∞. Ası́ que, si H0 es verdadera, entonces F̂n (x) debe ser una buena
aproximación a la verdadera cdf, F (x), para n grande.
Principal pregunta: ¿La distribución empı́rica realmente respalda el supuesto
de que H0 es verdadera?
La prueba K-S rechaza H0 si

D ≡ máx |F (x) − F̂n (x)| > Dα,n ,

x∈R

donde α es el nivel de significancia, y Dα,n es un cuantil K-S que depende de

la cdf hipotética F (x)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 19 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Ejemplo pequeño
Probemos
iid
H0 : X1 , X2 , . . . , Xn ∼ U (0, 1).
Bajo el supuesto U (0, 1) para F (x), el estadı́stico K-S se simplifica:

D ≡ máx |F (x) − F̂n (x)| = máx |x − F̂n (x)|.

x∈R 0≤x≤1

Puede ser fácil ver que el máximo sólo puede ocurrir cuando x es igual a una de
las observaciones X1 , X2 , . . . , Xn , i.e., en uno de los puntos de salto de F̂n (x).
De hecho, en el i-ésimo punto de salto, F̂n (x) aumenta de i−1 n a n
i

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 20 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Ejemplo pequeño (cont.)

Antes de dar un algoritmo sencillo para calcular D, primero definamos los puntos
ordenados, X(1) ≤ X(2) ≤ · · · ≤ X(n) . Por ejemplo, si X1 = 4, X2 = 1 y X3 = 6,
entonces X(1) = 1, X(2) = 4 y X(3) = 6.
Luego, calculamos
   
i i−1
D+ ≡ máx − X(i) , D− ≡ máx X(i) − ,
1≤i≤n n 1≤i≤n n

y finalmente resulta que D = máx(D+ , D− )

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 21 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Ejemplo numérico

Xi 0,039 0,706 0,016 0,198 0,793

X(i) 0,016 0,039 0,198 0,706 0,793
i
n 0,2 0,4 0,6 0,8 1,0
i−1
n 0 0,2 0,4 0,6 0,8
i
n − X(i) 0,184 0,361 0,402 0,094 0,207
X(i) − i−1
n 0,016 – – 0,106 –

Por lo tanto, D + = 0,402, D − = 0,106,y luego D = 0,402. Si

vamos a una tabla K-S para la distribución uniforme, tenemos
Dα;n = D0,05;5 = 0,563. Ası́ que, fallamos en rechazar la
uniformidad

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 22 / 32
Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Observaciones
K-S es conservador en el sentido de que es difı́cil rechazar H0

Se puede aplicar fácilmente K-S a otras distribuciones

Muchos otras pruebas GOF: Anderson-Darling, Cramér-von Mises,

Shapiro-Wilk (S-W), etc.

S-W es especialmente apropiada para probar normalidad. También se pueden

usar técnicas gráficas como los gráficos Q-Q para evaluar normalidad

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 23 / 32
Contenido

1 Pruebas de Bondad de Ajuste

2 Problemas Comunes en Modelamiento del Input

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 24 / 32
Problemas Comunes

Siempre existe algún problema del que nadie quiere hablar, pero cada proyecto de
simulación los tiene. Uno pensarı́a que después de toda la teorı́a que hemos
analizado, siempre podrı́amos encontrar buenas distribuciones que se ajusten a
nuestros datos. Pero, ¡no exactamente!
Aquı́ hay algunos casos con los cuales hay que tener cuidado:
1 Ausencia o insuficiencia de datos
2 Datos que no se parecen a ninguna de las distribuciones usuales
3 Datos no estacionarios (distribución que cambia dinámicamente)
4 Datos multivariados y correlacionados

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 25 / 32
Ausencia o Insuficiencia de Datos

Este problema aparece más a menudo de lo que uno esperarı́a, ¡incluso en la era
del big data! Literalmente, podrı́a no haber datos disponibles, o los datos que se
tienen son horribles. Entonces, ¿qué hacemos?
Para ser honesto, no hay grandes opciones, pero aquı́ hay algunas de ellas:
1 Entrevistar a los denominados “expertos”
I Intentar de obtener al menos valores mı́nimos, máximos y “más probables”
para luego probar distribuciones como la uniforme, triangular o Pert
I Obtener cuantiles
I Al menos discutir la naturaleza de las observaciones

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 26 / 32
Ausencia o Insuficiencia de Datos

2 Si se tiene alguna idea acerca de la naturaleza de las variables aleatorias, tal

vez se pueda hacer una buena estimación de la distribución:
I ¿Discreta o continua?
I ¿Las observaciones son éxitos/fallos?–Pensar en Bernoulli, binomial,
geométrica, binomial negativa
I ¿Las observaciones adhieren a los supuestos de Poisson? Entonces, Poisson (si
se están contando llegadas) o exponencial (tiempos entre llegadas)
I ¿Las observaciones son promedios o sumas? Entonces, quizás normal
I ¿Las observaciones están acotadas? Entonces, pensar en beta
I ¿Confiabilidad o tiempos de trabajo? Quizás gamma, Weibull, lognormal, etc.
I ¿Se puede pensar en algo más a partir de las caracterı́sticas fı́sicas subyacentes
a la variable aleatoria?

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 27 / 32
Mezcla de Distribuciones
Aquı́ hay una mezcla de gaussianas forzada (color cyan). La mayorı́a de packages
no pueden ajustarlas correctamente, y hacen ajustes como el de la gaussiana
magenta

0.20 Datos reales

Gaussiana ajustada
0.15
Densidad

0.10

0.05

0.00
10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0
x

Se puede intentar modelar como una mezcla de distribuciones razonables

Más fácil: se pueden tomar muestras de la distribución empı́rica o de una
versión suavizada de la empı́rica. Esta es una forma de bootstrapping
Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 28 / 32
Datos No Estacionarios

Las tasas de llegada cambian en el tiempo–pensar en restaurantes, bancos, tráfico

en carreteras, call centers, demandas estacionales, etc. Esta variabilidad debe ser
tomada en cuenta. ¡Recordar GIGO!
Sugerencia: proceso de Poisson no homogéneo (Unidad 4 - Parte 1). Simio
puede fácilmente lidiar con esto
Se necesita modelar la tasa correctamente
En Simio se puede especificar una tasa de llegada constante distinta para
cada periodo por separado (Rate Tables)

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 29 / 32
Datos Multivariados/Correlacionados
Hasta ahora, hemos hablado de la importancia del supuesto iid. Pero, ¿todos los
fenómenos tienen que ser iid? ¿Qué pasa si los datos poseen correlación serial en
el tiempo? O, ¿si los datos son multivariados?
Un ejemplo clásico de datos altamente correlacionados espacio-temporalmente son
los caudales afluentes de una red hidroeléctrica como muestra la figura

Embalse Colbún
80 Bocatoma laguna Maule
Caudal afluente (m3/s)

60

40

20

0
abr-1989 jun-1993 ago-1997 oct-2001 dic-2005 feb-2009 abr-2014 jun-2018

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 30 / 32
Datos Multivariados/Correlacionados
Algunos otros ejemplos:
Variable aleatoria multivariada: la edad y peso de una persona están
correlacionadas
Ejemplos de correlación serial:
I Tasas de desempleo mensuales
I Entradas a un sitio web pueden estar correlacionadas si hay algo interesante ahı́
I Una pieza muy dañada puede requerir más servicio de lo normal
I Si un servidor se cansa, sus tiempos de servicio pueden ser más largos
Entonces, ¿qué necesitamos hacer?
Identificar correlación serial o si los datos son multivariados
Proponer modelos apropiados. Algunos ejemplos:
I Normal multivariada
I Modelos de mezcla de gaussianas (GMM)
I Modelos de series de tiempo, e.g., AR(p), VAR(p), ARMA(p, q), etc.
Estimar parámetros relevantes (más fácil decir que hacer). Ejemplos:
I Normal multivariada: medias y varianzas marginales, más covarianzas
I Series de tiempo: para modelos AR(p) estimar coeficientes Φ es fácil. Para
modelos más complicados deben estimarse utilizando algún software
Validar para ver si el modelo estimado es realmente bueno
Alternativa: muestrear desde distribución empı́rica (si hay datos suficientes)
Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 31 / 32
 Universidad Andrés Bello
Facultad de Ingenierı́a
Ingenierı́a Civil Industrial

Unidad 4: Modelación del Input, Verificación,

Calibración y Validación - Parte 2
ACIN213 - Simulación de Sistemas Industriales

Universidad Andrés Bello (UNAB) Unidad 4 - Parte 2 ACIN213 - Simulación de Sist. Ind. 32 / 32