número PI (método Monte Carlo)

PARTE 1. EXPERIMENTO

Muchos, incluyendo a mí, han despreciado estadística todo este tiempo hasta que ha llegado el tema de IA que en mayoría de casos utiliza herramientas de estadística. Esto ha traído un cambio radical al mundo. Antes, para calcular los recorridos y trayectos más cortos deberíamos inventar algoritmos muy complejos y poco claros tipo SIMPLEX u otros. Ahora google maps utiliza datos que provienen de los propios conductores para hacer comparaciones e inferencias tan precisas que desaparece necesidad de aplicar algoritmos de tipo determinista que poco se ajustan a problemas de vida real.

Un ejemplo simple. Para calcular el número PI lo primero que se ocurre es coger una cuerda, hacer con ella una circunferencia, en el suelo, medir el radio y luego dividir longitud de la cuerda entre longitud del radio. Sabiendo que longitud de la circunferencia es 2PIr al dividir entre r obtenemos el doble del número PI. No nos va a dar valor exacto porque no podemos crear figuras perfectas de esta forma.

Para aproximarnos mucho mejor podemos usar algoritmo de Pitágoras pero hoy vamos hacerlo de otra forma.

Sabemos que área de una circunferencia es PIr^2. Vamos a dibujar una circunferencia inscrita en un cuadrado.

Área del cuadrado es (2r)^2

Entonces área de circunferencia es PIr^2/(2r)^2 = PI/4 veces menos que el cuadrado. Si hacemos un experimento, tirando con ojos cerrados un garbanzo en este campo, la probabilidad que el garbanzo cae sobre la circunferencia es de PI/4

Solo nos falta repetir el experimento varias veces. Bastantes veces. En el caso si el garbanzo cae dentro de las circunferencia lo vamos a llamar favorable (F) y el número total de estos casos lo vamos a llamar (f), el número total de intentos lo vamos a denotar por n. Esta proporción f/n debe aproximarse a PI/4. Así que el resultado final lo vamos a multiplicar por 4 y obtenemos una aproximación a PI. Por ley de números grandes a mayor número de intentos la aproximación sería mejor. Vamos a diseñar los experimentos con Python.

Creamos un campo 1x1 la verdad que tamaño no importa ya que se ha visto que el radio se “simplifica” en el desarrollo.

Observa que si trabajamos con una cuarta parte de imagen la proporción sería la misma. Pero esto ayuda a evitar situar el origen de ejes de coordenadas en el centro de la imagen y crear semirrectas negativas.  

No vamos a crear listas de puntos favorables. Mejor trabajamos con el radio que se calcula usando teorema de Pitágoras r = sqrt(x^2+y^2). Si lo único que necesitamos es conocer si el garbanzo queda en la circunferencia, entonces solo nos interesa un parámetro. El radio máximo en este caso es 1. Por esto podemos decidir si el garbanzo cae en zona favorable comprobando que r < 1

Repetimos experimento varias veces y obtenemos resultados. En cada iteración restamos del valor real el valor de la fracción. Éste sería el error de aproximación. Veamos en la gráfica que a mayor número de iteraciones este error disminuye. 

Éste método tiene muchas ventajas ya que nos proporciona protocolo para cálculo de áreas complejas sin tener que inventar algoritmos complicados. Se apoya en el teorema central del límite y en otros resultados matemáticos rigurosos por lo cual es un método fiable.

La desventaja de este método es el resultado inexacto y no solo por el tema de “fallos” probabilísticos. Hemos usado en el desarrollo una figura o bien dibujada a mano para tirar garbanzo encima de ella, o bien una circunferencia construida con herramientas informáticas. En ambos casos la figura no es perfecta. De hecho el ordenador no opera con los números reales. Son infinitos y el ordenador tiene memoria limitada. Por lo cual no puede trazar una circunferencia perfecta.

Pero ¿Qué más da si el método es aplicable a la vida real?

 PARTE 2. CONTROL DE ERROR

Ahora vamos a formalizar nuestro estudio y dotarle de un intervalo de confianza. Para asegurarnos el valor de PI con un porcentaje alto de probabilidad.

Volvemos a la interpretación gráfica para diseñar el experimento como proporción. Sabemos que la proporción de éxito es de PI/4 ¿Qué tamaño de muestra n debe tomarse para que el intervalo de confianza de 99% sería menor de 0,001?

Se sabe que la de proporción tiene distribución normal de media p = PI/4 y desviación típica sqrt(p*q/n)

p = PI/4
q = 1 – PI/4

si quiero un intervalo de eficacia 99% debo recurrir a la tabla normal. 

Aciertos en 99% da lugar a fallos 1% y la mitad es 0,5% que es 0,005 en términos de la probabilidad. Si vamos a la tabla normal con datos mayor de 0.5 tenemos que optar por el “complementario” 1 - 0,005 = 0,995

Aquí está:

z = 2.57

IC99 = [p-z·sqrt(p*q/n); p+z·sqrt(p*q/n)]

Si quiero que longitud total del intervalo no supera 0,001 entonces su amplitud no debe superar 0,0005 entonces z·sqrt(p*q/n) < 0,0005. De aquí solo tenemos que despejar la n

Sqrt(p*q/n)< 0,0005/z
p*q/n<(0,0005/z)^2
p*q<((0,0005/z)^2)*n
n> (p*q)/(0,0005/z)^2 = 4452967.790506449

Resulta que tenemos que tomar muestra mayor de 4452968 para asegurar intervalo tan pequeño a 99%
Probamos:

nos da como resultado pi = 3.1421353128969263

¿Está dentro de intervalo?

PI - pi  =  -0.0005426593071331531 No está en el intervalo. Pero casi.

El código:

import random
import math
import matplotlib.pyplot as pltÇ
n = 4452968; # cantidad de experimentos que vamos a usar
N= [];
f = 0;
L = [];
i = 1;
while i <= n:
    x = random.random();
    y = random.random();
    if pow(x,2)+pow(y,2) <= 1:
        f = f +1;
    L.append((f/i)*4-math.pi)
    i = i + 1;
    N.append(i)
# print(L)
   
 
fig, ax = plt.subplots()
plt.plot([0, n], [0, 0], color='green')
ax.plot(N, L)
plt.show()
print((f/n)*4)