viernes, 3 de mayo de 2013

Viajando a otras dimensiones con Monte Carlo

Cálculo de PI mediante el método Monte-Carlo

Pero ¿se puede afinar más en la búsqueda del número PI por este método?
 La respuesta es un SÍ rotundo.
 Una característica importante a la hora de ejecutar un método Monte Carlo es realizar una buena elección del generador de números aleatorios, cosa que no hice. El cálculo de PI es un ejemplo de integración mediante método Monte Carlo: calculamos el área del sector del círculo y, gracias a la relación entre el área de ese sector y del cuadrado que lo contiene, obtenemos la aproximación de PI.
El área del cuadrado será
A_{cuadrado} = r^2
y el área del sector será
A_{sector} = 1/4 \pi r^2
Si dividimos el área del círculo por el área del sector, tendremos:
\frac{A_{sector}}{A_{cuadrado}} = \frac{\pi r^2}{4 r^2} = \pi/4
Si realizamos la división del número total de puntos que tenemos dentro del círculo respecto al número total de puntos que se han señalado en toda la área del cuadrado, tendremos una aproximación del valor de π
\frac{puntos dentro circulo}{puntos totales} \approx \pi/4
En este tipo de problemas, un buen generador de números aleatorios será aquel que consiga cubrir el espacio estudiado lo más homogéneamente posible. 
De esta forma, habrá más probabilidad de que la razón entre el total de puntos lanzados y los puntos dentro de nuestro objeto sea similar a la razón real entre el volumen de nuestro campo de pruebas y el objeto circunscrito que deseamos medir.
Investigando el tema, descubrí un generador de números pseudoaleatorios llamado secuencia de Sobol (más información aquí) y que distribuye puntos de una forma muy homogénea a lo largo de las dimensiones que elijamos. A continuación, un ejemplo de una distribución de puntos obtenida uniformemente y otra obtenida mediante una secuencia Sobol:
Comparativa de la generación de números aleatorios mediante el generador uniforme de python y la secuencia Sobol
Comparativa de la generación de 1000 puntos aleatorios mediante el generador uniforme de python y la secuencia Sobol
Lo siguiente fue comprobar si realmente se notaba la diferencia al realizar el cálculo de PI usando un método u otro. Aquí tenéis la comparación del cálculo de PI con un generador uniforme de números aleatorios y con un generador de secuencias Sobol:
Cálculo de pi mediante distintos tipos de generadores de números aleatorios.
Cálculo de pi mediante distintos tipos de generadores de números aleatorios.
El cálculo se realizó para distintos números de puntos lanzados, hasta un máximo de 100.000 puntos. Se puede ver con claridad la mejora que produce el uso de una secuencia Sobol (azul) respecto a un generador de números aleatorios uniforme. 
La convergencia hacia el valor de PI es muy evidente cuando usamos secuencias Sobol.
Para ver el código fuente usado puedes ir a la parte final.

¿Por qué conformarnos solo con dos dimensiones?

Grosso modo lo que nos cuenta es el particular comportamiento del volumen de una esfera de radio 1 a medida que aumentamos el número de dimensiones. 
Podemos ver cómo en principio el volumen aumenta con el número de dimensiones para luego decrecer, tendiendo hacia cero a medida que aumentamos n. 
La fórmula para el volumen de una n-esfera de radio 1 es:
V = \frac{\pi^{n/2}}{\Gamma \left ( \frac{n}{2}+1 \right )}
La explicación es bastante lógica; para que un punto del espacio de n-dimensiones se encuentre dentro de la n-esfera de radio 1 correspondiente, debe cumplir que su distancia al origen sea menor o igual a 1, y a medida que aumentamos las dimensiones llega un momento en el que es muy difícil encontrar puntos que cumplan esta condición, ya que para cada dimensión tenemos que sumar el cuadrado de un nuevo término.
La fórmula para saber si un punto está dentro de una n-esfera es:
\sum_{i=1}^n x_i^2 \leq 1
La rabia del asunto es que es imposible visualizar una n-esfera de más de tres dimensiones dentro de un cubo, a su vez, de n-dimensiones. 
Quizás sea imposible imaginarlo, pero una computadora sí que puede ayudarnos a hacer algo similar: podemos calcular aproximadamente el volumen de una n-esfera de la misma forma que hicimos antes con el área de un círculo.
 Así que toca lanzar n-puntos a cascoporro para meter el dedo en la yaga y ver si todo esto es cierto.
Para cada dimensión, el cálculo del volumen de la n-esfera de radio 1 será el siguiente:
V_n = V_{n cubo} \frac{puntos\ dentro\ n-esfera}{puntos\ totales}
Siendo el volumen del n_cubo: V_{n cubo} = 2^n ya que cada arista del cubo tendrá longitud 2: [-1,1].
El código para realizar el cálculo sería el siguiente:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#Calculo del volumen de una n-esfera. Como parametro recibe el numero de dimensiones y de puntos a generar, como salida devuelve el volumen
def calc_n_sphere_sobol(n_dimensions, n_points):
    #numero de puntos dentro de la n-esfera
    hits = 0
    #inicializador de la secuencia sobol
    sout = random.randint(1,10000)
    for i in range(0, n_points):
        #generar punto n-dimensional
        p, sout = sobol_lib.i4_sobol ( n_dimensions, sout )
        #transforma el numero del intervalo [0,1] al intervalo [-1,1]
        p = (p*2)-1
        j = 0
        s = 0
        #sumar cuadrado de las componentes
        while (jand s<=1):
            s+=p[j]**2
            j+=1
        #comprobar si el punto esta contenido en la n-esfera
        if (s<=1):
            hits += 1
    #devolver volumen de la n-esfera
    return  (2**n_dimensions)* float(hits) / n_points
Y los resultados comparados con el volumen dado por la fórmula:
Cálculo del espacio ocupado por una n-esfera de radio 1 para distintas dimensiones
Cálculo del espacio ocupado por una n-esfera de radio 1 para distintas dimensiones
Voi-lá, podemos ver cómo el método Monte Carlo se acerca a la función que describe el volumen de la n-esfera. Aquí están los datos obtenidos para 1.000.000 de puntos lanzados para las distintas dimensiones:
dimensionesvolumen n-cubovolumen n-esferaMonte CarloPuntos dentro
122,0000002,0000001000000
243,1415933,141700785425
384,1887904,188096523512
4164,9348024,936320308520
5325,2637895,263712164491
6645,1677135,17465680854
71284,7247664,73408036985
82564,0587124,06604815883
95123,2985093,2890886424
1010242,5501642,5169922458
1120481,8841041,837056897
1240961,3352631,314816321
1381920,9106290,950272116
El número de puntos que cumplen las condiciones de la n-esfera disminuye con el número de dimensiones, pero a la vez aumenta el volumen del n-cubo en el que está circunscrita.
 La relación máxima entre los puntos “acertados” y los lanzados respecto al volumen del n-cubo se da cuando llegamos a cinco dimensiones, y a partir de ahí el volumen de la n-esfera empezará a converger hacia cero.
Otra consideración a tener en cuenta es que cada vez que aumentamos una dimensión, el volumen del n-cubo aumenta. Como consecuencia, la muestra de 1.000.000 de puntos es menos significativa, por lo que los resultados empeoran si no aumentamos el número de puntos aleatorios utilizados.
Publicado por en