Estimating π with Monte Carlo

Overview

A course assignment for my PhD in Artificial Intelligence. Throw random darts at the unit square, count how many land inside the quarter disc, and you have an estimate of π. The more interesting part is what happens around that estimate. The error shrinks like 1/√n, and the estimator's sampling distribution is normal, exactly as the Central Limit Theorem predicts. A single reproducible Python script produces every figure, data file and the PDF report.

Technical breakdown

The maths

A point drawn uniformly in the unit square is a hit when x² + y² ≤ 1. The quarter disc has area π/4, so the hit probability is π/4 and the estimate is π̂ₙ = 4 × hits / n.

Each hit is a Bernoulli trial with p = π/4. That makes the estimator unbiased, with variance π(4 − π)/n, a standard deviation of about 1.6422/√n. Because π̂ₙ is an average of n i.i.d. bounded terms, the CLT says it's approximately Normal(π, π(4 − π)/n). That gives two testable consequences. Each extra correct decimal digit costs 100× more samples. And at a fixed n, repeating the experiment gives a bell curve centred on π whose width is predicted, not fitted.

The code

def estimate_pi(n, rng=None):
    """Monte Carlo estimate of pi from n random points in the unit square."""
    rng = np.random.default_rng() if rng is None else rng
    return 4.0 * count_inside(int(n), rng) / int(n)

count_inside generates points in blocks of 1,000,000 and tests each block with one vectorised NumPy comparison. There's no Python-level loop over points, and peak memory stays around 16 MB whether n is 10³ or 10⁷. Every random draw comes from one seeded numpy.random.default_rng instance threaded through the whole run, so the default seed reproduces every number and figure byte for byte. The only dependencies are NumPy and matplotlib.

Results

For convergence, I took one estimate at each of 25 sample sizes from 10 to 10⁷. The estimates fall inside the theoretical 95% envelope π ± 1.96·σ/√n about 95% of the time, and on log-log axes the absolute error tracks the σ/√n line with a slope of −1/2.

For the sampling distribution, I repeated the whole experiment R = 500 times at each n:

Bias is close to zero at every n, the observed spread matches the CLT within a few percent, and skew and excess kurtosis are near zero. Every Kolmogorov–Smirnov distance to the exact CLT limit is below the 5% critical value of 0.061. The underlying draws are Bernoulli, about as far from normal as a random variable gets, yet their average is already indistinguishable from a Gaussian at n = 1000. Normality is a property of the average, not of the data.