Simulación de Monte Carlo en Visual Basic

Simulación de Monte Carlo

Un poco de Teoria

La simulación de Monte Carlo es una técnica cuantitativa que hace uso de la estadística y los ordenadores para imitar, mediante modelos matemáticos, el comportamiento aleatorio de sistemas reales no dinámicos (por lo general, cuando se trata de sistemas cuyo estado va cambiando con el paso del tiempo, se recurre bien a la simulación de eventos discretos o bien a la simulación de sistemas continuos).
La clave de la simulación MC consiste en crear un modelo matemático del sistema, proceso o actividad que se quiere analizar, identificando aquellas variables (inputs del modelo) cuyo comportamiento aleatorio determina el comportamiento global del sistema. Una vez identificados dichos inputs o variables aleatorias, se lleva a cabo un experimento consistente en:

1. Generar –con ayuda del ordenador- muestras aleatorias (valores concretos) para dichos inputs, y
2. Analizar el comportamiento del sistema ante los valores generados. 

Tras repetir n veces este experimento, dispondremos de n observaciones sobre el comportamiento del sistema, lo cual nos será de utilidad para entender el funcionamiento del mismo –obviamente, nuestro análisis será tanto más preciso cuanto mayor sea el número n de experimentos que llevemos a cabo.
Veamos un ejemplo sencillo:
En la tabla inferior se muestra un análisis histórico de 200 días sobre el número de consultas diarias realizadas a un sistema de información empresarial (EIS) residente en un servidor central.
La tabla incluye el número de consultas diarias (0 a 5) junto con las frecuencias absolutas (número de días que se producen 0, 1, ..., 5 consultas), las frecuencias relativas (10/200 = 0,05, ...), y las frecuencias relativas acumuladas.

Consultas EIS	Frecuencias Absolutas (dias)	Frecuencia Relativa	Frecuencia Relativa Acumulada
0	10	0.05	0.05
1	20	0.10	0.15
2	40	0.20	0.35
3	60	0.30	0.65
4	40	0.20	0.85
5	30	0.15	1.00
Total	200	1.00

Podemos interpretar la frecuencia relativa como la probabilidad de que ocurra el suceso asociado, en este caso, la probabilidad de un determinado número de consultas (así, p.e., la probabilidad de que se den 3 consultas en un día sería de 0,30), por lo que la tabla anterior nos proporciona la distribución de probabilidad asociada a una variable aleatoria discreta (la variable aleatoria es el número de consultas al EIS, que sólo puede tomar valores enteros entre 0 y 5).
Supongamos que queremos conocer el número esperado (o medio) de consultas por día.
La respuesta a esta pregunta es fácil si recurrimos a la teoría de la probabilidad: Denotando por X a la variable aleatoria que representa el número diario de consultas al EIS,sabemos que: E[X]=0*0.05+1*0.10+...+5*0,15=2,95
Por otra parte, también podemos usar simulación de Monte Carlo para estimar el número esperado de consultas diarias (en este caso se ha podido obtener el valor exacto usando teoría de probabilidad, pero ello no siempre será factible). Veamos cómo:
Cuando se conozca la distribución de probabilidad asociada a una variable aleatoria discreta, será posible usar la columna de frecuencias relativas acumuladas para obtener los llamados intervalos de números aleatorios asociados a cada suceso. En este caso, los intervalos obtenidos son:

• [0.00,0.05) para el suceso 0 
• [0.05,0.15) para el suceso 1 
• [0.15,0.35) para el suceso 2 
• [0.35,0.65) para el suceso 3 
• [0.65,0.85) para el suceso 4 
• [0.85,1.00) para el suceso 5 

El gráfico siguiente nos muestra cada una de las probabilidades sobre el número de consultas. En él, se aprecia claramente la relación existente entre probabilidad de cada suceso y el área que éste ocupa.

Esto significa que, al generar un número pseudo-aleatorio con el ordenador (proveniente de una distribución uniforme entre 0 y 1), estaremos llevando a cabo un experimento cuyo resultado, obtenido de forma aleatoria y según la distribución de probabilidad anterior, estará asociado a un suceso. Así por ejemplo, si el ordenador nos proporciona el número pseudo-aleatorio 0.2567, podremos suponer que ese día se han producido 2 consultas al EIS. A continuación simularemos una serie de numeros pseudo-aleatorios y veremos cual es el resultado:

Aleatorio	Clase
0.74594	4
0.02769	0
0.41357	3
0.08281	1
0.92104	5
0.45407	3
0.79709	4
0.07949	1
0.48303	3
0.38276	3
0.99765	5
Promedio	2.909

En este caso, hemos obtenido un valor estimado que se acerca bastante con el valor real anteriormente calculado vía la definición teórica de la media. Sin embargo, debido a la componente aleatoria intrínseca al modelo, normalmente obtendremos valores “cercanos” al valor real, siendo dichos valores diferentes unos de otros (cada simulación proporcionará sus propios resultados). 

