Simulación – Método Monte Carlo -
1
Feller: Pero… qué hacen ustedes? Santos: Simulamos, Feller. Simulamos ...
Introducción
¿Qué es Simular?
2
Algunas definiciones De acuerdo a la RAE, simular es la acción de “representar algo, fingiendo o imitando lo que no es”. Esta definición no aplica con exactitud a las técnicas de simulación que estudiamos. Los modelos de simulación representan un sistema o una realidad de manera simplificada con la finalidad de entender su comportamiento, determinar resultados potenciales, evaluar estrategias y tomar decisiones.
A diferencia del análisis de sensibilidad o la utilización de escenarios, la simulación permite explorar un número alto de escenarios y condensar sus resultados en un gráfico o probabilidad de ocurrencia de un resultado. 3
Simulación En la actualidad, casi cualquier situación puede analizarse a través
de la simulación. Desde un proyecto de inversión y su respectiva evaluación, realizar análisis de riesgo o situaciones que se pueden describir en el contexto de una fila de espera (muchísimas situaciones se pueden
considerar como una forma de fila de
espera, recordar el ejemplo del estacionamiento de la clase pasada).
En resumen, la simulación es un experimento estadístico y por lo tanto, sus resultados se deben interpretar con las pruebas
estadísticas adecuadas. 4
Simulación vs. Escenarios En el análisis de escenarios se plantean posibles “estados de la naturaleza” respecto a una variable incierta. Si las variables a analizar son varias, la creación de escenarios puede ser engorrosa.
Por ejemplo, si analizamos dos variables y planteamos cuatro resultados posibles para cada una, obtenemos 16 escenarios posibles. Si agregáramos una variable más el número trepa a 64.
La simulación de las variables críticas (definidos los rangos de valores que pueden adoptar y la distribución de probabilidad acumulada) permite analizar las combinaciones posibles de los diferentes comportamientos de la variables críticas, estimando un resultado final con un determinado nivel de confianza. 5
Simulación – Método Monte Carlo El método de simulación Monte Carlo está dirigido hacia la estimación de parámetros estocásticos o determinísticos con base en el muestreo aleatorio.
El método Monte Carlo permite obtener soluciones de problemas matemáticos o físicos por medio del análisis de distribuciones de variables aleatorias. Se realizan las pruebas aleatorias utilizando para tal fin números aleatorios. En el método Monte Carlo el elemento tiempo no es un factor pertinente.
6
Por qué Monte Carlo? El método lleva este nombre debido a que la ruleta es un generador
simple de números aleatorios. Mientras que hoy en día la capital mundial del juego posiblemente sea Las Vegas, durante mucho tiempo lo fue Mónaco y su Casino de Monte Carlo. Como herramienta de investigación, proviene del trabajo realizado en el desarrollo de la bomba atómica durante WW2 en el Laboratorio Nacional de Los Álamos en EE.UU. Este trabajo conllevaba la simulación
de
problemas
probabilísticos
de
hidrodinámica
concernientes a la difusión de neutrones en el material de fisión. Esta difusión posee un comportamiento eminentemente aleatorio. 7
Simulación Montecarlo - Ejemplo (Taha) Veamos un ejemplo simple: estimación del área de un círculo de ecuación (x – 1)2 + (y – 2)2 = 25 Con r = 5 cm y centro en (x, y) = (1, 2). El procedimiento para estimar el área consiste en encerrar el circulo en un cuadrado cuyo lado sea igual al diámetro del círculo, como en la figura.
La estimación del área del círculo se basa en la hipótesis de que todos los puntos dentro del cuadrado tienen igual probabilidad de presentarse en un muestreo aleatorio. Siendo N los puntos que componen la muestra y M los que resultan estar dentro del círculo: Estim. del área del círculo = Área del cuadrado . M / N 8
Simulación Montecarlo - Ejemplo (Taha) Se determina el punto aleatorio (x, y) de la siguiente manera (modelo): x = - 4 + 10 . R1 y = - 3 + 10 . R2
Siendo R1 y R2 una pareja de números aleatorios. Luego, se debe verificar: (x – 1)2 + (y – 2)2 < 25 para cada uno de los pares de puntos aleatorios que conforman la muestra. El tamaño de la muestra finalmente dependerá del error admisible en la determinación del parámetro y en el costo asociado. Con las herramientas computacionales actuales (softs comoSimulAR), esto pierde relevancia. 9
Tipos de simulación La ejecución actual de la simulación suele basarse en la idea de usar el muestreo conjuntamente con el método Monte Carlo. Existen dos clases de simulación: Modelos Continuos, que manejan sistemas cuyo comportamiento cambia continuamente con el tiempo. Un ejemplo característico es el estudio de la dinámica demográfica mundial. Modelos discretos, relacionados principalmente con el estudio de líneas de espera, cuyo objetivo es determinar medidas como la longitud de la cola y el tiempo de espera promedio. Esas medidas sólo cambian cuando entra o sale un cliente del sistema. En los lapsos intermedios nada sucede en el sistema, desde el punto de vista de recolección de datos estadísticos. Los instantes en los que suceden los cambios, en puntos discretos en el tiempo, dan lugar al nombre de simulación de evento discreto. 10
Ocurrencia de eventos en simulación ¿Cómo
determina la simulación el momento de ocurrencia de los eventos?
Los
eventos de llegada están separados por el intervalo entre llegadas sucesivas y los eventos de salida son función del tiempo de servicio en la instalación. Estos
tiempos pueden ser determinísticos (por ejemplo, un tren llega a la estación cada 5 minutos) o probabilísticos (por ejemplo, la llegada aleatoria de clientes a un banco, dada por 1/λ ). Si
el tiempo entre eventos es determinísta, la determinación de sus tiempos es ocurrencia directa. Si es probabilístico se usa un procedimiento especial para muestrear a partir de la distribución correspondiente de probabilidades.
11
Muestreo a partir de distribuciones de probabilidad Cuando
el intervalo t entre eventos sucesivos es estocástico, la simulación
requiere de generar muestras aleatorias sucesivas (t = t1, t2, …) a partir de una distribución de probabilidades f(t).
Se
utilizan diferentes métodos para generar estas muestras aleatorias,
cada uno de acuerdo a la complejidad del tipo de distribución de densidad de probabilidades. El más común se denomina Método de la transformada inversa, y es de utilidad para distribuciones exponencial, uniforme y con
algunas consideraciones, las de Poisson, Erlang y Normal (llamado método de convolución).
Todos
los métodos se basan en el uso de números aleatorios uniformes
(0,1) independientes e idénticamente distribuidos. 12
Método de la transformada inversa Se
trata de obtener una muestra aleatoria x partiendo de la función de
densidad de probabilidad f(x) (continua o discreta).
Con
este método se determina primero una expresión cerrada de la función
de densidad acumulada F(x), en donde 0 ≤ F(x) ≤ 1 para todos los valores definidos de y.
El
número aleatorio R es una variable obtenida de una distribución
uniforme (0,1). Suponiendo que F-1 es la inversa de F, el método de la transformada inversa consta de dos pasos:
13
Generar el número aleatorio (0,1), R.
Calcular la muestra x = F-1(R).
Método de la transformada inversa Como
se puede observar en el gráfico, el valor aleatorio uniforme (0,1) R1,
se proyecta desde la escala vertical F(x) para obtener el valor de muestra x1 en la escala horizontal. Esto tanto para la distribución aleatoria continua como discreta.
El
método resulta válido dado que la variable aleatoria z = F(x) está
uniformemente distribuida en el intervalo 0 ≤ z ≤ 1 . 14
MTI : Distribución exponencial El
tiempo t entre llegadas de clientes a una instalación se representa con una distribución exponencial con media E{t} = 1/λ unidades de tiempo, cuya función de densidad de probabilidad es:
f(t) = λe-λt , t > 0 Como
se mencionó, el objetivo es determinar una muestra aleatoria t a partir de f(t). La función de densidad acumulada es: 1
F(t) = ∫ λe -λx dx = 0
Al
1 – e -λt
, t>0
igualar R = F(t) se despeja t y se obtiene:
t = – (1 / λ). ln(1 – R)
(1)
1 – R es el complemento de R, que es un número aleatorio entre 0 y 1, (1 – R) también lo es. Por lo tanto, puede reemplazar ln(1 – R) por ln(R) Como
15
MTI : Distribución exponencial En
términos de simulación, el resultado indica que las llegas están distanciadas t unidades de tiempo. ejemplo, si λ = 4 clientes / hora, y R = 0,9 , el tiempo que transcurre hasta la próxima hasta la siguiente llegada de un cliente se calcula: Por
t1 = – (1 / 4). ln(1 – 0,9) =
0,576 horas ~ 34,5 minutos
que los números aleatorios R, R1, R2, … usados para obtener muestras sucesivas se obtienen aleatoriamente de una distribución uniforme (0,1) a través de distintos métodos (tablas de números aleatorios, métodos aritméticos – por ejemplo, el congruente multiplicativo – , o simplemente en excel – función “ALEATORIO()” –). Recordar
16
Método de convolución En
este método se logra expresar la muestra deseada como la suma
estadística (convoluciones) de otras variables aleatorias más fáciles de muestrear.
Por
ejemplo, las muestras para distribuciones de Erlang y Poisson se
pueden obtener a partir de muestras de variables aleatorias exponenciales independientes e idénticamente distribuidas.
17
Distribución de Erlang Las
variables aleatorias Erlang y Poisson se definen como la suma estadística de m variables aleatorias exponenciales independientes e idénticamente distribuidas. Sea y la variable aleatoria de m-Erlang, entonces: y = y1 + y2 + … + ym donde yi , i = 1, 2, … m son variables aleatorias exponenciales independientes e idénticamente distribuidas cuya función de densidad se define como: En
f(yi) = λe-λy , yi > 0 , i = 1, 2, … m Y
por lo tanto, la i-ésima muestra de distribución exponencial sería: yi = – (1 / λ). ln(Ri)
18
, i = 1, 2, … m
Distribución de Erlang Entonces
la muestra de Erlang se calcula de la siguiente manera: y = – (1 / λ). {ln(R1) + ln(R2) + … + ln(Rm)} y = – (1 / λ). ln(R1 R2 … Rm)}
Por
(2)
ejemplo, si m = 3 y λ = 4 eventos por hora, obtenidos 3 números aleatorios se obtiene el resultado: (0,0589)(0,6733)(0,4799) = 0,0190 y luego y = – (1 / 4). ln(0,019) = 0,991 hora
19
MTI : Distribución Poisson Como
vimos al analizar la teoría de colas, si la distribución del tiempo entre la ocurrencia de eventos sucesivos es exponencial, la distribución de la cantidad de eventos por unidad de tiempo debe ser Poisson. supone que la distribución de Poisson tiene un valor medio de λ eventos por unidad de tiempo. Entonces, el tiempo entre eventos es exponencial, con una media de 1/λ unidades de tiempo. Se
Esto
quiere decir que una muestra de Poisson n se efectuara durante t unidades de tiempo si y sólo si: Período hasta que suceda n ≤ t < Período hasta que suceda n + 1
t1 + t2 + t3 + … + tn ≤ t < t 1 + t2 + t3 + … + tn+1 , n > 0 0 ≤ t 0
0 ≤ t ≤ – (1 / λ). ln(R1)
,n =0
Lo que se reduce a:
R 1R 2 … R n
≥ e– λt ≥ R1R2 … Rn+1
0 ≥ e– λt ≥ R1
21
,n>0 ,n =0
MTI : Distribución Poisson Un
ejemplo de la implementación del proceso de muestreo, supongamos que λ = 4 eventos por hora y que se quiere obtener una muestra para un período t = 0,5 horas. El
resultado es:
e– λt = e– 2 = 0,1353
Siendo
el número aleatorio (obtenido por alguno de los métodos mencionados) R1 = 0,0852, se obtiene que:
e– λt > R1 Por
22
lo que la muestra corresponde a n = 0.
Simulación (Ejemplo TAHA) El
tiempo de llegada de clientes a una peluquería tiene una media de 15 minutos (distribución exponencial). En el local hay un solo peluquero y tarda de 10 a 15 minutos (distribución uniforme) para terminar un corte de pelo. Los clientes son atendidos de acuerdo al orden de llegada (FIFO). El
objetivo de la simulación es calcular las siguientes medidas de desempeño: Nivel
de ocupación de la instalación de servicio (peluquero).
Cantidad
Tiempo
Como
promedio de clientes en espera
de espera promedio de un cliente en la cola
ya se mencionó, la lógica del modelo de simulación de eventos discretos se puede describir en términos de las acciones asociadas con sus eventos de llegada y eventos de salida.
23
Simulación (Ejemplo TAHA) Evento
de llegada - Pasos
1.
Se genera y guarda cronológicamente la hora de llegada del siguiente cliente ( = hora de simulación actual + tiempo entre llegadas).
2.
Si la instalación (peluquero) está inactiva:
3.
24
A.
Se inicia el servicio y se declara ocupada la instalación. Se actualizan las estadísticas de utilización de la instalación.
B.
Se genera y guarda cronológicamente la hora de salida del cliente ( = hora de simulación actual + tiempo de servicio).
Si la instalación está ocupada el cliente ingresa a la fila de espera y se actualiza la estadística de la cola.
Simulación (Ejemplo TAHA) Evento
de salida - Pasos
1.
Si la cola está vacía se declara inactiva la instalación. Se actualizan las estadísticas de utilización de la instalación.
2.
Si la cola no esta vacía:
25
A.
Se ingresa un cliente de la cola a la instalación. Se actualizan las estadísticas de la cola y de utilización de la instalación.
B.
Se genera y guarda cronológicamente la hora de salida del cliente ( = hora de simulación actual + tiempo de servicio).
Simulación (Ejemplo TAHA) Según
los datos del problema el tiempo entre llegadas tiene distribución exponencial y media de 15 minutos y el tiempo de servicio una distribución uniforme entre 10 y 15 minutos. Si
p y q representan muestras aleatorias de tiempos de llegada y servicio, entonces se obtiene: p = – 15 ln(R) minutos, 0 ≤ R ≤ 1 q = 10 + 5R minutos, 0 ≤ R ≤ 1 Vamos
a representar como T la simulación del tiempo y vamos a suponer que el primer cliente llega en T = 0 y que la instalación comienza vacía.
26
Simulación (Ejemplo TAHA) Llegada T
= 0.
Llegada T
del cliente 2
= 0 + p1 = 0 + [ –15 ln(0,0589)] = 42,48 minutos
Salida T
del cliente 1
del cliente 1
= 0 + q1 = 0 + (10 + 5 . 0,6733) = 13,37 minutos
Entonces
el cliente 2 ingresa a servicio inmediatamente, dado que la instalación estuvo ociosa entre los 13,37 y 42,48 minutos. Llegada T
27
del cliente 3
= 42,48 + p2 = 42,48 + [ –15 ln(0,4799)] = 53,49 minutos
Simulación (Ejemplo TAHA) Salida T
del cliente 2
= 42,48 + q2 = 42,48 + (10 + 5 . 0,9486) = 57,22 minutos
lo tanto, el cliente 3 debió esperar en la cola (57,22 – 53,49 = 3,73) minutos para ser atendido. Por
Llegada T
= 53,49 + p3 = 53,49 + [ –15 ln(0,6139)] = 60,81 minutos
Salida T
del cliente 4
del cliente 3
= 57,22 + q3 = 57,22 + (10 + 5 . 0,5933) = 70,19 minutos
lo tanto, el cliente 4 debió esperar en la cola (70,19 – 60,81 = 9,38) minutos para ser atendido. Por
28
Simulación (Ejemplo TAHA) Llegada T
= 60,81 + p4 = 60,81 + [ –15 ln(0,9341)] = 61,83 minutos
Salida T
del cliente 5
del cliente 4
= 70,19 + q4 = 70,19 + (10 + 5 . 0,1782) = 81,08 minutos
lo tanto, el cliente 5 debió esperar en la cola (81,08 – 61,83 = 19,25) minutos para ser atendido. Por
Salida T
29
del cliente 5
= 81,08 + q5 = 81,08 + (10 + 5 . 0,3473) = 92,82 minutos
Simulación (Ejemplo TAHA) En
los gráficos se muestra un resumen de los cambios en la variación de la cola y de la utilización de la instalación, en función del tiempo de simulación.
30
Simulación (Ejemplo TAHA) De
los datos obtenidos
tenemos que: Tiempo
de espera en la cola:
3,73 + 9,38 + 19,25 = 32,36 Tiempo
promedio de espera:
32,36 / 5 = 6,47 Longitud
promedio de la cola:
32,36 / 92,82 = 0,349 clientes Utilización
promedio de la instalación:
(13,37 + 54,34) / 92,82 = 0,686 peluquero
31