Modelos y Sumilación

Método de Monte Carlo

¿Qué es el método Monte Carlo?

El método de Montecarlo es un método no determinativo o estadístico numérico, usado para aproximar expresiones matemáticas complejas y costosas de evaluar con exactitud.El método se llamó así en referencia al Casino de Monte Carlo (Principado de Mónaco) por ser “la capital del juego de azar”, al ser la ruleta un generador simple de números aleatorios. El nombre y el desarrollo sistemático de los métodos de Monte Carlo datan aproximadamente de 1944 y se mejoraron enormemente con el desarrollo de la computadora.

(http://es.wikipedia.org/wiki/M%C3%A9todo_de_Montecarlo)

¿Cual es su origen?

La invención del método de Monte Carlo se asigna a Stanislaw Ulam y a John von Neumann. Ulam ha explicado cómo se le ocurrió la idea mientras jugaba un solitario durante una enfermedad en 1946. Advirtió que resulta mucho más simple tener una idea del resultado general del solitario haciendo pruebas múltiples con las cartas y contando las proporciones de los resultados que computar todas las posibilidades de combinación formalmente. Se le ocurrió que esta misma observación debía aplicarse a su trabajo de Los Álamos sobre difusión de neutrones, para la cual resulta prácticamente imposible solucionar las ecuaciones íntegro-diferenciales que gobiernan la dispersión, la absorción y la fisión.

(http://es.wikipedia.org/wiki/M%C3%A9todo_de_Montecarlo)

Técnicas de validación estadística

En cualquier trabajo de modelado y muy especialmente en simulación, el proceso de plantear hipótesis, construir el modelo y validarlo es un proceso cíclico. Muchas veces las partes individuales de un modelo parecen representar la realidad, pero cuando se consideran en conjunto, resulta en un pobre reflejo de la conducta del sistema en general.

Existen dos formas de validar un modelo:

a) Permitir que el usuario chequee que la simulación se desarrolla como debe. El usuario no tiene por qué entender el código, pero sí debe poder entender el diagrama de actividades y debe participar activamente en el planteo de los objetivos del trabajo y por ende en la lógica y detalles de la simulación. Es importante brindarle resultados visuales del comportamiento de las colas, entidades y el uso de recursos, que le permitan ver si la simulación se comporta en forma similar al sistema real.

b) Brindar estadísticas que confirmen que la simulación produce resultados similares a los del sistema real. Esto necesita de una recolección de datos adicional acerca de promedios de largos de colas, tiempo de ocupación de los servidores y tiempos de espera, los que se confrontarán con los obtenidos mediante la simulación.

