Calculando cosas mientras jugamos a la ruleta

Hoy vamos a tratar una serie de simulaciones que deben su nombre al lujoso casino de Monte Carlo y a la capacidad de la ruleta de comportarse como un generador de números aleatorios. Vamos a ver varios ejemplos que ilustran el cálculo aproximado de áreas y del número $latex pi$ mediante una simulación Monte Carlo.
En primer lugar vamos a calcular un área considerando el siguiente problema: Tenemos un recuadro que mide 10 km de lado y dentro de este recuadro tenemos un polígono de lados irregulares que, por la razón que sea, no sabemos cuanto ocupa, es decir, cuál es su área.
Para resolver este problema de forma aproximada vamos a usar shapely, además de los omnipresentes numpy y matplotlib. Empezamos importando todo lo que necesitamos:
[sourcecode language=”python”]
from shapely.geometry import Polygon, Point
import numpy as np
import matplotlib.pyplot as plt
[/sourcecode]
Generamos los dos polígonos, el cuadrado que sí sabemos lo que mide (poligono_externo) y el polígono con área desconocida (polígono interno):
[sourcecode language=”python”]
poligono_externo = Polygon(((0,0),(0,10),(10,10),(10,0),(0,0)))
poligono_interno = Polygon(((3,2),(5,7),(9,6),(6,5),(7,1),(3,2)))
[/sourcecode]
Si dibujamos ambos polígonos tendremos algo parecido a lo siguiente:
[sourcecode language=”python”]
plt.ion()
plt.xlim(-1,11)
plt.ylim(-1,11)
plt.plot(poligono_externo.boundary.xy[0], poligono_externo.boundary.xy[1])
plt.plot(poligono_interno.boundary.xy[0], poligono_interno.boundary.xy[1])
[/sourcecode]

Ahora creamos una función que será la que se encarga de hacer el cálculo:
[sourcecode language=”python”]
def simulacion_mc(num_pruebas):
xvect = np.random.rand(num_pruebas)*10.
yvect = np.random.rand(num_pruebas)*10.
cont = 0
for x, y in zip(xvect, yvect):
if poligono_interno.contains(Point(x, y)):
cont += 1
return poligono_externo.area * cont / num_pruebas
[/sourcecode]
¿Qué hemos hecho en la anterior pieza de código? Hemos generado unos valores aleatorios para la variable x y para la variable y de la posición. Comprobamos si cada uno de los puntos generados se encuentran dentro del polígono de área desconocida y los contamos. Finalmente devolvemos el cálculo del área aproximada. El área se calcula aprovechando que la probabilidad de que un punto esté dentro del área del polígono interior será proporcional al área de ese polígono:
$latex frac{Area_{text{poligono exterior}}}{Area_{text{poligono interior}}}proptofrac{text{puntos contenidos en el poligono exterior}}{text{puntos contenidos en el poligono interior}} = frac{text{num pruebas}}{text{cont}} $
Por tanto, si despejamos el Área del polígono interior de la anterior fórmula tenemos que el área aproximada la podemos obtener de la siguiente forma:
$latex Area_{text{poligono interior}} = Area_{text{poligono exterior}}frac{text{cont}}{text{num pruebas}} $
Donde num_pruebas es el número de puntos (posiciones) que usamos en la simulación (todos contenidos dentro del polígono exterior) y cont es el número de puntos que han ‘caído’ en el polígono interior.
Una simulación de 1000 puntos lo que hace es lo que podemos ver en la siguiente gráfica:
[sourcecode language=”python”]
num_pruebas = 1000
xvect = np.random.rand(num_pruebas)*10.
yvect = np.random.rand(num_pruebas)*10.
plt.xlim(-1,11)
plt.ylim(-1,11)
plt.plot(poligono_externo.boundary.xy[0], poligono_externo.boundary.xy[1])
plt.plot(poligono_interno.boundary.xy[0], poligono_interno.boundary.xy[1])
plt.plot(xvect, yvect, ‘r.’)
[/sourcecode]

El número de puntos dentro del polígono interno sería cont en la función que hemos definido anteriormente mientras que el número de puntos dentro del polígono externo serían todos los puntos usados, es decir, num_pruebas.
Vamos a hacer varias simulaciones con diferentes números de puntos y representar el valor del área obtenida en función del número de puntos usado:
[sourcecode language=”python”]
## Aviso, este cálculo puede tardar bastante…
resultados = [simulacion_mc(i) for i in np.arange(0, 10000, 10) + 10]
plt.plot(np.arange(0,10000,10) + 10, resultados)
plt.xlabel(u’número de puntos aleatorios usados’)
plt.ylabel(u’Área del polígono interior’)
plt.show()
[/sourcecode]
El resultado sería el siguiente:

