¿Cómo calcular el valor de Pi usando programación?
En esta clase, ejecutaremos un algoritmo para calcular el estimado del número Pi mediante simulaciones y técnicas estadísticas. Es fascinante descubrir cómo mediante un simple uso de geometría y probabilidad, podemos determinar un valor aproximado de este número irracional tan crucial en matemáticas.
¿Cómo crear una función que lance agujas?
Para empezar, se requiere una función en Python que simule el lanzamiento de agujas sobre un plano para calcular Pi. Aquí están los pasos básicos:
Definir la función aventar_agujas: Esta función recibe el número de agujas como parámetro.
defaventar_agujas(numero_agujas): adentro_del_circulo =0for _ inrange(numero_agujas): x = random.random()* random.choice([-1,1]) y = random.random()* random.choice([-1,1]) distancia = math.sqrt(x **2+ y **2)if distancia <=1: adentro_del_circulo +=1return4* adentro_del_circulo / numero_agujas
Utilizar funciones aleatorias: Con la biblioteca random se generan coordenadas x y y en el rango de [-1,1].
Calcular la distancia: Usar el Teorema de Pitágoras para determinar si la aguja cae dentro del círculo de radio 1.
¿Cómo utilizar la estadística para mejorar la estimación?
Usando conceptos estadísticos como la media y la desviación estándar, mejoramos la precisión de nuestra estimación.
Importar funciones estadísticas: Utilizar una librería previamente creada con funciones de media y desviación estándar.
from estadisticas import media, desviacion_estandar
Crear una función de estimación: Comienza ejecutando simulaciones múltiples veces para obtener un valor promedio de estimación.
Ejecutar la función principal: Llamar esta función con valores como precision=0.01 y numero_intentos=1000 para obtener un resultado.
Esta clase no solo explica un método eficiente para estimar Pi, sino también ofrece un enfoque práctico del uso de estadísticas en simulaciones. A través de la práctica, el estudiante irá familiarizando el uso de aleatoriedad y estadística en programación. ¡Esperamos que hayas disfrutado y aprendido cómo estas técnicas pueden aplicarse para resolver problemas reales! No olvides compartir tus dudas y comentarios en la sección correspondiente.
Solución del ejercicio de estimación de pi con apoyo grafico para entender que esta sucediendo durante la simulación, en único dato a considerar es la precisión ya que entre más cifras significativas se usen el tiempo de computo se hace exponencialmente largo.
import os
import random as random
import math as math
import statistics as std
from bokeh.plottingimport figure, output_file, show
os.system('cls')tries =10puntos =100def estimar_pi(puntos): in_circle_x =[] in_circle_y =[] out_circle_x =[] out_circle_y =[] pi_array =[]for i inrange(puntos): pos_x = random.uniform(-1,1) pos_y = random.uniform(-1,1)if math.sqrt(pos_x**2+ pos_y**2)<=1: in_circle_x.append(pos_x) in_circle_y.append(pos_y)else: out_circle_x.append(pos_x) out_circle_y.append(pos_y) estimate_pi =(4*len(in_circle_x))/ puntos
return(estimate_pi, in_circle_x, in_circle_y, out_circle_x, out_circle_y) #return()def crear_muestra(tries): pi_array =[]for i inrange(tries): pivalor,in_circle_x, in_circle_y, out_circle_x, out_circle_y =estimar_pi(puntos) pi_array.append(pivalor)return(pi_array, in_circle_x, in_circle_y, out_circle_x, out_circle_y)deviation =1 # Valor inicial para poder empezar el ciclo whilepresicion =0.1 # Precision en cifras significativas para determinar el
# valor estimado de pi(CUIDADO!!!!)iteration =1while deviation >=(presicion /1.96): pi_array, in_circle_x, in_circle_y, out_circle_x, out_circle_y =crear_muestra(tries) deviation = std.stdev(pi_array) variance = std.variance(pi_array) mean = std.mean(pi_array)print(f'------------------ Iteracion number: {(iteration)} ------------------')print(f'Standard deviation: {round(deviation,5)}, Variance : {round(variance,5)}, pi estimated: {round(mean,5)}')print(f'Numero de intentos {tries}, Numero de puntos {puntos}\n\n')#print(f) puntos *=10 tries *=10 iteration +=1output_file("line.html")plot =figure(plot_width=600, plot_height=600)plot.circle(in_circle_x, in_circle_y, size=5, color="red", alpha=0.5)plot.circle(out_circle_x, out_circle_y, size=5, color="navy", alpha=0.5)show(plot)
Espectacular Carlos
Gracias!!
Goolge, facebook y amazon preguntan este problema de calcular pi en sus entrevistas de trabajo, aquí dejo un vídeo que lo explica.
https://www.youtube.com/watch?v=pvimAM_SLic
Muy buen aporte!, lo explica gráficamente a la perfección
Excelente aporte, gracias Juan
Un 95% de confianza significa que el 95% de mis valores estarán dentro del rango que yo voy a definir con la precisión, en este caso el rango es pi +- 0,01. En la tabla de la distribución normal obtenemos que ese intervalo se encuentra en el rango de media +- 1,96sigma, por lo que 1,96 sigma debe ser igual a la precisión. Lo que quizás hace que el código sea confuso es que se le asigna a sigma el valor de precisión, lo cual lo único que permite es entrar al ciclo while, ya que luego es inmediatamente reasignado, por lo que se le podría asignar cualquier número mayor a precisión / 1,96.
Excelente explicación. Para complementar algo pequeño: por lo tanto el programa debe hacer tantas iteraciones sean necesarias (es decir arrojar y arrojarr agujas) para que sigma se haga tan pequeño que al multiplicarlo por 1.96 me de igual a la precisión que estoy buscando.
Gracias por la explicación, me ayuda a reforzar aún más.
Adjunto mi solución y las gráficas de las simulaciones por si a alguien le sirve:
Para el que no pudo entender del todo el ultimo while
Explicacion:
Mientras que nuestra varianza sea mayor a nuestra confiaza en el calculo, se ejecutara el while
en el ejemplo, la confianza es de 95% porque da 0.0196, osea que nuestra aproximacion de PI va a estar 0.0196 alejado del PI real
while sigma >= presicion /1.96:```
muchas gracias, me ha servido tu aporte
Tuve que repasar un poco pitágoras, acá les dejo una pequeña demostración visual de el cálculo de la distancia desde el centro:
hice algo parecido, gracias por compartir
Comparto una implementación diferente y un poco mas sencilla del calculo de pi :
# una forma mucho mas corta de implementarlo
def pi(n, batch=1000): t =0for i inrange(n // batch): p = np.random.rand(batch,2) p =(p * p).sum(axis=1) t +=(p <=1).sum()return4* t / n
Una explicación mas detallada la pueden ver aquí en mi GitHub.
Claro, es mucho más bonito usar numpy y un array de numeros aleatorios
Me costo un poco comprender lo de la precisión. Pero según entiendo esta precisión corresponde simplemente a definir cual va a ser mi 2*sigma y así se logre ocupar la cantidad de agujas necesarias para cumplir con ese nivel de precisión. Creo que eso falto explicar un poco mejor a demás que en el vídeo se habla de la precisión (0,01) como un porcentaje y creo que eso es incorrecto ya que es adimensional al representar simplemente un margen del numero pi.
para complementar habría que decir que la idea entonces es que el 95% de los resultados se encuentre entre la media mas/menos 0,01 ( la precisión ) . No lo recuerdo bien pero para esto creo que se requieren de mas de 100000 agujas.
Estadísticamente probable != Válido
Hola ⚡
Hice esta animación para poder visualizar como va cambiando el valor de pi y como van cayendo los puntos de manera aleatoria en el cuadro 👀, se ve bien cool!
!Simulación
Si quieren echarle un ojo al código, esta en este enlace
Me perdí un poco en el programa en la línea
while sigma >= precision / 1.96
minuto 9:40
la explicación no la entendí ya que en la regla empírica estabamos usando 2 desviaciones estandar para el 95% y no comprendo el cambio a 1.96 ¿Por qué ese número y no 1.98 o 1.97?
la desviavicion estandar multiplicada por 2 nos da 95% y un poquito mas. para que nos devuelva 95% exacto, lo esta multiplicando por 1.96
Aquí se menciona una explicación un tanto más técnica y elaborada de por qué el profesor utilizar el 1.96.
Hola,
Después de casi 10 minutos corriendo el código, obtuve los siguientes resultados:
Est = 3.14233, sigma = 0.05133, agujas = 1000
Est = 3.14048, sigma = 0.03577, agujas = 2000
Est = 3.14154, sigma = 0.02563, agujas = 4000
Est = 3.14119, sigma = 0.01807, agujas = 8000
Est = 3.1413, sigma = 0.01324, agujas = 16000
Est = 3.1421, sigma = 0.00923, agujas = 32000
Est = 3.1412, sigma = 0.00648, agujas = 64000
Est = 3.14151, sigma = 0.00473, agujas = 128000
Hola comunidad, ya tengo un mini proyecto en mente para calcular la probabilidad de digamos "Sacar 3 cartas +4 en tu mano en un juego de Uno"; y así varios ejercicios.
.
Pero saben de algún sitio donde pueda encontrar ejercicio para practicar Simulaciones de Montecarlo? Si tienen ideas o propuestas también son válidas :) y podemos mirarlas juntos
a mi también me interesa un sitio para poder practicar lo aprendido.
tengo este repositorio de github donde la idea es subir varias simulaciones con el tiempo, puedes hacerme un pull request si gustas
https://github.com/Sgewux/Simulations
Según fórmula:
A_cuadrado = 4R
A_circulo = piR
Igualando ambas ecuaciones (despejando R)
A_cuadrado/4 = A_circulo/pi = R
Despejando pi:
pi = 4*A_circulo/A_cuadrado
Puesto que las agujas adentro del circulo son proporcionales al area del circulo y el numero de agujas es proporcional al area de cuadrado pordemo reemplazarlos
Recuerda que trabajamos con un cuadrado de 2 X 2 de area 4. Lo que implica un circulo de radio 1.
Entendí todo el proceso del programa, y entendi la parte de estimar pi, pero lo que no entiendo exactamente es:
De donde sale la formula presicion/1.96 ?
Se que 1.96 es la cantidad de "desviacion estandar" pero no logro entender de donde salió o cual es su representación grafica
Nosotros queremos obtener el valor de pi, pero hay que recordar que este número es irracional, es decir que tiene infinitos decimales, y calcular todos es imposible. Así que definimos una precisión que sería el criterio de parada de nuestro algoritmo, pero queremos que el 95% de los datos tengan una precisión menor a 0.01, es decir que los datos se encuentren entre (media - 1.96 * sigma y media + 1.96 * sigma) o se puede ver también así (media -0.01 y media +0.01). por lo tanto {0.01 = 1.96 * sigma}.
Tal vez con esta imagen te quede más claro.
La región sombreada serian el 95% de los datos y para asegurarnos de eso aumentamos el número de agujas hasta llegar a ese criterio de parada.
Muchas gracias bro
Si, había entendido la formula pero nunca el porque, y pues eso es lo mas importante, entender de donde fue el origen
Me quedo muy claro con tu explicación, gracias por el tiempo
Saludos
https://github.com/karlbehrensg/programacion-dinamica-y-estocasticaCálculo de PI
Calcularemos PI con puntos al azar esparcidos en un plano cartesiano utilizando los scripts de desviación estándar y media que creados anteriormente. Queremos tener un 95% de certeza, entonces para ello realizaremos el cálculo para 1/2 del área de un circulo, optimizando nuestros recursos.
import random
import math
from estadisticas import desviacion_estandar, media
defaventar_agujas(numero_de_agujas): adentro_del_circulo =0for _ inrange(numero_de_agujas): x = random.random()* random.choice([-1,1]) y = random.random()* random.choice([-1,1]) distancia_desde_el_centro = math.sqrt(x**2+ y**2)if distancia_desde_el_centro <=1: adentro_del_circulo +=1# La variable adentro_del_circulo representa 1/4 del área del círculo,# y como solo utilizaremos 1/2 vamos a multiplicarlo por 2.return(2* adentro_del_circulo)/ numero_de_agujas
defestimacion(numero_de_agujas, numero_de_intentos): estimados =[]for _ inrange(numero_de_intentos): estimacion_pi = aventar_agujas(numero_de_agujas) estimados.append(estimacion_pi) media_estimados = media(estimados) sigma = desviacion_estandar(estimados)# La variable media_estimados tiene los resultados sobre 1/2 del área del# círculo. Para obtener la estimación de PI completo lo vamos a multiplicar por 2.print(f'Est={round(media_estimados,5)*2}, sigma={round(sigma,5)}, agujas={numero_de_agujas}')return(media_estimados, sigma)defestimar_pi(precision, numero_de_intentos): numero_de_agujas =1000 sigma = precision
while sigma >= precision /1.96: media, sigma = estimacion(numero_de_agujas, numero_de_intentos) numero_de_agujas *=2return media
if __name__ =='__main__': estimar_pi(0.01,1000)