Si usamos pocas observaciones, los valores que obtendríamos no serían estimaciones tan buenas al valor real. Por el contrario, es de esperar que si hubiésemos usado 1000 (o mejor aún 10000) observaciones, los valores que obtendríamos en estarían todos muy cercanos al valor real.

Ahora implementemos esta Simulación

Primero veamos que necesitamos para implementar esta simulación:

1. Definir las "FRECUENCIAS" para cada "CLASE"
2. Generar los "RANGOS" en base a las "FRECUENCIAS ACUMULADAS" de las "CLASES"
3. En base a un "VALOR ALEATORIO" necesitamos evaluar en que "RANGO" cae, para así saber que "CLASE" le corresponde
4. Retornar como resultado la "CLASE" a la que le corresponde este "VALOR ALEATORIO"
5. Obtener el "VALOR ESPERADO" en base a una serie de "OBSERVACIONES"

Ahora, se han resaltado algunos puntos, y pasaremos a ver uno de ellos, el Rango. para eso definiremos la siguiente clase:

Public Class Rango
    Sub New(i As Nullable(Of Double), f As Nullable(Of Double))
        If i.HasValue Then Inicio = i
        If f.HasValue Then Final = f
    End Sub
    Property Inicio As Nullable(Of Double) = Double.NegativeInfinity
    Property Final As Nullable(Of Double) = Double.PositiveInfinity
    Public Function Dentro(v As Double) As Boolean
        Return Inicio <= v And v < Final
    End Function
    Public Overrides Function ToString() As String
        Return String.Format("De {0,6:N10} a {1,6:N10}", Inicio, Final)
    End Function
End Class

Si nos fijamos, en el constructor, pasamos el inicio y el final para el rango, pero usamos Nullable(Of Double) esto para permitir rangos orientados hacia el infinito positivo o negativo, según sea el caso.
Esto se comprueba en la definición de las propiedades
La función Dentro es la que evalúa si un valor esta dentro de este rango.
Por ultimo se ha sobreescrito la función ToString, con el fin de facilitar la depuración y ver el rango sin necesidad de expandir el la depuración

Ahora vamos a ver la Clase MonteCarlo:

Public Class Montercarlo
    Public Shared Property TamañoMuestra As Integer = 1000
    Private Shared _aleatorio As New Random
    Private _rangos As Dictionary(Of Integer, Rango)
    Private _frecuencias As Dictionary(Of Integer, Double)

    Public Sub New()
        _frecuencias = New Dictionary(Of Integer, Double)
    End Sub

    Public Sub AddFrecuencia(clase As Integer, frecuencia As Double)
        _frecuencias.Add(clase, frecuencia)
        CalcularRangos()
    End Sub

    Private Sub CalcularRangos()
        Dim _frecOrd = (From f In _frecuencias Order By f.Value Descending).ToList
        Dim _frecSum = _frecOrd.Sum(Function(k) k.Value)
        _rangos = New Dictionary(Of Integer, Rango)
        Dim v1 As Decimal = 0
        Dim v2 As Decimal = 0
        For Each f In _frecOrd
            v1 = v2
            v2 += f.Value / _frecSum
            _rangos.Add(f.Key, New Rango(v1, v2))
        Next
    End Sub

    Public Function GetUno() As Integer
        Return GetId(GetAleatorio())
    End Function

    Private Function GetId(v As Double) As Integer
        Return (From r In _rangos Where r.Value.Dentro(v)).FirstOrDefault.Key
    End Function

    Public Function GetEsperado() As Integer
        Dim muestra As New List(Of Double)
        For i = 1 To TamañoMuestra
            muestra.Add(GetUno())
        Next
        Return Aggregate a In muestra Into Average(a)
    End Function

    Private Function GetAleatorio() As Double
        Return _aleatorio.NextDouble)
    End Function
End Class

Ahora, tenemos varias cosas que resaltar:

Private _rangos As Dictionary(Of Integer, Rango)
Private _frecuencias As Dictionary(Of Integer, Double)

Aqui tenemos dos diccionarios privados, uno para los Rangos y otro para las Frecuencias. Y por que diccionarios? pues porque facilita la busqueda por Clave. La idea es que conforme vamos ingresando Frecuencias, recalcularemos los rangos, ya que no se sabe a priori los rangos. Por eso se implemento la funcion AddFrecuncia:

Public Sub AddFrecuencia(clase As Integer, frecuencia As Double)
    _frecuencias.Add(clave, frecuencia)
    CalcularRangos()
End Sub

