Using Monte Carlo to estimate π

When you have to solve a computationally intensive problem, you can try to solve it using a Monte Carlo method. We draw a random sample and use the information to estimate a parameter we are interested in. In the following, we estimate the value of the circle number $\pi$. The area of the circle is $A_{\circ} = \pi r^2$. The area of a square is $A_{\square} = (2r)^2$, while $r$ is the radius of the circle, and therefore $2r$ is the length of the edge of the square. The ratio between these areas multiplied by four is $\pi$. This ratio describes the probability of picking a random point $(x,y)$ within the square that is also in the circle. When we create a sample of $N$ random points within the square, $\frac{\pi}{4}$ of them need to be in the circle. Therefore, we can estimate the number $\pi \approx 4\frac{N_{\circ}}{N_{\square}}$.

n_samples = int(2e06)

def estimate_pi_On(n_samples) -> float:
    inside_circle = 0
    for _ in range(n_samples):
        x, y = np.random.uniform(-1, 1, 2)
        if x**2 + y**2 <= 1:
            inside_circle += 1
    return 4 * inside_circle / n_samples

The code above uses $2,000,000$ samples to estimate $\pi$ as $3.141652$. However, the function estimate_pi_On() takes $10.3$ seconds. By vectorizing the function, the estimation takes only $0.09$ seconds.

def estimate_pi_O1(n_samples) -> float:
    x = np.random.uniform(-1, 1, n_samples)
    y = np.random.uniform(-1, 1, n_samples)
    inside_circle = (x**2 + y**2) <= 1
    return 4 * np.sum(inside_circle) / n_samples