\documentclass[a4paper,10pt]{article}
\usepackage[latin1]{inputenc}
\usepackage[spanish]{babel}
\usepackage[dvips]{epsfig}
\usepackage{mathrsfs}
\usepackage{amsmath}
\usepackage{pst-all}
\usepackage{amssymb}
\usepackage{fullpage}
\usepackage{graphicx}


\begin{document}

\title{\textbf{Módulo Montecarlo\\
Software MES \\
CGE Consultores }}
\author{René González Leiva \\
\emph{aplikadito@gmail.com}\\ \\
Estudiante Ingeniería Civil Informática \\
Universidad Técnica Federico Santa María}
\date{Julio, 2008}
\maketitle

\section{Introducción}

La simulación Montecarlo básicamente permite obtener muestreos
aleatorios de situaciones reales basado en modelos de comportamiento
de las variables que participan en cierto problema.

La implementación a realizar n el software MES de la simulación
Montecarlo deberá permite asignar distintas distribuciones de
probabilidad a las variables:

\begin{itemize}
  \item Tiempo entre Fallas (TBF)
  \item Tiempo de Reparación (TTR)
  \item Costo de Reparación (CR)
\end{itemize}

Con estas variables, se calcularán los siguientes indicadores:
\begin{itemize}
  \item Disponibilidad
  \item Costo de Mantención
  \item Costo de la Falta
\end{itemize}

Los resultados serán probabilísticos, por lo que se utilizan
histogramas para presentar la información.

\section{Distribuciones en Módulo Montecarlo para MES}
El módulo Montecarlo para MES implementará las siguientes
distribuciones de probabilidad:

\begin{itemize}
  \item Exponencial $Exp(\lambda)$
  \item Weibull $Weibull(\alpha,\beta)$
  \item Normal $N(\mu,\sigma^2)$
  \item LogNormal $Ln(\mu,\sigma^2)$
  \item Triangular $Tr(a,b,c)$
  \item Uniforme $U(a,b)$
  \item Dirac $Dirac(p)$
\end{itemize}

Se definen varios métodos comunes a todas las distribuciones. Se
pueden destacar los siguientes métodos:
\begin{itemize}
  \item Un método que implemente la función de densidad de la
  distribución $f(x)$
  \item Un método que implemente la función de probabilidad acumulada
  de la distribución $F(x)$
  \item Un método que implemente la función inversa de la función de
  distribución de probabilidad acumulada $F^{-1}(x)$
\end{itemize}

\section{Obtención de números aleatorios}
Sea una variable X que se puede modelar con una exponencial de media
0,007, $Exp(0,007)$. Para obtener valores de la variable X acorde
con su comportamiento exponencial, se realiza lo siguiente:
\begin{enumerate}
  \item Se obtiene un número aleatorio uniformemente distribuido
  entre 0 y 1 (utilizando la clase \textbf{Random} y el método \textbf{nextDouble()})
  \item Luego, con ese número, se llama al método que obtiene la
  función inversa de la función de distribucion acumulada.
  \item El resultado de la inversa es un número dentro del rango
  definido por la distribución exponencial asociada.
\end{enumerate}

Por ejemplo, suponiendo que el número aleatorio obtenido fue $0.4$.
Sabemos que la función de distribución acumulada de una exponencial
es

$$ F(x) =  1 - e^{-(\lambda x)} $$

La inversa se puede obtener directamente

$$ F(x)^{-1} =  - \frac{e^{1-y}}{\lambda} $$

donde $y$ es un número entre 0 y 1. De forma genérica

$$ F(x)^{-1} =  - \frac{e^{1- U(0,1)}}{\lambda} $$

Ahora, simplemente se reemplazan los valores: $y = 0.4; \lambda =
0.007$

$$ F(x)^{-1} = 260.3 $$