Como se aprecia, pasamos la clase y la frecuencia, y la clase actua como Key en el diccionario y la frecuencia como Value. Inmediatamante añadida la frecuencia recalculamos los rangos

    Private Sub CalcularRangos()
        Dim _frecOrd = (From f In _frecuencias Order By f.Value Descending).ToList
        Dim _frecSum = _frecOrd.Sum(Function(k) k.Value)
        _rangos = New Dictionary(Of Integer, Rango)
        Dim v1 As Decimal = 0
        Dim v2 As Decimal = 0
        For Each f In _frecOrd
            v1 = v2
            v2 += f.Value / _frecSum
            _rangos.Add(f.Key, New Rango(v1, v2))
        Next
    End Sub

Aqui calculamos los rangos. si se fijan hay una variable _frecord, que contendra las frecuncias ordendas; en base a una consulta "LinQ to Objects". Ordenamos las frecuencias como una optimización del Metodo Monte Carlo para algunas aplicaciones (veremos para que en un post futuro). Tambien Hay una variable _frecsum, que calcula, en base a una consulta LinQ la sumatoria de las frecuencias, para poder calcular luego las frecuencias acumuladas, requisito indispensable para calcular los rangos. La sumatoria se esta haciendo con una función Lambda Bueno, luego hay un simple bucle donde vamos agregando objetos "Rango" con un Key igual a la clase. El resultado es el siguiente:

Como vemos, ya estan generados los rangos, pero cabe destacar que estos mismos estan reordenados en forma descendente en base a su frecuencia, esto ayuda en casos en que se necesita aplicar un modificador a la variable aleatoria y asi dar prioridad a los mas frecuentes, como dije, esto lo veremos en otro post.

Ahora veamos como generar numeros aleatorios:

Private Shared _aleatorio As New Random
Private Function GetAleatorio() As Double
    Return _aleatorio.NextDouble)
End Function

Primero definimos una variable estatica del tipo Random, y nos valemos de su método NextDouble, que nos devuelve un aleatorio entre 0 y 1 Y ahora veremos como obtener una "Observacion":

    Public Function GetUno() As Integer
        Return GetId(GetAleatorio())
    End Function

    Private Function GetId(v As Double) As Integer
        Return (From r In _rangos Where r.Value.Dentro(v)).FirstOrDefault.Key
    End Function

La funcion "GetUno" llama a la funcion "GetId" pasandole un aleatorio.
La funcion "GetId" busca en el diccionario de rangos, usando una consulta LinQ, uno donde el valor aleatorio se encuentre "dentro" y no devuelve el "key" del primer resultado
Ahora ya sabemos como obtener una observación, solo nos faltaria obtener el "valor esperado", y esto se hace con:

    Public Shared Property TamañoMuestra As Integer = 1000
    Public Function GetEsperado() As Integer
        Dim muestra As New List(Of Double)
        For i = 1 To TamañoMuestra
            muestra.Add(GetUno())
        Next
        Return Aggregate a In muestra Into Average(a)
    End Function

Como vemos, hemos presupuesto un tamaño demuestra de 1000 observaciones, y simplemente con un bucle vamos añadiendo muestras a una lista, para luego, mediante una consulta LinQ de agregación calculamos la media y retornamos este valor.
Ahora simulemos el caso que vimos en la parte teorica:

Imports GaSKSoft.Predictor.Montecarlo
Module Module1

    Sub Main()
        Dim m As New Montercarlo
        m = New Montercarlo
        m.AddFrecuencia(0, 10)
        m.AddFrecuencia(1, 20)
        m.AddFrecuencia(2, 40)
        m.AddFrecuencia(3, 60)
        m.AddFrecuencia(4, 40)
        m.AddFrecuencia(5, 30)
        Console.WriteLine("Frecuencias:")
        Console.WriteLine("0: 10")
        Console.WriteLine("1: 20")
        Console.WriteLine("2: 40")
        Console.WriteLine("3: 60")
        Console.WriteLine("4: 40")
        Console.WriteLine("5: 30")
        Console.WriteLine(
            "El valor esperado para una muestra de {0} observaciones es:{1}",
            Montercarlo.TamañoMuestra,
            m.GetEsperado())
        Console.ReadKey()
    End Sub

End Module

Y el resultado es:

Como vemos, con una muestra de 1000 observaciones, el resultado es mas preciso.

Concluciones

El metodo de Monte Carlo nos ayuda a simular situaciones en las que las probabilidades de varios eventos son conocidos, y se requiere obtener, en base a estas, un valor esperado, tomando como base una "nube" de observaciones.

En este post tambien se ha utilizado Consultas LinQ to Objects, y demuestra que reduce significativamente el proceso de consultar listas, incluso se le pueden calcular agregados del tipo suma, promedio, producto, cuenta, etc. en una sola linea!!!!!

Finalizando

Espero que el post les guste y les sirva para alguna implementación futura.