Vemos que los valores obtenidos están en torno a 15 y pico y que hay menos dispersión a medida que aumentamos el número de puntos usados en la simulación. Si ahora usamos la propiedad area de la clase Polygon podemos ver el área que tiene el polígono interno:
[sourcecode language=”python”]
print poligono_interno.area
[/sourcecode]
El resultado, si hemos hecho todo bien, sería 15.5 km². Parecido a los resultados de nuestras simulaciones!!!
Todo el código junto (sin la generación de gráficas) quedaría de la siguiente forma:
[sourcecode language=”python”]
from shapely.geometry import Polygon, Point
import numpy as np
import matplotlib.pyplot as plt
## Definimos dos polígonos, uno del que conoceremos su área
## y un segundo de área desconocida (a priori)
poligono_externo = Polygon(((0,0),(0,10),(10,10),(10,0),(0,0)))
poligono_interno = Polygon(((3,2),(5,7),(9,6),(6,5),(7,1),(3,2)))
## Función que calcula nuestra área aproximada
def simulacion_mc(num_pruebas):
xvect = np.random.rand(num_pruebas)*10.
yvect = np.random.rand(num_pruebas)*10.
cont = 0
for x, y in zip(xvect, yvect):
if poligono_interno.contains(Point(x, y)):
cont += 1
return poligono_externo.area * cont / num_pruebas
## Lanzamos varias simulaciones
## Aviso, este cálculo puede tardar bastante…
resultados = [simulacion_mc(i) for i in np.arange(0, 10000, 10) + 10]
[/sourcecode]
Ahora dejo un código para calcular $latex pi$ usando 10000 puntos aleatorios. Si alguien no lo entiende que pregunte en los comentarios pero se trata de la misma metodología usada anteriormente (pistas: 1, 2):
[sourcecode language=”python”]
import numpy as np
import matplotlib.pyplot as plt
## Función que calcula nuestra área aproximada
def simulacion_mc(num_pruebas):
xvect = np.random.rand(num_pruebas)
yvect = np.random.rand(num_pruebas)
cont = 0
for x, y in zip(xvect, yvect):
if np.sqrt(x * x + y * y) <= 1:
cont += 1
return 4. * cont / num_pruebas
print simulacion_mc(10000)
[/sourcecode]
Espero que a alguien le resulte útil.
P.D.: Escribiendo las fórmulas en latex de más arriba no encontraba como escribir el símbolo $latex propto$ y lo he encontrado gracias a ésta web (algún día hablaremos sobre metodologías para hacer lo que se hace en esa web :-))

8 pensamientos sobre “Calculando cosas mientras jugamos a la ruleta”

  1. Muy, pero muy interesante el artículo, aunque no entiendo bien como se podría aplicar a un sistema real, en cambio veo que las herramientas que usas y lo fácil que se implementa está llamativo. Para ser mas específico, la librería de geometría me ha gustado mucho.
    Sobre el cálculo de pí, no me termina de convencer lo de obtener a partir del área, ya que la suma perimetral mediante póligonos infinitesimales me parece mucho más eficiente.
    Sin duda, es un gran artículo ¡felicidades!
    ha.. sobre latex… yo generalmente uso TeXstudio y éste tiene entre sus símbolos a varpropto y propto, me imagino que ese símbolo es el que usas, te recomiendo.

  2. Muy buen artículo!
    Me pareció interesante lo fácil que implementás, la librería de geometría está muy buena.
    Sobre el cálculo de pi, no me termina de convencer el cálculo mediante puntos aleatorios y el 4. como factor de escala, ¿es por los 4 cuadrantes verdad?.
    Prefiero una aproximación de pi mediante suma perímetral de polígonos infinitesimales
    De todas formas el artículo está muy intersante aunque no encuentro (en este momento) donde se podría aplicar, tal vez en un análisis por imágenes y comparaciones aunque esto sería mucho mas trabajoso (muchos cálculos) comparado con análisis de señales digitales.
    Nota1: ha…por cierto está conviertiendo un símbolo < en "<" en el último código 😀
    Nota 2: sobre latex, te recomiendo TeXstudio, que dentro sus símbolos incluye ese que necesitabas! 😉

    1. Hola.
      El 4 es por los 4 cuadrantes, sí. Solo uso uno de los 4 cuadrantes por lo que el área de circunferencia usada solo es un cuarto y por tanto es pi * r^2 / 4.
      Un ejemplo práctico de uso. En mi trabajo usamos máquinas que tienen unas tasas de fallos que no conocemos. Pero para los planes financieros necesitamos estimar el coste que suponen esos fallos para saber si el negocio tendrá un retorno u otro y si sería rentable o no. Esos fallos no suponen el mismo coste si suceden, por ejemplo, a las 02.00h o a las 15.00h. Mediante simulaciones Monte Carlo podemos hacer una estimación de lo que podría ser una situación ‘normal’ de esos fallos.
      Nota 1: Tengo que corregir lo del <, que wordpress me la ha jugado ahí. Gracias.
      Nota 2: No podemos usar TeXstudio ni mathjax dentro de wordpress.com, por eso ponemos las fórmulas como las ponemos.

  3. Pingback: Práctica de Montecarlo: La panadería de Pierre | Jobliz's Log

  4. Pingback: Números aleatorios en Python con NumPy y SciPy « Pybonacci

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

seventy − = sixty three