Con esto, se realiza un mapeo de una función uniforme entre 0 y 1 a
una exponencial con media 0.007. Y en general, este procedimiento es
similar para toda distribución que posea una función inversa de
probabilidad acumulada definida. Este método se conoce como
\textbf{Algoritmo de la Transformada Inversa} y se puede definir, de
forma general, de la siguiente manera:\\

\emph{Sea U una variable aleatoria uniforme en (0,1). Para cualquier
función de distribución continua F, invertible, la variable
aleatoria X definida como }

$$X = F^ {-1}(U)$$

\emph{tiene distribución $F$ [$F^{-1}$ se define como el valor de x
tal que $F(x) = u$]}


\section{Obtención de números aleatorios para Normal y LogNormal}

\subsection{Usando método nextGaussian() de la clase Random}
Dado que no es posible aplicar el método anterior a las
distribuciones Normal y LogNormal, se debe utilizar un generador de
números aleatorios normalmente distribuido (utilizando la clase
\textbf{Random} y el método \textbf{nextGaussian()}) de manera
estándar. Para obtener un número aleatorio respecto a una
$Normal(\mu, \sigma^2)$, simplemente, se realiza un cambio de
variables

$$ n = N(0,1)$$
$$ z = \sigma \cdot n + \mu$$

donde $z$ es un número $N(\mu,\sigma)$-mente distribuido.

Para obtener un valor aleatorio asociado a una LogNormal con media
$\mu$ y desviación estandar $\sigma$, primero, utilizan las
siguientes variables auxiliares

$$ \mu_{aux} = \lg \bigg( \frac{\mu^2}{ \sqrt{ \mu^2 + \sigma^2} } \bigg) $$
$$ \sigma_{aux} = \sqrt{ \lg \Bigg( 1 + \bigg( \frac{\sigma}{\mu} \bigg)^2 \Bigg) } $$

Para después realizar el cambio de variables

$$z = e^{\sigma_{aux} \cdot n + \mu_{aux} }$$

donde $z$ es un número $Log(\mu,\sigma)$-mente distribuido. Este es
el método que actualmente se está usando.

\subsection{Método de búsqueda recursiva}
Este método no se está utilizando actualmente debido a que su
performance no es buena. No obstante, su implementanción no ha sido
borrada debido a que en el futuro eventualmente podría ser usada. Se
adjunta a continuación su explicación

Primero, se define una distribucion Normal Estandar, N(0,1) que
sirve de base para los calculos de las distribuciones de
probabilidad acumulada y de la inversa. Las razones se detallan a
continuación.

Una distribución Normal no posee una función de distribucion
acumulada definida. El cálculo para la obtención de esta función se
realiza utilizando métodos de integración numérica, o sea, que
utilizan la idea básica de Riemann, dividir el dominio de la función
en varios intervarlos y calcular el area de estos intervalos con
polígonos limitados en su altura por los valores de la función, ya
sea por exceso o por defecto. El método se Simpson es uno de los
mejores métodos de integracion numérica que hay. Crea polígonos un
poco mas complejos a lo largo de la curva y calcula con bastante
precisión el area bajo ella. La única dificultad de este método es
que necesita de los límites superior e inferior de integración, los
cuales no existen para el caso particular de una Normal. Para ello,
se utiliza la Normal Estandar, para la cual se conocen limites
superiores e inferiores que son factibles de utilizar. En este caso,
se utiliza como limite inferior -5 y limite superior 5, dado que las
probabilidades de que los números estén fuera del intervalo
\textbf{[-5,5]} es ínfima. Por lo tanto, se utilizan estos límites
para decir que una Normal Estandar tiene un dominio que sólo existe
dentro del intervalo ya definido. Ahora, es posible utilizar el
método de Simpson para integrar una Normal Estandar.

Si se desea obtener la función de probabilidad acumulada para una
$Normal(\mu,\sigma^2)$ para una valor x (es decir, se desea saber
cuanto vale la integral desde $-\infty$ hasta el punto x) se puede
utilizar el siguiente cambio de variable:

