Implementación del Cálculo de PI

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.plotting import figure, output_file, show
os.system('cls')

tries = 10
puntos = 100

def estimar_pi(puntos):
    in_circle_x = []
    in_circle_y = []
    out_circle_x = []
    out_circle_y = []
    pi_array =[]

    for i in range(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 in range(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 while
presicion = 0.1     # Precision en cifras significativas para determinar el 
                    # valor estimado de pi (CUIDADO!!!!)

iteration = 1

while 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 += 1

output_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)

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

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.

Adjunto mi solución y las gráficas de las simulaciones por si a alguien le sirve:

# -*- coding: utf-8 -*-
import random as rd
import matplotlib.pyplot as plt

def pi_simulated(points):
    square_x, square_y = [], []
    circle_x, circle_y = [], []
    circle_points = 0

    for _ in range(points):
        square_x.append(1 - 2 * rd.random())
        square_y.append(1 - 2 * rd.random())
    square_coord = list(zip(square_x, square_y))

    # are the points inside the circle?
    for i in range(points):
        if ((square_coord[i][0]) ** 2 + (square_coord[i][1]) ** 2) < 1:
            circle_x.append(square_coord[i][0])
            circle_y.append(square_coord[i][1])
            circle_points += 1
        else:
            circle_x.append(0)
            circle_y.append(0)
            circle_points += 0
    circle_coord = list(zip(circle_x, circle_y))
    computed_pi = round((circle_points / points) * 4, 4)

    # plot
    fig, ax = plt.subplots(figsize=(5, 5))
    series1 = ax.scatter(square_x, square_y, color='b', s=5)
    series2 = ax.scatter(circle_x, circle_y, color='r', s=5)
    plt.ylim(-1, 1)
    plt.xlim(-1, 1)
    ax.set(xlabel='X', ylabel='Y',
           title=f'Pi Number with {points} points. Computed Pi: {computed_pi}')
    ax.legend([series1, series2], ('Square', 'Circle'), loc='upper right', facecolor='white', framealpha = 1)
    ax.grid()
    plt.show()

    return computed_pi


def main(points):
    analytical_pi = 3.1416
    error = analytical_pi
    while error > 0.001:
        pi_sim = pi_simulated(points)
        error = abs(pi_sim - analytical_pi)
        print(pi_sim, error)
        points *= 2

if __name__ == '__main__':
    points = int(input('Numero de puntos'))
    #print(pi_simulated(points))
    main(points)

Yo pensando que era tirar las agujas en horizontal y en realidad era como tirar dardos jajajajaja

Esta linea de codigo …

while sigma >= precision / 1.96:

Lo deduje de la formula que entrego el profesor en la clase de la regla empírica:

Alguien llego a la misma conclusión o la explicación es otra?

Les dejo una explicación para quienes no entendieron lo del intervalo de confianza

https://www.youtube.com/watch?v=2wugQGs1GNY

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:```

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:

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 = 0
    for i in range(n // batch):
        p = np.random.rand(batch, 2)
        p = (p * p).sum(axis=1)
        t += (p <= 1).sum()
    return 4 * t / n

Una explicación mas detallada la pueden ver aquí en mi GitHub.

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.

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!

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?

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

https://github.com/karlbehrensg/programacion-dinamica-y-estocastica
Cá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

def aventar_agujas(numero_de_agujas):
    adentro_del_circulo = 0

    for _ in range(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


def estimacion(numero_de_agujas, numero_de_intentos):
    estimados = []
    for _ in range(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)

def estimar_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 *= 2

    return media

if __name__ == '__main__':
    estimar_pi(0.01, 1000)

Vamos a la consola y ejecutamos nuestro programa.

python3 calculo_pi.py   # Ejecutamos nuestro script.

# Y estos serán nuestros resultados.
Est=3.14234, sigma=0.02594, agujas=1000
Est=3.13966, sigma=0.01795, agujas=2000
Est=3.143, sigma=0.01272, agujas=4000
Est=3.14076, sigma=0.00949, agujas=8000
Est=3.14142, sigma=0.00678, agujas=16000
Est=3.14154, sigma=0.00457, agujas=32000

De hecho en el trading existe un indicador que se denomina metodos estocasticos la cual nos indica el punto donde se presenta una sobre compra o una sobreventa de un activo financiero

si buscamos en la tabla de la distribución normal Z(0.95) osea un nivel de confianza del 95%, encontramos que Z(0.95)=1.95. Ahora hay que entender que la normal la media es cero y la desviación es uno.
Lo que hacemos es normalizar nuestra precisión y eso lo hacemos dividiendo por 1.96

Como ya vieron, el valor de 1.96 proviene de la tabla Z. Dado que necesitamos un 95% de confianza, se busca en la tabla el valor más próximo a 0.95. Para quienes no conocemos y sabemos bien sobre esto, encontré una forma de obtener el valor (1.96) a partir del porcentaje de confianza deseado:
Si necesitamos 95 de grado de confianza, devuelve 1.96
Si necesitamos 95 de grado de confianza, devuelve 1.645

from scipy.stats import norm

confianza = 0.95
z = round(norm.interval(0.95, loc=0, scale=1)[1], 3)
print(z)

Mi compu literalmente murio a las 64mil estimaciones con un total de media de 3.1416
Nice algoritmo

Sacar los números aleatorios de -1 a 1 también se puede con:
x = random.random()*2-1

Una explicación más sintetizada.

https://www.youtube.com/watch?v=pvimAM_SLic

import math

suma = 0
for valor in range(1, 100000000):
    suma = suma + (1 / valor ** 2)

pi = math.sqrt(suma * 6)
print(f'El valor de pi es : {pi}')

Una solución determinista, propuesta por Euler para resolver el problema de Basilea.

Este es un pequeño Script que yo hice. No tiene gráficas, pero calcula con varias simulaciones el valor.

import random

def aventar_agujas(numero):
    distancias = []
    for _ in range(numero):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        distancia = (x**2 + y**2)**(1/2) #distancia del centro
        #print(f'X : {x}, Y: {y},  dist = {distancia}')
        distancias.append(distancia)
    return distancias

def calcular_areas(distancias):
    agujas_dentro = 0
    agujas_fuera = 0
    for distancia in distancias:
        if distancia < 1:
            agujas_dentro += 1
        elif distancia >= 1:
            agujas_fuera += 1
    area_circ = (4 * agujas_dentro) / (agujas_fuera+agujas_dentro)
    return area_circ

def main(num_simulaciones, agujas):
    valores_pi = []
    for simulacion in range(num_simulaciones):
        distancias = aventar_agujas(agujas)
        area_circulo = calcular_areas(distancias)
        valores_pi.append(area_circulo)
        print(f'PI en simulacion {simulacion + 1} = {area_circulo}')
    promedio = sum(valores_pi) / len(valores_pi)
    print(f'Valor promedio de pi = {promedio}')

if __name__ == '__main__':
    num_simulaciones = int(input('Cuantas veces quieres correr la sim? : '))
    main(num_simulaciones, 10000000)

En la consola se ve algo así:

Cuantas veces quieres correr la sim? : 5
PI en simulacion 1 = 3.141802
PI en simulacion 2 = 3.1426052
PI en simulacion 3 = 3.1424452
PI en simulacion 4 = 3.14152
PI en simulacion 5 = 3.1411352
Valor promedio de pi = 3.1419015200000002

Este tipo de aproximaciones estocásticas a problemas de la vida real se pueden tener cuando no puedo medir de algún modo ciertos parámetros de un sistema, más si puedo apoyarme en los efectos que estos generan sobre otros sistemas. Por ejemplo simular el comportamiento de un robot en un ambiente que jamás he medido físicamente, sin embargo del cual tengo alguna información en fotografías que me permite introducirla como insumos para las labores o acciones de dicho robot.

Notas:

• random.random() → genera un numero aleatorio entre 0 y 1
• Una solución estadísticamente válida no es una solución válida, la solución estadísticamente valida depende de la desviación estándar (la cantidad de muestras que tomamos), mientras que una solución valida dependerá del planteamiento del problema (algoritmo)
• La precision dependerá de el valor sigma: p=constante_confiabilidad*sigma y de la confiabilidad que queramos obtener (1sigma = 68%, 1.96sigma=95%, 2sigma=95.45%…). Entonces para obtener una precision de 0.01 necesitamos un sigma de 0.01/1.96=0.0051

Codigo:

import random
import math
from estadisticas import desviacion_estandar, media

def aventar_agujas(numero_agujas):
    adentro_del_circulo = 0

    for _ in range(numero_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

    return(4 * adentro_del_circulo) / numero_agujas

def estimacion(numero_agujas, numero_de_intentos):
    estimados = []
    for _ in range(numero_de_intentos):
        estimacion_pi = aventar_agujas(numero_agujas)
        estimados.append(estimacion_pi)

    media_estimados = media(estimados)
    sigma = desviacion_estandar(estimados)
    print(f'Est={round(media_estimados, 5)}, sigma={round(sigma, 5)}, agujas={numero_agujas}')

    return (media_estimados, sigma)

def estimar_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 *= 2

    return media


if __name__ == '__main__':
    estimar_pi(0.01, 1000)```

import random
import math
from estadisticas import desviacion_estandar, media

def aventar_agujas(numero_de_agujas):
	adentro_del_circulo = 0

	for _ in range(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

	return (4 * adentro_del_circulo) / numero_de_agujas

def estimacion(numero_de_agujas, numero_de_intentos):
	estimados = []
	for _ in range(numero_de_intentos):
		estimacion_pi = aventar_agujas(numero_de_agujas)
		estimados.append(estimacion_pi)

	media_estimados = media(estimados)
	sigma = desviacion_estandar(estimados)
	print(f'Est = {round(media_estimados,)}, sigma={round(sigma,5)}, agujas={numero_de_agujas}')

	return (media_estimados, sigma)

def estimar_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 *= 2

	return media

if __name__ == '__main__':
	estimar_pi(0.01,1000)

Solucion sencilla del ejercicio se usa otro metodo para graficar
Los metodos mas importantes son estimar pi y crear muestra, se puede modificar para no usar el metodo graficar, incluso solo usar el metodo estimar_pi

import os
from random import uniform
from math import sqrt
from statics import Estadisticos
from bokeh.plotting import figure, output_file, show
from bokeh.models import ColumnDataSource
os.system("cls")

def estimar_pi(puntos):
    coor_points = [(uniform(-1,1), uniform(-1,1)) for _ in range(puntos)]
    in_circle = 0
    for i in coor_points:
        if sqrt(i[0]**2 + i[1]**2)<=1:
            in_circle+=1
    valor_pi = (4* in_circle) /puntos
    graficar(coor_points)
    #print(valor_pi)
    return valor_pi


def graficar(arreglo):
    in_circle = {"x_invalues": [i[0] for i in arreglo if sqrt(i[0]**2 +i[1]**2)<=1],
        "y_invalues": [i[1] for i in arreglo if sqrt(i[0]**2 +i[1]**2)<=1],
       }
    out_circle = { 
        "x_outvalues":[i[0] for i in arreglo if sqrt(i[0]**2 +i[1]**2)>1],
        "y_outvalues":[i[1] for i in arreglo if sqrt(i[0]**2 +i[1]**2)>1]
    }
    in_source = ColumnDataSource(data=in_circle)
    out_source = ColumnDataSource(data=out_circle)
    p = figure()
    p.circle(x='x_invalues', y='y_invalues', source=in_source, size=5, color="red")
    p.circle(x='x_outvalues', y='y_outvalues', source=out_source, size=5, color="navy")
    show(p)

def crear_muestra(tries):
    #Valores de pi
    pi_array = []
    for i in range(tries):
        puntos = int(input("Cuantos puntos tendra su muestra: "))
        pi_array.append(estimar_pi(puntos))


crear_muestra(4)

🫠 Una explicación un poco más profunda que me ayudó a entender el código:
.
Lo que estamos buscando es una simulación que, a determinado número de dardos e intentos, tenga una población de valores de pi de los cuales el 95% (casi todos) tengan tan sólo una diferencia de máximo 0.01 con su media. Es decir, tendríamos valores de pi con una precisión de 0.01 con un grado de confianza del 95%.

• Recordemos que ese 95% de la población puede ser representado como una diferencia de 1.96 veces la desviación estándar (σ) respecto a la media.
• Si buscáramos que el 99.7% de la población tuviera ese mismo máximo de 0.01 de diferencia con su media, estaríamos hablando de una diferencia de 3σ respecto a la media.
  .

Cada simulación de la función estimacion nos devuelve su media y su sigma (desviación estándar). Así mismo, se intuye que a mayor número de agujas, habrá menor variabilidad de los valores de π, pues ambas áreas serían más exactas.
Como mencionaba al inicio, buscamos esa simulación ideal con una precisión de 0.01 dentro de una distancia de 1.96σ de la media, y la pregunta sería ¿cómo reconocer esa simulación?, pues con su desviación estándar.
Y es que la desviación estándar de esa simulación multiplicada por 1.96 debe dar como resultado la precisión que buscamos para esa población (0.01). Es decir, 1.96σ = 0.01, o, despejando sigma, σ = 0.01 / 1.96. Es por eso que en el bucle while del código observamos que el código solo se detendrá cuando una simulación tenga una desviación estándar menor a la descrita anteriormente, que equivale aproximadamente a 0.00510204.

• El motivo de esta “ecuación” es que a 1.96 sigmas de la media debemos tener un valor igual a nuestra precisión, o 0.01 más de la media, por lo que a 1 sigma de la media deberíamos de tener un valor igual a aproximadamente la mitad de la precisión (si redondeamos ese 1.96 a 2 jajaja).

Si buscaramos que un 99.7% de la población tuviera una precisión de 0.001, la desviación estándar (σ) multiplicada por 3 debe ser igual a la precisión, o despejando sigma y escribiendo de otra forma: σ = 0.001 / 3. O escrito en el código sería while sigma >= precision / 3: # usar más agujas jaja
.
😉 ❤️ Espero que a alguien que esté tan perdido como yo le sea útil.

Resultado de la implementación de cálculo del número pi, al mirar las frecuencias de aparición de cada número es posible generar la gráfica de distribución normal. Al aumentar el número de intentos solo aumenta la cantidad de puntos, por otro lado al aumentar la cantidad de agujas lanzadas la gráfica se hace mas estrecha debido a que baja la desviacion estandar.

Una manera sencilla de calcular PI con random

https://gist.github.com/alirioangel/b55640838a75bb892e9d3a0c11116089

Dejo link de un tutorial que hice explicando de donde surge el 1.96 del 95% de confianza y más. espero les ayude.
https://platzi.com/tutoriales/1835-programacion-estocastica/9585-nivel-de-confianza-y-valor-de-z/

al dividir por 1.96 significa que el 95% de los datos estarán incluidos en el rango en el cual la media sea pi ± 0,01, adjunto los resultados del código para mostrar un poco más cómo funciona. 😃

<con 1000 agujas y una desviación de  0.01
 la media es 3.14332 con una desviación estándar de 0.05177
 con 2000 agujas y una desviación de  0.05177
 la media es 3.14305 con una desviación estándar de 0.03528
 con 4000 agujas y una desviación de  0.03528
 la media es 3.14303 con una desviación estándar de 0.02564
 con 8000 agujas y una desviación de  0.02564
 la media es 3.14105 con una desviación estándar de 0.01801
 con 16000 agujas y una desviación de  0.01801
 la media es 3.14113 con una desviación estándar de 0.0135
 con 32000 agujas y una desviación de  0.0135
 la media es 3.14186 con una desviación estándar de 0.00911
 con 64000 agujas y una desviación de  0.00911
 la media es 3.1413 con una desviación estándar de 0.00665
 con 128000 agujas y una desviación de  0.00665
 la media es 3.14162 con una desviación estándar de 0.00455>

Hola a todos, gracias a vuestros aportes me inspire para hacer una representacion grafica del metodo, es simplemente un grafico en el que muestro las agujas dentro y fuera del circulo, espero que le sirva a alguien, este es el codigo:

<import random
from varianza import desviacion_estandar
from bokeh.plotting import figure, output_file, show



numero_agujas = int(input('con cuantas agujas quieres correr la simulacion?? '))  
in_circle_x = [ ]
in_circle_y = [ ]
out_circle_x = [ ]
out_circle_y = [ ]
dentro_del_circulo = 0
    

for _ in range(numero_agujas):
    x = random.uniform(1,-1)
    y = random.uniform(1,-1)

    if (x**2 + y**2) ** 0.5 < 1:
        in_circle_x.append(x)
        in_circle_y.append(y)
        dentro_del_circulo += 1
    else:
        out_circle_x.append(x)
        out_circle_y.append(y)

    #pi = (4 * dentro_del_circulo) / numero_agujas
output_file("line.html")
plot = figure(plot_width=600, plot_height= 600)
plot.circle(in_circle_x, in_circle_y, size=5, color='blue', alpha=0.5)
plot.circle(out_circle_x,out_circle_y,size=5, color='red',alpha=0.5)
show(plot)


>

Investigando un poco la diferencia entre conclusión Estadisticamente valida vs conclusión valida, es que con las muestras, usamos las desviaciones estandar como un intervalo es decir, nosotros podemos asegurar con un 95% de probabilidad que el valor de pi esta entre la media, que en este caso me dio 3.14 y ± 1.96*desviaciones estandar

 [(3.14 - (1.96*desviacion estandar )), (3.14 + (1.96*desviacion estandar)]

Y la “conclusión valida” sería decir que el valor de pi es 3.1415…

Genial esta clase, muy buena, ya me voy sintiendo como un data scientist!

https://colab.research.google.com/drive/18Gj_5zP36Lxzt77MYKSdnoznoyR_DuPs?usp=sharing Collab de la clase

Este fue mi código para el ejercicio:

# calcular pi mediante una simulacion

from random import uniform, random, choice
from bokeh.plotting import figure, output_file, show
import os

os.system('clear')

# los puntos oscilaran tanto en X y Y de 0 ~ 10
def generarPuntos(cantidad):
  # definimos puntos como una lista para los puntos
  puntos = []
  # generamos la cantidad de puntos requerida
  for _ in range(cantidad):
    # generamos coordenadas
    # x = uniform(-1, 1)
    # y = uniform(-1, 1)
    x = random() * choice([-1, 1])
    y = random() * choice([-1, 1])
    # generamos un punto
    punto = (x, y)
    # registramos el punto
    puntos.append(punto)
  
  # y regresamos los puntos
  return puntos

# calculamos la distancia de un punto al origen
def distanciaAlOrigen(punto):
  # usaremos la formula de calcular distancia entre dos puntos
  # punto p o punto 2
  x2 = punto[0]
  y2 = punto[1]
  # calculamos la distancia
  distancia = (x2 ** 2 + y2 ** 2) ** .5
  #Nota: el origen no se usa en el calculo debido a que sus coordenadas son 0,0

  # y retornamos la distancia redondeada a enteros
  return distancia

# comprobamos que un punto este dentro de la circunferencia 
def estaDentro(punto):
  # obtenemos la distancia del punto
  distancia = distanciaAlOrigen(punto)
  # retornamos una comprobacion de si el punto esta dentro
  return distancia <= 1 if True else False

# funcion para graficar
def graficar(x_dentro, y_dentro, x_fuera, y_fuera, pi):
  # nombre del html
  output_file("calculo_pi.html")
  # creamos la grafica
  grafica = figure(plot_width=800, plot_height=800, title=f'Resultado de Pi: {pi}')
  # dibujamos los puntos de dentro
  grafica.circle(x_dentro, y_dentro, size=5, color="red", alpha=0.5)
  # dibujamos puntos de fuera
  grafica.circle(x_fuera, y_fuera, size=5, color="navy", alpha=0.5)
  # y renderizamos la grafica
  show(grafica)

# funcion principal
def main(cantidad):
  # generamos los puntos
  puntos = generarPuntos(cantidad)
  
  # registramos las X y Y de los puntos para contar y graficar
  x_dentro = []
  y_dentro = []
  x_fuera = []
  y_fuera = []

  # leemos y comparamos los puntos
  for punto in puntos:
    # puntos que estan dentro
    if estaDentro(punto):
      x_dentro.append(punto[0])
      y_dentro.append(punto[1])
    # puntos que estan fuera
    else:
      x_fuera.append(punto[0])
      y_fuera.append(punto[1])
  
  # vemos como estuvieron los puntos
  print(f'Puntos -- Dentro: {len(x_dentro)}, Fuera: {len(x_fuera)}')
  # calculamos pi: area del cuadrado * puntos del circulo / puntos del cuadrado
  pi = (4 * len(x_dentro)) / cantidad
  # graficamos los puntos que tenemos
  graficar(x_dentro, y_dentro, x_fuera, y_fuera, pi)
  # y lo retornamos
  return pi


if __name__ == "__main__":
  # pedimos una cantidad n de puntos a simular
  cantidad = int(input('Con cuantos puntos deseas simular: '))
  # calculamos pi
  pi = main(cantidad)
  # imprimimos el resultado
  print(f'Resultados: Pi = {pi}')

Muy Loco!

Muy bueno este ejemplo David

La razón por la que David usó normalidad, se debe a que la media muestral se distribuye normalmente (Teorema del Límite Central).

Hola, ayuda.
entiendo de donde sale el 1,96 pero no entiendo por qué se debe hacer la dividir

while sigma >= precision / 1.96:
        media, sigma = estimacion(numero_de_agujas, numero_de_intentos)
        numero_de_agujas *= 2

Amigos, a mi me quedo algo confuso la parte de niveles de confianza y distribución normal. Algo recordaba de las clases de la u y me puse a investigar.

Entiendo que para el caso de significancia, cuando esta es 0.01 o 1% le corresponde un valor Z(0.01/2 = 0.005) estandarizado asociado al nivel de riesgo de entre 2.56 y 2.57.
Para el caso de significancia al 0.05 o 5% le corresponde un valor Z(0.05/2 = 0.025) asociado al nivel de riesgo de 1.96.

Esto lo puede revisar en una tabla de distribución normal y ayudandome el siguiente enlace.

Sigo sin entender porque en la funcion divide a la presicion entre el 1.96(que entiendo es para conseguir una confianza del 95% y le corresponde a una significancia de 0.05)

Agradecere alguien pueda aclarar mi duda, gracias 😃

def estimar_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 *= 2
    
    return media

https://www.youtube.com/watch?v=cflk_8eLR7A

muy buen programa, implementacion de las funciones que ya aviamos hecho

En la ingeniería estructural calcular la energía que disipan los elementos de una estructura como vigas y columnas se resume a determinar el área de un dominio. En las siguiente imagen me base en el contenido del curso para resolver el problema mediante el método de Montecarlo.

¿Este un curso de estadística o de python? ¿O de ninguno? Porque el profesor se salta muchos conceptos estadísticos cuando hace los ejercicios y no explica muy bien la abstracción del código que realiza.
Tiene las mejores intenciones pero realmente no sabe realizar un curso de manera óptima para quienes son nuevos en el tema, no se puede depender de la comunidad para que el curso sea exitoso.

Encontre un error dentro del archivo de la clase, ya que aparece:

return (2 * adentro_del_circulo) / numero_de_agujas

Pero realmente es:

return (4 * adentro_del_circulo) / numero_de_agujas

Obtuve estos resultados:

Est=3.1406, sigma=0.05205, agujas = 1000
Est=3.14126, sigma=0.03632, agujas = 2000
Est=3.14058, sigma=0.02619, agujas = 4000
Est=3.14122, sigma=0.01793, agujas = 8000
Est=3.14223, sigma=0.01278, agujas = 16000
Est=3.14166, sigma=0.00904, agujas = 32000
Est=3.14158, sigma=0.00659, agujas = 64000
Est=3.14187, sigma=0.00444, agujas = 128000

Aqui mi codigo para graficar el area del circulo con r = pi, con puntos aleatorios

https://github.com/srjulianf/AREA-DE-PI-CON-AGUJAS-ALEATORIAS/blob/13f30d5989ef68c5e4bbcb88580a8fb9abe7d19b/bokeh_plot.png

https://github.com/srjulianf/AREA-DE-PI-CON-AGUJAS-ALEATORIAS/blob/2b8eccff75dceb4b236ac5c4a6ef7d8639d402e8/grafica_py.py

import random
import math
from bokeh.models.mappers import ColorMapper
from bokeh.plotting import figure, show

def graficar(x1, y1, x2 , y2 ):
    grafica = figure(title= "Area de pi con agujas" )
    grafica.scatter(x1,y1,color ="blue")
    grafica.scatter(x2,y2, color= "red")
    show(grafica)


def aventar_agujas(numero_agujas):
    dentro_del_circulo = 0
    X_blue = []
    Y_blue = []
    X_red = []
    Y_red = []

    for _ in range(numero_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) #radio 1
        if distancia_desde_el_centro <= 1:
            X_blue.append(x)
            Y_blue.append(y)
        else:
            X_red.append(x)
            Y_red.append(y)

        
    graficar(X_blue,Y_blue, X_red, Y_red)
       # if distancia_desde_el_centro <= 1:
        #    dentro_del_circulo += 1



if __name__ == '__main__' :
    aventar_agujas(10000)

La simulación por cada que sigma y el valor de PI cambian

La inferencia estadística se usa para aceptar o rechazar hipótesis

para una implementacion mas corta:
me puse como reto implementarlo con matrices de numpy

import numpy as np
from matplotlib import pyplot as plt

def estimacion_de_pi(nro_intentos, nro_agujas):
    puntos = np.random.rand(nro_intentos, nro_agujas, 2) * 2 - 1 # genera puntos aleatorias con valores del -1 a 1
    distancias_al_origen = np.linalg.norm(puntos, axis=2)
    pts_circulo = np.sum(distancias_al_origen <= 1, axis=1) # las que entraron al circulo
    pi_estimations = (pts_circulo / nro_agujas) * 4
    pi = pi_estimations.mean()
    sigma = pi_estimations.std()
    return pi, sigma

def estimar_pi(nro_intentos, precicion=0.04):
    sigma = precicion
    nro_agujas = 1000
    while sigma >=  precicion / 1.96:
        pi, sigma = estimacion_de_pi(nro_intentos, nro_agujas)
        nro_agujas *= 2
        print(f'Est: {np.around(pi, 5)}, desviacion estandar {np.around(sigma, 5)} con {nro_agujas} agujas')

if __name__ == "__main__":
	nro_intentos = 1000
	estimar_pi(nro_intentos)

Quiero recordarles que cierren la consola o cancelen la ejecución del programa. O les va a pasar esto…
![](

Una conclusión válida es super diferente a una conclusión estadísticamente válida.

Lo que nuestro algoritmo y simulación hace, no es como tal querer generar un conclusión válida, sino generar a través de una reducción de sigma por cada iteración una conclusión estadísticamente válida.

El número PI, es tan maluco de calcular que necesitamos computadoras.

No hay necesidad de usar como tal el, round, para redondear los números, desde la cadena de formato, podemos hacer esto:

print(f"tu variable {xd:.4d") 

El 4 es el valor que puede variar según tu desición

Entendemos cómo funcionan las cosas, apartir de su expresión más pequeña.

Excelente aporte que se hace con bokeh, asi podremos ver como es que funciona graficamente,

Les comparto mi solución de la práctica. La realicé con 15000 y 50000 puntos.

import random
from bokeh.plotting import figure, show, output_file

def sim(num_pts, iter, num_iter, plotting=False):
    in_circle_x = []
    in_circle_y = []
    in_circle_pts = 0
    out_circle_x = []
    out_circle_y = []
    out_circle_pts = 0

    for _ in range(num_pts):
        x = random.random() * random.choice([-1, 1])
        y = random.random() * random.choice([-1, 1])

        radio = (x**2 + y**2)**0.5
        if radio <= 1:
            in_circle_x.append(x)
            in_circle_y.append(y)
            in_circle_pts += 1
        else:
            out_circle_x.append(x)
            out_circle_y.append(y)
            out_circle_pts += 1
    
    if plotting:
        output_file("aprox_pi_montecarlo.html")
        fig = figure(plot_width = 800,
        plot_height = 800,
        title=f"Aproximación de Pi por medio de Montecarlo: {num_pts} puntos. Iteración {iter}/{num_iter}")
        fig.circle(in_circle_x, in_circle_y, color='red', alpha=0.5)
        fig.circle(out_circle_x, out_circle_y, color='green', alpha=0.5)
        show(fig)
    
    return 4*in_circle_pts/(out_circle_pts+in_circle_pts)


def avg(nums):
    return sum(nums)/len(nums)

def std_deviation(nums):
    temp_sum = 0
    average = avg(nums)
    for num in nums:
        temp_sum += (num-average)**2
    return temp_sum/len(nums)


def main(num_pts, num_iter, plotted_iters):
    sims_pi = []

    for i in range(num_iter):
        plot = False
        if i in plotted_iters:
            plot = True
        simu = sim(num_pts, i, num_iter, plot)
        sims_pi.append(simu)
        #print(simu)
    
    print(f"La media de las simulaciones es: {avg(sims_pi)}")
    print(f"Desviación estándar: {std_deviation(sims_pi)}")
    

if __name__ == '__main__':
    num_pts = int(input("Número de puntos: "))
    num_iter = int(input("Número de iteraciones: "))
    plotted_iters = [num_iter/2]

    main(num_pts, num_iter, plotted_iters)

from random import uniform
# from bokeh.plotting import figure, show

def is_inside_circle(root=(0, 0), coordinate = (1, 1), circle_radio=1.0):
    x, y = coordinate
    root_x, root_y = root
    coordinate_distance_from_root = ((root_x-x)**2 + (root_y-y)**2)**.5
    if coordinate_distance_from_root < circle_radio: return True
    return False

def put_coordinate_in_field(coordinate, rectangle_field, circle_field, circle_radio):
    if is_inside_circle(coordinate=coordinate, circle_radio=circle_radio):
        circle_field.append(coordinate)
    else:
        rectangle_field.append(coordinate)

def throw_needle(x_range=(0.0, 1.0), y_range=(0.0, 1.0)):
    x = uniform(x_range[0], x_range[1])
    y = uniform(y_range[0], y_range[1])
    return (x, y)

# def graph_field(field, color="red", p=None):
#     xs = [i[0] for i in field]
#     ys = [i[1] for i in field]

#     if p == None:
#         p = figure(plot_width=500, plot_height=500)
#     p.circle(xs, ys, size=10, color=color, alpha=0.5)
#     return p

def estimate_pi(circle_field, total_neddles):
    neddles_in_circle = len(circle_field)
    pi = 4*neddles_in_circle / total_neddles
    return pi

def main(needles, circle_radio):
    circle_field = []
    rectangle_field = []
    for _ in range(needles):
        coordinate = throw_needle(
            x_range=(-circle_radio, circle_radio),
            y_range=(-circle_radio, circle_radio)
        )
        put_coordinate_in_field(
            coordinate,
            rectangle_field, 
            circle_field,
            circle_radio
        )
    
    pi = estimate_pi(circle_field, needles)
    print(f"Estimated value for pi with {needles} neddles is {pi}")

    # p = graph_field(rectangle_field, color="red")
    # p = graph_field(circle_field, color="blue", p = p)
    # show(p)

if __name__ == '__main__':
    number_of_needles = int(input("number of needles: "))
    circle_radio = float(input("radio of the circle: "))
    main(number_of_needles, circle_radio)

Acá comparto mi solución, si quieren descomenten los comentarios para graficar con bokeh. Cabe destacar que al ser aleatorio pues no va a converger como uno quisiera.

Ya aplicando la media y la desviación estándar, y un numero n de intentos conseguí tres decimales de precisión

Resultados: Pi = 3.141252, desviacion estadar de 0.008508166429965982

este es el código, ya comentado

# calcular pi mediante una simulacion

from random import uniform, random, choice
from bokeh.plotting import figure, output_file, show
import os
from varianza import media, desviacion_estandar

os.system('clear')

# los puntos oscilaran tanto en X y Y de 0 ~ 10
def generarPuntos(cantidad):
  # definimos puntos como una lista para los puntos
  puntos = []
  # generamos la cantidad de puntos requerida
  for _ in range(cantidad):
    # generamos coordenadas
    # x = uniform(-1, 1)
    # y = uniform(-1, 1)
    x = random() * choice([-1, 1])
    y = random() * choice([-1, 1])
    # generamos un punto
    punto = (x, y)
    # registramos el punto
    puntos.append(punto)
  
  # y regresamos los puntos
  return puntos

# calculamos la distancia de un punto al origen
def distanciaAlOrigen(punto):
  # usaremos la formula de calcular distancia entre dos puntos
  # punto p o punto 2
  x2 = punto[0]
  y2 = punto[1]
  # calculamos la distancia
  distancia = (x2 ** 2 + y2 ** 2) ** .5
  #Nota: el origen no se usa en el calculo debido a que sus coordenadas son 0,0

  # y retornamos la distancia redondeada a enteros
  return distancia

# comprobamos que un punto este dentro de la circunferencia 
def estaDentro(punto):
  # obtenemos la distancia del punto
  distancia = distanciaAlOrigen(punto)
  # retornamos una comprobacion de si el punto esta dentro
  return distancia <= 1 if True else False

# funcion para graficar
def graficar(x_dentro, y_dentro, x_fuera, y_fuera, pi):
  # nombre del html
  output_file("calculo_pi.html")
  # creamos la grafica
  grafica = figure(plot_width=800, plot_height=800, title=f'Resultado de Pi: {pi}')
  # dibujamos los puntos de dentro
  grafica.circle(x_dentro, y_dentro, size=5, color="red", alpha=0.5)
  # dibujamos puntos de fuera
  grafica.circle(x_fuera, y_fuera, size=5, color="navy", alpha=0.5)
  # y renderizamos la grafica
  show(grafica)

# funcion principal
def calcularPi(cantidad):
  # generamos los puntos
  puntos = generarPuntos(cantidad)
  
  # registramos las X y Y de los puntos para contar y graficar
  x_dentro = []
  y_dentro = []
  x_fuera = []
  y_fuera = []

  # leemos y comparamos los puntos
  for punto in puntos:
    # puntos que estan dentro
    if estaDentro(punto):
      x_dentro.append(punto[0])
      y_dentro.append(punto[1])
    # puntos que estan fuera
    else:
      x_fuera.append(punto[0])
      y_fuera.append(punto[1])
  
  # vemos como estuvieron los puntos
  print(f'Puntos -- Dentro: {len(x_dentro)}, Fuera: {len(x_fuera)}')
  # calculamos pi: area del cuadrado * puntos del circulo / puntos del cuadrado
  pi = (4 * len(x_dentro)) / cantidad

  # y lo retornamos junto con los puntos para graficar
  return pi, [x_dentro, y_dentro, x_fuera, y_fuera]

def main(cantidad, intentos):
  # estimaciones de pi, y los puntos para graficar
  valores_pi = []
  puntos = []
  
  # hacemos las simulaciones requeridas
  for _ in range(intentos):
    # registramos pi, los puntos solo asumiran el valor de la ultima estimacion
    pi, puntos = calcularPi(cantidad)
    valores_pi.append(pi)
  
  # calculamos la desviacion de pi
  des = desviacion_estandar(valores_pi)

  # asumimos que pi sea la media de las estimaciones
  media_pi = media(valores_pi)

  # graficamos la ultima simulacion
  graficar(puntos[0], puntos[1], puntos[2], puntos[3], media_pi)

  # y retornamos la media y la desviacion
  return media_pi, des


if __name__ == "__main__":
  print('Recomendados 10,000 puntos')
  # pedimos una cantidad n de puntos a simular
  cantidad = int(input('Con cuantos puntos deseas simular: '))
  intentos = int(input('Cuantas veces deseas correr la simulacion: '))
  # calculamos pi
  pi, des = main(cantidad, intentos)
  # imprimimos el resultado
  print(f'Resultados: Pi = {pi}, desviacion estadar de {des}')

Muy Interesante.

```python
import random
import math as math
from estadisticas import desviacion_estandar, media
from bokeh.plotting import figure, output_file, show

def aventar_agujas(numero_de_agujas):
    
    en_circulo_x = []
    en_circulo_y = []
    fuera_circulo_x = []
    fuera_circulo_y = []
        
    for _ in range(numero_de_agujas):
        x = random.random() * random.choice([-1, 1])  #.random = valores entre 0.0 - 1.0; choice([-1,1]) = multiplica por el +1 o el -1
        y = random.random() * random.choice([-1, 1])

        distancia_desde_el_centro = math.sqrt(x**2 + y**2) # hipotenusa de la cooredenada x,y

        if distancia_desde_el_centro <= 1:
            en_circulo_x.append(x)
            en_circulo_y.append(y)
        else:
            fuera_circulo_x.append(x)
            fuera_circulo_y.append(y)
    
    datos = list(zip(en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y))
    estimacion_pi = (4 * len(en_circulo_x)) / numero_de_agujas
   
    return  estimacion_pi, datos     # dentro_de_circulo = son las agujas que cayeron dentro del circulo / numero_de_agujas = numero 
                                                      # retorna el valor de pi segun la formula pi= 4*agujas_encirculo/(agujas_cuadrado * r^2) donde r = 1

def estimacion (numero_de_agujas, numero_de_intentos): # numero_agujas = agujas lanzadas = agujas dentro del area del cuadrado
    estimados = []
    for _ in range(numero_de_intentos):                     # bucle del numero de intentos
        
        estimacion_pi, datos = aventar_agujas(numero_de_agujas)  # estimacion de pi cada vez
        estimados.append(estimacion_pi)                           

    media_estimados = media(estimados)                      # media
    sigma = desviacion_estandar(estimados)                  # desviacion standard
    
    print(f'Est={round(media_estimados, 5)}, sigma = {round(sigma, 5)}, numero de agujas= {numero_de_agujas}')
    
    return(media_estimados, sigma, datos)

def estimar_pi(precision, numero_de_intentos, numero_de_agujas):
    sigma = precision   

    while sigma >= precision /1.96:
        media, sigma, datos = estimacion(numero_de_agujas, numero_de_intentos)
        numero_de_agujas +=1000
    
    return media, datos

if __name__== '__main__':

    media, datos = estimar_pi(0.01,1000,1000)
    en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y = zip(*datos)
    output_file("line.html")
    plot = figure(plot_width=600, plot_height=600)
    plot.circle(en_circulo_x, en_circulo_y, size=5, color="red", alpha=0.5)
    plot.circle(fuera_circulo_x, fuera_circulo_y, size=5, color="navy", alpha=0.5)
    show(plot)

Muy interesante.

import random
import math as math
from estadisticas import desviacion_estandar, media
from bokeh.plotting import figure, output_file, show

def aventar_agujas(numero_de_agujas):
    
    en_circulo_x = []
    en_circulo_y = []
    fuera_circulo_x = []
    fuera_circulo_y = []
        
    for _ in range(numero_de_agujas):
        x = random.random() * random.choice([-1, 1])  #.random = valores entre 0.0 - 1.0; choice([-1,1]) = multiplica por el +1 o el -1
        y = random.random() * random.choice([-1, 1])

        distancia_desde_el_centro = math.sqrt(x**2 + y**2) # hipotenusa de la cooredenada x,y

        if distancia_desde_el_centro <= 1:
            en_circulo_x.append(x)
            en_circulo_y.append(y)
        else:
            fuera_circulo_x.append(x)
            fuera_circulo_y.append(y)
    
    datos = list(zip(en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y))
    estimacion_pi = (4 * len(en_circulo_x)) / numero_de_agujas
   
    return  estimacion_pi, datos   

def estimacion (numero_de_agujas, numero_de_intentos): # numero_agujas = agujas lanzadas = agujas dentro del area del cuadrado
    estimados = []
    for _ in range(numero_de_intentos):                     # bucle del numero de intentos
        
        estimacion_pi, datos = aventar_agujas(numero_de_agujas)  # estimacion de pi cada vez
        estimados.append(estimacion_pi)                           

    media_estimados = media(estimados)                     
    sigma = desviacion_estandar(estimados)                  
    
    print(f'Est={round(media_estimados, 5)}, sigma = {round(sigma, 5)}, numero de agujas= {numero_de_agujas}')
    
    return(media_estimados, sigma, datos)

def estimar_pi(precision, numero_de_intentos, numero_de_agujas):
    sigma = precision   

    while sigma >= precision /1.96:
        media, sigma, datos = estimacion(numero_de_agujas, numero_de_intentos)
        numero_de_agujas +=1000
    
    return media, datos

if __name__== '__main__':

    media, datos = estimar_pi(0.01,1000,1000)
    en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y = zip(*datos)
    output_file("line.html")
    plot = figure(plot_width=600, plot_height=600)
    plot.circle(en_circulo_x, en_circulo_y, size=5, color="red", alpha=0.5)
    plot.circle(fuera_circulo_x, fuera_circulo_y, size=5, color="navy", alpha=0.5)
    show(plot)
    

![](

Estos fueron mis resultados del calculo de pi:
(

import random
import math
from estadisticas import desviacion_estandar, media

def aventar_agujas(numero_de_agujas):
adentro_del_circulo = 0

for _ in range(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

return (4* adentro_del_circulo) / numero_de_agujas

def estimacion(numero_de_agujas, numero_de_intentos):
estimados = []
for _ in range(numero_de_intentos):
estimacion_pi = aventar_agujas(numero_de_agujas)
estimados.append(estimacion_pi)

media_estimados = media(estimados)
sigma = desviacion_estandar(estimados)
print(f'Est= {round(media_estimados, 5)}, sigma = {round(sigma, 5)}, agujas = {numero_de_agujas}')

return (media_estimados, sigma)

def estimar_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 *= 2

return media

if name == “main”:
estimar_pi(0.01, 1000)

import random
import math
from estadistica import estadisticas_robadas, media


def aventar_agujas(numero_de_agujas):
    adentro_del_circulo = 0

    for _ in range(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

    return (2 * adentro_del_circulo) / numero_de_agujas


def estimacion(numero_de_agujas, numero_de_intentos):
    estimados = []
    for _ in range(numero_de_intentos):
        estimacion_pi = aventar_agujas(numero_de_agujas)
        estimados.append(estimacion_pi)

    media_estimados = media(estimados)
    sigma = desviacion_estandar(estimados)
    print(f'Est={round(media_estimados, 5)}, sigma={round(sigma, 5)}, agujas={numero_de_agujas}')

    return (media_estimados, sigma)

def estimar_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 *= 2

    return media

if __name__ == '__main__':
    estimar_pi(0.01, 1000)```

Mi programa solo llega hasta los 1000 intentos 😕

Muy chevere esta clase!

Eh aquí una explicación de lo obtención de pi:

Este video expone 3 métodos, el método de Montecarlo se expone en el minuto 2:29
https://www.youtube.com/watch?v=DQ5qqHukkAc

Además que el video trae consigo 3 programas en c++ para cada método

Fuente:
QuantumFracture

Eh aqui el codigo de
QuantumFracture https://docs.google.com/document/d/1Ud_FCbVrxJoowz7rW2-FsJ3VTSoDjUQbaqvFkAxquFs/edit

Me parece una forma muy creativa de utilizar los recursos que tenemos para resolver problemas desde un punto de vista completamente diferente

Alguien sabe por que me salen esos valores, y porque el numero de agujas no se llega a multiplicar por 2 cada vez y rompe ahi?

Comparto mis resultados
Est = 3.13934, desviacion = 0.05272, agujas = 1000
Est = 3.13993, desviacion = 0.03569, agujas = 2000
Est = 3.14089, desviacion = 0.02603, agujas = 4000
Est = 3.14194, desviacion = 0.01883, agujas = 8000
Est = 3.14142, desviacion = 0.01306, agujas = 16000
Est = 3.14165, desviacion = 0.00925, agujas = 32000
Est = 3.1419, desviacion = 0.00652, agujas = 64000
y un link de una explicación desde otra perspectiva
https://www.youtube.com/watch?v=JjfrNc-G-zA

mis resultados de la simulacion:
Est=3.14128, sigma=0.05186, agujas=1000
Est=3.1411, sigma=0.03598, agujas=2000
Est=3.14104, sigma=0.02668, agujas=4000
Est=3.14258, sigma=0.01848, agujas=8000
Est=3.14169, sigma=0.01301, agujas=16000
Est=3.14189, sigma=0.00943, agujas=32000

Un poco aparte quise validar la aleatoriedad de la librería random con un código sencillo y efectivamente las pruebas siempre fueron acertadas. Valide la posibilidad de que en un rango de 10 saliera un numero y la probabilidad para cada uno siempre fue la misma. Luego valide si la probabilidad de que en un rango de 100 un numero aleatorio salga menor que 50, mayor que 50 e igual que 50… la suma de los múltiples intentos que hice siempre fue 1, luego probé con un rango mas grande y el resultado fue igual!!! alguna vez escuche que la aleatoriedad en computación no se había podido implementar al 100% alguien sabe porque es? o si simplemente fue un dato errado?

Aqui mi código:

import random

def main(numero_de_intentos):
    mayor = 0
    menor = 0
    igual = 0
    for _ in range(numero_de_intentos):
        numero = random.randrange(1000)
        if numero > 500:
            mayor +=1
        if numero < 500:
            menor +=1
        if numero == 500:
            igual +=1


    return menor, mayor, igual

  


if __name__ == "__main__":
    numero_de_intentos = 150000
    contado  = main(numero_de_intentos)
    probabilidad_menor = contado[0] / numero_de_intentos
    probabilidad_mayor = contado[1] / numero_de_intentos
    probabilidad_igual = contado[2] / numero_de_intentos
    total = probabilidad_igual + probabilidad_mayor + probabilidad_menor
    print(f'la probabilidad que sea menor a 500 es de {probabilidad_menor}')
    print(f'la probabilidad que sea mayor a 500 es de {probabilidad_mayor}')
    print(f'la probabilidad que sea igual a 500 es de {probabilidad_igual}')
    print(f'suma de las probabilidades {total}')

Resultados de simulación: Tiempo de procesamiento 9m 43s

Est=3.14153, sigma=0.0525, agujas=1000
Est=3.14114, sigma=0.03636, agujas=2000
Est=3.14278, sigma=0.02586, agujas=4000
Est=3.14121, sigma=0.01798, agujas=8000
Est=3.14192, sigma=0.0126, agujas=16000
Est=3.14167, sigma=0.00888, agujas=32000
Est=3.14171, sigma=0.00645, agujas=64000
Est=3.14141, sigma=0.00433, agujas=128000

Esto fue lo que a mi me dio:
Est=1.57139, sigma=0.02634, agujas=1000
Est=1.57087, sigma=0.01815, agujas=2000
Est=1.57069, sigma=0.01288, agujas=4000
Est=1.5709, sigma=0.00942, agujas=8000
Est=1.57089, sigma=0.00645, agujas=16000
Est=1.57068, sigma=0.00466, agujas=32000

Alguna sugerencia ?

Para simular la demanda de energía, en la universidad un compañero realizo esto para un proyecto, lo hizo en excel

por que la division precision / 1.96 ?

Entiendo el de donde salio el 1.96 , pero no entiendo la division

Este método es tan poderoso que se usa para calcular las integrales de camino de Feynman. Muy bien explicado, gracias !

no era más fácil no multiplicar por el random.choice y sencillamente multiplicar por 4?

Entendido y puesto en práctica.

Estimado = 3.14271 , sigma = 0.05179 -- Agujas : 1000
Estimado = 3.14168 , sigma = 0.03652 -- Agujas : 2000
Estimado = 3.14016 , sigma = 0.0272 -- Agujas : 4000
Estimado = 3.14117 , sigma = 0.01847 -- Agujas : 8000
Estimado = 3.14128 , sigma = 0.01275 -- Agujas : 16000
Estimado = 3.14135 , sigma = 0.0088 -- Agujas : 32000
Estimado = 3.14147 , sigma = 0.00659 -- Agujas : 64000
Estimado = 3.14182 , sigma = 0.0047 -- Agujas : 128000

acá mis resultados

Tip: En vez de utilizar el random Choice, si ustedes quieren generar números aleatorios entre 2 números, llamemoslos a y b, y no conocen alguna función que haga esto.

Pueden generarlos usando la siguiente fórmula:

Aleatorio entre a y b = (b - a) * u(0,1) + a.

Donde u(0,1) es una función que genera números aleatorios entre 0 y 1 (random.random()).

En nuestro ejemplo quedaría:

Aleatorio entre -1 y 1 = (1 - (-1)) * random.random() + (-1)
Aleatorio entre -1 y 1 = (1 + 1) * random.random() -1
Aleatorio entre -1 y 1 = 2 * random.random() -1

Espero les sea útil.

Les dejo un programa para calcular Pi de forma determinista por el método de busqueda binaria que se vió en cursos anteriores:

import math
import time

def f(x):
    return math.sin(x)

decimales = int(input('¿Con cuantos decimales quieres aproximar a pi?: '))
start_time = time.time()
resultado = 2
lim_inf = 2
lim_sup = 4
tolerancia = float('0.' + '0'*decimales + '1')
start_time = time.time()
iteraciones=0
resultado_anterior = 0
while abs(f(resultado)) > tolerancia:

    resultado_anterior = f(resultado)
    resultado = (lim_inf + lim_sup) / 2
    print(f'sin({resultado}) = {f(resultado)}')
    if f(resultado) > 0:
        lim_inf = resultado
    else:
        lim_sup = resultado
    iteraciones += 1
    if f(resultado) == resultado_anterior:
        print('\nNo es posible iterar mas')
        break
    
print(f'Aproximacion de Pi: {resultado}\n Iteraciones: {iteraciones}')
# print(start_time)
# print(time.time())
print('--- %s seconds ---' % (time.time() - start_time))

     


if __name__ == "__main__":
    pass

Tabla de dsitribución normal para allar con más o menos preción valores aleatorios:

Excelente esta clase!

Este es mi còdigo, lo hice antes de ver la clase.

import random

def dentro_circulo(puntos_aleatorios, puntos):
    count = 0
    for i in range(puntos-1):
        dot_x = puntos_aleatorios[i][0]
        dot_y = puntos_aleatorios[i][1]
        if (dot_x**2 + dot_y**2)**0.5 < 1:
            count += 1
    return count


def crear_puntos_aleatorios(puntos):
    coord_puntos = []
    for _ in range(puntos-1):
        coord_x = (random.random() - 0.5)*2
        coord_y = (random.random() - 0.5)*2
        coord_puntos.append((coord_x,coord_y))
    #print(coord_puntos)
    #print(coord_puntos[0][0]) lo usé para ver los numeros que generaba
    #print(coord_puntos[0][1])
    
    return coord_puntos


def calculo_pi(puntos, intentos):
    
    CalcPI =[]

    for i in range(intentos):
        puntos_aleatorios = crear_puntos_aleatorios(puntos)
        puntos_in = dentro_circulo(puntos_aleatorios, puntos)   
        CalcPI.append(4*puntos_in/puntos)
    
    return sum(CalcPI)/(len(CalcPI))


if __name__ == "__main__":
    num_puntos = 1000000 #1M
    intentos = 10 #1K
    
    print(f'Para {num_puntos} puntos y {intentos} intentos, el calculo de Pi es: {calculo_pi(num_puntos, intentos)}')
>

Se me ocurre que en una carrera de caballos se puede utilizar estos métodos; sin embargo, me da qué pensar que un caballo gane también depende de otras variables como la alimentación y los ejercicios.

Asi quedo.

Muy Interesante…

```python
import random
import math as math
from estadisticas import desviacion_estandar, media
from bokeh.plotting import figure, output_file, show

def aventar_agujas(numero_de_agujas):
    
    en_circulo_x = []
    en_circulo_y = []
    fuera_circulo_x = []
    fuera_circulo_y = []
        
    for _ in range(numero_de_agujas):
        x = random.random() * random.choice([-1, 1])  #.random = valores entre 0.0 - 1.0; choice([-1,1]) = multiplica por el +1 o el -1
        y = random.random() * random.choice([-1, 1])

        distancia_desde_el_centro = math.sqrt(x**2 + y**2) # hipotenusa de la cooredenada x,y

        if distancia_desde_el_centro <= 1:
            en_circulo_x.append(x)
            en_circulo_y.append(y)
        else:
            fuera_circulo_x.append(x)
            fuera_circulo_y.append(y)
    
    datos = list(zip(en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y))
    estimacion_pi = (4 * len(en_circulo_x)) / numero_de_agujas
   
    return  estimacion_pi, datos     # dentro_de_circulo = son las agujas que cayeron dentro del circulo / numero_de_agujas = numero 
                                                      # retorna el valor de pi segun la formula pi= 4*agujas_encirculo/(agujas_cuadrado * r^2) donde r = 1

def estimacion (numero_de_agujas, numero_de_intentos): # numero_agujas = agujas lanzadas = agujas dentro del area del cuadrado
    estimados = []
    for _ in range(numero_de_intentos):                     # bucle del numero de intentos
        
        estimacion_pi, datos = aventar_agujas(numero_de_agujas)  # estimacion de pi cada vez
        estimados.append(estimacion_pi)                           

    media_estimados = media(estimados)                      # media
    sigma = desviacion_estandar(estimados)                  # desviacion standard
    
    print(f'Est={round(media_estimados, 5)}, sigma = {round(sigma, 5)}, numero de agujas= {numero_de_agujas}')
    
    return(media_estimados, sigma, datos)

def estimar_pi(precision, numero_de_intentos, numero_de_agujas):
    sigma = precision   

    while sigma >= precision /1.96:
        media, sigma, datos = estimacion(numero_de_agujas, numero_de_intentos)
        numero_de_agujas +=1000
    
    return media, datos

if __name__== '__main__':

    media, datos = estimar_pi(0.01,1000,1000)
    en_circulo_x, en_circulo_y, fuera_circulo_x, fuera_circulo_y = zip(*datos)
    output_file("line.html")
    plot = figure(plot_width=600, plot_height=600)
    plot.circle(en_circulo_x, en_circulo_y, size=5, color="red", alpha=0.5)
    plot.circle(fuera_circulo_x, fuera_circulo_y, size=5, color="navy", alpha=0.5)
    show(plot)

![](![calculo_pi.png](https://static.platzi.com/media/user_upload/calculo_pi-1246e6d5-7e69-49a3-b439-032ea4233de4.jpg)

Gran ejemplo.

Realicé la actividad pero me dí cuenta que con demasiadas iteraciones para obtener mejores resultados al procesador le cuesta bastante y mi computador se traba, así que utilicé el procesamiento en paralelo, adjunto el codigo.

import random
from statistics import stdev, mean
from multiprocessing import Queue, Process

class WorkerProcess(Process):
    def __init__(self, args):
        Process.__init__(self, args=args)
        self.n = args[0]
        self.q = args[1]
        self.numProcess = args[2]
    
    def run(self):
            self.q.put(simulation(self.n,int(self.n/self.numProcess)))

def randomPoints(size):
    x = []
    y = []

    for _ in range(size):
        x_ = random.uniform(-1,1)
        y_ = random.uniform(-1,1)
        x.append(x_)
        y.append(y_)

    return x,y

def distance(x,y):
    return (x**2+y**2)**0.5

def simulation(size,n):
    piAprox = []
    for _ in range(n):
        x, y  = randomPoints(size)
        inside = []
        outside = []

        for x,y in zip(x,y):
            inside.append((x,y)) if distance(x,y)<=1 else outside.append((x,y))

        piAprox.append(4*len(inside)/size)
    pi = mean(piAprox)
    return pi

if __name__ == "__main__":
    numPoints = int(input('Cuantos puntos quieres lanzar en cada iteracion? '))
    numProcess = int(input('Cuantos procesos quieres? '))
    iteraciones = int(input(f'Cuantas veces quieres correr la simulacion? (Debe ser multiplo de {numProcess}) '))
    if iteraciones%numProcess!=0:
        raise TypeError('El numero de simulaciones no es multiplo de la cantidad de procesos.')

    q = Queue()
    for _ in range(iteraciones):
        pro = WorkerProcess(args=(iteraciones, q, numProcess))
        pro.start()

    piAprox = []
    for _ in range(iteraciones):
        piAprox.append(q.get())

    print(f'pi = {mean(piAprox)}')```

wait, no entendi en el codigo, por que se generan varias estimaciones, en que parte esta?

Estoy un poco confundido, segun yo para un nivel de confianza del 95% deberíamos estar dividiendo la precision de 0.05/1.96

y para un nivel de confianza de 99% sería de 0.01/2.58, o estoy equivocado?