(http://www.fing.edu.uy/inco/cursos/simulacion/archivos/clases/clase11web.pdf)

Cadenas de Markov

Una cadena de markov consta de unos estados E1 E2 E3 E4…..En. que inicialmente en un tiempo 0 o paso 0 se le llama estado inicial, además de esto consta de una matriz de transición que significa la posibilidad de que se cambie de estado en un próximo tiempo o paso.

Matriz de transmisión:

Una matriz de transición para una cadena de Markov de n estado es una matriz de n X n con todos los registros no negativos y con la propiedad adicional de que la suma de los registros de cada columna (o fila) es 1. Por ejemplo: las siguientes son matrices de transición.

Ejemplo:

En un país como Colombia existen 3 operadores principales de telefonía móvil como lo son tigo, Comcel y movistar (estados).

Los porcentajes actuales que tiene cada operador en el mercado actual son para tigo 0.4 para Comcel 0.25 y para movistar 0.35. (estado inicial)

Se tiene la siguiente información un usuario actualmente de tigo tiene una probabilidad de permanecer en tigo de 0.60, de pasar a Comcel 0.2 y de pasarse a movistar de 0.2; si en la actualidad el usuario es cliente de Comcel tiene una probabilidad de mantenerse en Comcel del 0.5 de que esta persona se cambie a tigo 0.3 y que se pase a movistar de 0.2; si el usuario es cliente en la actualidad de movistar la probabilidad que permanezca en movistar es de 0.4 de que se cambie a tigo de 0.3 y a Comcel de 0.3.
Partiendo de esta información podemos elaborar la matriz de transición.

La suma de las probabilidades de cada estado en este caso operador deben ser iguales a 1

Po= (0.4 0.25 0.35) → estado inicial

También se puede mostrar la transición por un método grafico

Ahora procedemos a encontrar los estados en los siguientes pasos o tiempos, esto se realiza multiplicando la matriz de transición por el estado inicial y asi sucesivamente pero multiplicando por el estado inmediatamente anterior.

Como podemos ver la variación en el periodo 4 al 5 es muy mínima casi insignificante podemos decir que ya se a llegado al vector o estado estable.
(http://investigaciondeoperaciones2markov.blogspot.com/p/teoria-y-ejemplos.html)

Algoritmos de Metrópolis-Hastings

Este algoritmo construye una cadena de Markov apropiada definiendo las probabilidades de transición de la siguiente manera.

Sea $Q(\bmath{\theta}^* \vert \bmath{\theta})$ una distribución de transición (arbitraria) y definamos

$\begin{displaymath} \alpha(\bmath{\theta}^*, \bmath{\theta}) = \min \left\{ \fra... ...x}) \, Q(\bmath{\theta}^* \vert \bmath{\theta})}, 1 \right\}. \end{displaymath}$

Algoritmo. Dado un valor inicial $\bmath{\theta}^{(0)}$ , la -ésima iteración consiste en:

1. generar una observación $\bmath{\theta}^*$ de $Q(\bmath{\theta}^* \vert \bmath{\theta}^{(t)})$ ;

2. generar una variable $u \sim U(0,1)$ ;

3. si $u \leq \alpha(\bmath{\theta}^*, \bmath{\theta}^{(t)})$ , hacer $\bmath{\theta}^{(t+1)} = \bmath{\theta}^*$ ; en caso contrario, hacer $\bmath{\theta}^{(t+1)} = \bmath{\theta}^{(t)}$ .

Este procedimiento genera una cadena de Markov con distribución de transición

$\begin{displaymath} P(\bmath{\theta}^{(t+1)} \vert \bmath{\theta}^{(t)}) = \alp... ...t)}) \, Q(\bmath{\theta}^{(t+1)} \vert \bmath{\theta}^{(t)}). \end{displaymath}$

La probabilidad de aceptación $\alpha(\bmath{\theta}^*, \bmath{\theta})$ sólo depende de $p(\bmath{\theta} \vert \bmath{x})$ a través de un cociente, por lo que la constante de normalización no es necesaria.

Comentario. La versión original del algoritmo de Metropolis toma $Q(\bmath{\theta}^* \vert \bmath{\theta}) = Q(\bmath{\theta} \vert \bmath{\theta}^*)$ , en cuyo caso

$\begin{displaymath} \alpha(\bmath{\theta}^*, \bmath{\theta}) = \min \left\{ \fra... ...rt \bmath{x})}{p(\bmath{\theta} \vert \bmath{x})}, 1 \right\}. \end{displaymath}$

Dos casos particulares utilizados comúnmente en la práctica son:

$\diamond$ Caminata aleatoria. Sea $Q(\bmath{\theta}^* \vert \bmath{\theta}) = Q_1(\bmath{\theta}^* - \bmath{\theta})$ , donde $Q_1(\cdot)$ es una densidad de probabilidad simétrica centrada en el origen. Entonces

$\begin{displaymath} \alpha(\bmath{\theta}^*, \bmath{\theta}) = \min \left\{ \fra... ...rt \bmath{x})}{p(\bmath{\theta} \vert \bmath{x})}, 1 \right\}. \end{displaymath}$

$\diamond$ Independencia. Sea $Q(\bmath{\theta}^* \vert \bmath{\theta}) = Q_0(\bmath{\theta}^*)$ , donde $Q_0(\cdot)$ es una densidad de probabilidad sobre $\Theta$ . Entonces