$$ z =  \frac{x - \mu}{\sigma} $$

donde z es un valor mapeado para el dominio de una Normal Estandar.
Luego, se calcula la integral de una Normal Estandar desde $-\infty$
(que en nuestro caso particular hemos definido como \textbf{-5})
hasta el punto z usando el método de Simpson.

De forma idéntica, para una LogNormal se utiliza:

$$z = \frac{ln(x) - \mu}{\sigma} $$

y se calcula la integral para la Normal Estandar usando z.

Finalmente, es posible realizar el cálculo de la función inversa de
la distribución acumulada usando un algoritmo recursivo. Para ello,
se debe realizar un mapeo inverso. Primero, supongamos que se tiene
un número \textbf{p} aleatorio uniformemente distribuido entre 0 y
1. Lo que deseamos saber es con que valor \textbf{z} en la función
de distribución acumulada de una Normal Estandar es posible obtener
el valor \textbf{p}. En otras palabras, encontrar el valor
\textbf{z} tal que la integral entre $-\infty$ y \textbf{z} para una
Normal Estandar sea igual a \textbf{p}.

Este valor \textbf{z} se obtiene realizando una búsqueda binaria
entre los limites de la Normal Estandar anteriormente definidos de
la siguiente manera:
\begin{enumerate}
  \item Se comienza con limites definidos ( [-5,5] )
  \item Se evalua en la mitad de los limites actuales ( F(0) )
  \item Si el valor objetivo \textbf{p} es mayor que la evaluación
  en el punto medio, se corre el limite inferior hasta el punto medio y se
  vuelve a evaluar en la mitad de los nuevos limites.
  \item Si el valor objetivo \textbf{p} es menor que la evaluación
  en el punto medio, se corre el límite superior hata el punto medio
  y se vuelve a evaluar en la mitad de los nuevos limites.
  \item El algorimo para cuando se la diferencia entre el punto
  evaluado y la meta \textbf{p} es menor a cierta precisión.
\end{enumerate}

Esto es análogo al proceso de búsqueda binaria de, por ejemplo, un
nombre en un directorio telefónico, buscar iterativamente en las
mitades de las mitades.

De esta forma, se obtiene un valor \textbf{z} que cumple las
condiciones ya mencionadas. Este valor se mapea para el dominio de
una $Normal(\mu,\sigma^2)$ de la siguiente manera:

$$x = \sigma*z + \mu$$

donde \textbf{x} es el valor mapeado para una $Normal(\mu,\sigma^2)$
respecto de \textbf{z} en una Normal Estandar. Para una
$LogNormal(\mu,\sigma^2)$ el cambio de variable es:

$$x = e^{\sigma * z + \mu}$$

donde \textbf{x} es el valor mapeado para una
$LogNormal(\mu,\sigma^2)$ respecto de \textbf{z} en una Normal
Estandar.

A pesar de tener definido y resuelto el problema de la obtención de
inversas, es necesario definir cotas inferiores y superiores para
las distribuciones, debido al problema computacional de graficar las
funciones y además del requerimiento para el módulo de que las
distribuciones puedan ser truncadas por cotas definidas por el
usuario. Se deben definir cotas superiores e inferiores por defecto.
Éstas, arbitrariamente, se definen de la siguiente manera:
\begin{itemize}
  \item \textbf{Cota Inferior}: un valor \textbf{x\_inf} tal que la probabilidad desde $-\infty$ hasta
  aquel punto sea \textbf{0,1\%}. O sea, la funcion inversa de probabilidad
  acumulada de la distribucion para el valor \textbf{0,001}.
  \item \textbf{Cota Superior}: el valor \textbf{x\_sup} tal que la probabilidad desde $-\infty$ hasta
  aquel punto sea \textbf{99,9\%}. O sea, la funcion inversa de probabilidad
  acumulada de la distribucion para el valor \textbf{0.999}.
\end{itemize}

\end{document}
