¿Mi $\pi$ ¿funciona el algoritmo?

Question

Sólo me gustaría saber si este divertimento que he tenido funciona de verdad, y si funciona, ¿hay alguna forma significativa de mejorarlo? Pensé que un algoritmo basado en RNG para aproximar $\pi$ sonaba bien, así que este es mi intento de encontrar uno.

Paso 1: Sea $n = m = 0$ y elija un número entero positivo grande $N$ .

Paso 2: Seleccionar aleatoriamente un número entero $j \in [0,N^2)$ .

Paso 3: Aumentar $m$ por 1.

Paso 4: Si [ $\left\lfloor\frac{j}{N}\right\rfloor^2 + (j \mod N)^2$ ] $\leq N^2$ luego aumentar $n$ por 4.

Paso 5: $\pi \approx \frac{n}{m}$ . Vaya al paso 2.

Maximizar el número de iteraciones realizadas, así como la magnitud de $N$ debe optimizar el algoritmo.

La idea detrás del paso 4 es mapear los enteros $0$ a $N^2$ en un $N \times N$ cuadrícula; la probabilidad de que el cuadrado de la cuadrícula correspondiente se encuentre a una distancia de $N$ desde el origen es aproximadamente $\pi/4$ .

Answer 1

Su algoritmo funciona, pero siempre es bueno ver lo bien que lo hace y lo sensible que es a sus parámetros de entrada.

En primer lugar, obsérvese que para un $N$ y aumentando el número de iteraciones, el error se estanca:

Es fácil ver aquí que su algoritmo subestima sistemáticamente $\pi$ y hacerlo funcionar durante más tiempo no mejorará la situación.

Sin embargo, manteniendo constantes las iteraciones y cambiando $N$ conduce a una convergencia más rápida (pero el error sigue estancado):

Numéricamente, esto sugiere que tendrá que aumentar tanto el número de iteraciones y $N$ juntos. Intuitivamente, esto tiene sentido, a medida que el área aumenta necesitarás más iteraciones para llegar a la línea delimitadora.

Si desea aumentar la velocidad de convergencia, es posible que desee pensar en algunos de los otros parámetros: los incrementos $m \rightarrow 1, n \rightarrow 4$ probablemente debería depender de $N$ .

Buena suerte y disfrute experimentando con otros algoritmos de Monte Carlo. El código python para reproducir el primer gráfico está más abajo:

import numpy as np

def supbro(N, n_iters, n_samples):
    m = np.zeros(n_samples)
    n = np.zeros(n_samples)

    for _ in xrange(n_iters):
        j = np.random.randint(0, N**2, n_samples).astype(float)
        m += 1
        idx = (np.floor(j/N)**2 + (j % N)**2) <= N**2
        n[idx] += 4

    return n/m

N = 100
n_samples = 1000

import pylab as plt
import seaborn as sns
plt.figure(figsize=(8,4.5))

mu = []
for k in range(2, 5):
    error = np.pi-supbro(N, 10**k, n_samples)
    mu.append(error.mean())
    print "Error with N", 10**k, mu[-1]
    sns.distplot(error, label="n_iter=$10^{}$".format(k))

sns.rugplot(mu,color='k')
plt.legend()
plt.ylim(0,4)
plt.title("N={}, n_samples={}".format(N,n_samples))
plt.xlabel("supbro's algo error (from question 2653543)")
sns.despine()
plt.savefig("error_n_iter.png")