$\begin{displaymath} \alpha(\bmath{\theta}^*, \bmath{\theta}) = \min \left\{ \fr... ...omega(\bmath{\theta}^*)}{\omega(\bmath{\theta})}, 1 \right\}, \end{displaymath}$

con $\omega(\bmath{\theta}) = p(\bmath{\theta} \vert \bmath{x}) / Q_0(\bmath{\theta})$ .

En la práctica es común utilizar, después de una reparametrización apropiada, distribuciones de transición normales ó

de Student ligeramente sobredispersas, e.g.

$\begin{displaymath} Q(\bmath{\theta}^* \vert \bmath{\theta}) = N_d(\bmath{\theta... ...eta}})) \; \; \; \; \; \; \; \; \; \mbox{(caminata aleatoria)} \end{displaymath}$

$\begin{displaymath} Q_0(\bmath{\theta}^*) = N_d(\bmath{\theta}^* \vert \hat{\bm... ...{\theta}})) \; \; \; \; \; \; \; \; \; \mbox{(independencia)}, \end{displaymath}$

donde $\hat{\bmath{\theta}}$ y $\bmath{V}(\hat{\bmath{\theta}})$ denotan a la media y a la matriz de varianzas-covarianzas de la aproximación normal asintótica para $p(\bmath{\theta} \vert \bmath{x})$ , respectivamente, y $\kappa \geq 1$ es un factor de sobredispersión.

$\begin{Example} Consideremos el modelo jer\'arquico definido por \begin{displaym... ...bmath{\omega}) \propto p_0(\phi) \, p(\bmath{\omega} \vert \phi)$. \end{Example}$

(http://www.dpye.iimas.unam.mx/eduardo/MCB/node26.html)

Muestreo de GIBBS

Al igual que el algoritmo de Metropolis, el algoritmo de Gibbs permite simular una cadena de Markov $\bmath{\theta}^{(1)},\bmath{\theta}^{(2)},\ldots$ con distribución de equilibrio $p(\bmath{\theta} \vert \bmath{x})$ . En este caso, sin embargo, cada nuevo valor de la cadena se obtiene a través de un proceso iterativo que sólo requiere generar muestras de distribuciones cuya dimensión es menor que

y que en la mayoría de los casos tienen una forma más sencilla que la de $p(\bmath{\theta} \vert \bmath{x})$ .

Sea $\bmath{\theta} = (\bmath{\theta}_1,\ldots,\bmath{\theta}_k)$ una partición del vector $\bmath{\theta}$ , donde $\bmath{\theta}_i \in \Rex^{d_i}$ y $\sum_{i=1}^{k} d_i = d$ . Las densidades

$\begin{displaymath} \begin{array}{cr} p(\bmath{\theta}_1 \vert \bmath{\theta}_2,... ...theta}_1,\ldots,\bmath{\theta}_{k-1}, \bmath{x}) & \end{array}\end{displaymath}$

se conocen como densidades condicionales completas y en general pueden identificarse fácilmente al inspeccionar la forma de la distribución final $p(\bmath{\theta} \vert \bmath{x})$ . De hecho, para cada $i=1,\ldots,k$ ,

$\begin{displaymath} p(\bmath{\theta}_i \vert \bmath{\theta}_1,\ldots,\bmath{\the... ...}_{k}, \bmath{x}) \propto p(\bmath{\theta} \vert \bmath{x}), \end{displaymath}$

donde $p(\bmath{\theta} \vert \bmath{x}) = p(\bmath{\theta}_1,\ldots,\bmath{\theta}_k \vert \bmath{x})$ es vista sólo como función de $\bmath{\theta}_i$ .

Dado un valor inicial $\bmath{\theta}^{(0)} = (\bmath{\theta}_1^{(0)},\ldots,\bmath{\theta}_k^{(0)})$ , el algoritmo de Gibbs simula una cadena de Markov en la que $\bmath{\theta}^{(t+1)}$ se obtiene a partir de $\bmath{\theta}^{(t)}$ de la siguiente manera:

