Últimamente estoy demasiado matemático. No se si será por los vídeos de derivando o por la cantidad de algoritmos de IA que me estoy metiendo en vena (ambos, segurísimo). Hace un tiempo en un curso online una de las cosas que aprendí como ejemplo para manejar distribuciones aleatorias en Python fue a crear una simulación con el método de Montecarlo para averiguar el valor de PI.

Este método consiste en «simular» dardos disparados aleatoriamente a una diana rodeada a su vez por un cuadrado cuyo diámetro  es igual al lado del cuadrado.

Como imagino que sabrás si  tienes una base mínima de probabilidad, la probabilidad que  tenemos de caer dentro del círculo (vamos, de acertar en la diana) es la probabilidad de dar en el círculo dividido entre la probabilidad de dar en cualquier parte del cuadrado (que incluye al círculo). Es decir: si cae en el círculo está también dentro del cuadrado. No obstante, si cae fuera del círculo está dentro del cuadrado pero fuera del círculo.

Probabilidad de dar en diana montecarlo piY cuando hablamos de probabilidad hablamos de la cantidad de puntos en el espacio que ocupa cada uno de ellos, es decir, sus áreas. Por tanto:

Sabemos que el área del círculo es PIxR^2, o también: PIxRxR, que es lo mismo. Y el área del cuadrado es su lado al cuadrado. No obstante, como sabemos que el lado del cuadrado es dos veces su diámetro: 2 x R x 2 x R.

Y si despejamos (que no es difícil en este caso)…

Vale. Pues matemáticamente habando, esto es todo amigos. ¿Qué pasa con esto ahora? Bien, pensemos una cosa: esto significa que si yo repito este experimento y cuento las veces que caigo en la diana y lo multiplico por 4, me va a dar pi. ¿Lo pillas? No es broma. Reléelo si hace falta. Pero es verdad. Bueno vamos a ver, no te va a dar pi exacto a la primera. Consiste en aproximarse al valor por simulaciones ojo, esto debe quedar claro.

La teoría ya la sabes. Ahora la práctica. Vamos a programarlo en Python. Te recomiendo Jupyter. Si lo conoces ya sabes por qué lo recomiendo. Y si no lo conoces, búscalo y pruébalo: entonces ya sabrás por qué lo recomendaba.

Debemos seguir el siguiente proceso general:

  • Generamos dos números aleatorios entre 0 y 1
  • Calculamos X X + Y Y
    • Si el valor es < 1 estaremos dentro del círculo
    • Si el valor es > 1 estaremos enconces fuera del círculo
    • Calculamos el número total de veces (puntos) que están dentro del círculo y dividimos entre el total de intentos.  Así obtenemos la probabilidad (aproximada) de caer en el cículo.
  • Con esa probabilidad podremos aproximar el valor de PI, y repetiemos el experimento un número de veces determinadas para obtener diferentes aproximaciones. Nunca va a ser exacto.

Si ya conoces un poco de esto, te habrás dado cuenta de que estamos haciendo una pequeña trampa: solo estamos usando el cuadrante superior derecho del círculo y del cuadrado (un 25 %). Si quisiéramos usar toda la superficie, deberíamos genera números entre -1 y 1. OJO: siempre suponiendo que nuestro origen de coordenadas está en el centro del círculo.

Es importante que utilicemos una distribución normal para generar los números aleatorios, ya que todos deben ser igual de probables, por eso usamos «numpy.random.uniform».

El código final es el siguiente, es bastante fácil de entender:

 

import matplotlib.pyplot as plt
import numpy as np

n = 1000
numExperimentos = 1000
mediaPi = 0
listaValoresPi = []
for i in range(numExperimentos):
    valor = 0
    X = np.random.uniform(0, 1, n).tolist()
    Y = np.random.uniform(0, 1, n).tolist()
    for j in range(n):
        Z = np.sqrt(X[j] * X[j] + Y[j] * Y[j])
        if Z <= 1:
            #Estamo dentro del círculo
            valor += 1
    valorFlotante = float(valor)
    valorPi = valorFlotante * 4 / n
    listaValoresPi.append(valorPi)
    mediaPi += valorPi

    
piFinal = mediaPi / numExperimentos
print(piFinal)
plt.plot(listaValoresPi)

Podeis descargaros el código fuente en Jupyter desde este enlace a mi repositorio, donde iré guardando los demás ejemplos que vaya publicando.