Does my $\pi$ algorithm work?

Question

I'd just like to know if this bit of fun I had really does work, and if it does, is there a meaningful way to improve it? I thought an RNG-based algorithm for approximating $\pi$ sounded cool, so this is my attempt at finding one.

Step 1: Let $n = m = 0$, and choose some large positive integer $N$.

Step 2: Randomly select an integer $j \in [0,N^2)$.

Step 3: Increase $m$ by 1.

Step 4: If [ $\left\lfloor\frac{j}{N}\right\rfloor^2 + (j \mod N)^2$ ] $\leq N^2$ then increase $n$ by 4.

Step 5: $\pi \approx \frac{n}{m}$. Go to Step 2.

Maximizing the number of iterations performed, as well as the magnitude of $N$, should optimize the algorithm.

The idea behind Step 4 is to map the integers $0$ to $N^2$ onto an $N \times N$ grid; the probability that the corresponding grid square is within a distance of $N$ from the origin is roughly $\pi/4$.

Answer 1

Your algorithm works but it's always nice to see how well, and how sensitive it is to your input parameters.

First, note that for a fixed $N$, and increasing the number of iterations, the error stagnates:

It's easy to see here that your algorithm consistently underestimates $\pi$ and running it for longer won't improve the situation.

However, holding the iterations constant and changing $N$ leads to a more rapid convergence (but the error still stagnates):

Numerically, this suggests that you'll need to increase both the number of iterations and $N$ together. Intuitively, this makes sense, as the area increases you'll need more iterations to get to the bounding line.

If you'd like to increase the speed of your convergence, you may want to think about some of the other parameters: the increments $m \rightarrow 1, n \rightarrow 4$ should probably depend on $N$.

Good luck, and enjoy experimenting with other Monte-Carlo algorithms! The python code to reproduce the first plot is below:

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