generar una observación $\bmath{\theta}_1^{(t+1)}$ de $p(\bmath{\theta}_1 \vert \bmath{\theta}_2^{(t)},\bmath{\theta}_3^{(t)},\ldots,\bmath{\theta}_k^{(t)}, \bmath{x})$ ;

generar una observación $\bmath{\theta}_2^{(t+1)}$ de $p(\bmath{\theta}_2 \vert \bmath{\theta}_1^{(t+1)},\bmath{\theta}_3^{(t)},\ldots,\bmath{\theta}_k^{(t)}, \bmath{x})$ ;

$\vdots$

generar una observación $\bmath{\theta}_k^{(t+1)}$ de $p(\bmath{\theta}_k \vert \bmath{\theta}_1^{(t+1)},\bmath{\theta}_2^{(t+1)},\ldots,\bmath{\theta}_{k-1}^{(t+1)}, \bmath{x})$ .

La sucesión $\bmath{\theta}^{(1)},\bmath{\theta}^{(2)},\ldots$ así obtenida es entonces una realización de una cadena de Markov cuya distribución de transición está dada por

$\begin{displaymath} P(\bmath{\theta}^{(t+1)} \vert \bmath{\theta}^{(t)}) = \prod... ...eta}_{i+1}^{(t)},\ldots, \bmath{\theta}_{k}^{(t)}, \bmath{x}). \end{displaymath}$

Comentario. En ocasiones la distribución final implica cierta estructura de independencia condicional entre algunos de los elementos del vector $\bmath{\theta}$ . Es estos casos es común que muchas de las densidades condicionales completas se simplifiquen.

$\begin{Example} % latex2html id marker 3099Sea $\bmath{x} = (\bmath{z},\tilde{... ...ert \bmath{z})$, es decir, de la distribuci\'on que nos interesa. \end{Example}$

(http://www.dpye.iimas.unam.mx/eduardo/MCB/node27.html)

Algoritmo Templado Simulado

Simulated annealing (SA) o recocido simulado es un algoritmo de búsqueda meta-heurística para problemas de optimización global; el objetivo general de este tipo de algoritmos es encontrar una buena aproximación al valor óptimo de una función en un espacio de búsqueda grande. A este valor óptimo se lo denomina"óptimo global".

El método fue descrito independientemente por Scott Kirkpatrick, C. Daniel Gelatt y Mario P. Vecchi en 1983 y por Vlado Černý en 1985. El método es una adaptación del algoritmo Metropolis-Hastings, unmétodo deMontecarloutilizado para generar muestras de estados de unsistema termodinámico.

ALGORITMO DE TEMPLADO SIMULADO

Algoritmo de Wang-Landau

El algoritmo Wang-Landau es una extensión del método de Monte Carlo, propuesto por Fugao Wang y David P. Landau, que permite calcular la densidad de estados de un sistema dado sin tener ningún conocimiento previo sobre ella. Se dice entonces que la densidad de estados se calcula on the fly, o sea durante la simulación. El algoritmo es independiente de la temperatura (un problema común con otros algoritmos Monte Carlo) y es fácil de implementar.

Este algoritmo fue pensado inicialmente para muestrear el espacio de configuraciones de ciertos sistemas que no podían ser tratados mediante el algoritmo de Metropolis. El algoritmo de Metropolis realiza un muestreo usando el factor de Boltzmann para aceptar o descartar configuraciones. Metropolis funciona bien si todas las configuraciones posibles del sistema se encuentran dentro de un rango relativamente angosto de energías.

Modelos y Sumilación

jueves, 15 de agosto de 2013

Simulación discreta

martes, 13 de agosto de 2013

miércoles, 7 de agosto de 2013

Cuestionario

lunes, 5 de agosto de 2013

Material de ayuda

jueves, 1 de agosto de 2013