You’ve probably heard that March 14 is unofficially (or maybe officially in some circles) recognized as Pi Day. That’s because the digits 3-1-4 match the first three digits in the number pi (π = 3.14159265…). I’ve always said there’s no better way to celebrate made-up holidays than by doing math, so let’s get started.

One way to estimate pi is to analyze the area-ratio of a square and a circle inscribed within it. There’s a bit to unpack from that statement so let me add a visual aid:

Here we have a circle circumbscribed by a square. That’s a fancy way of saying the square’s overall width is the same as the circle’s overall width. In more formal terms, the circle’s radius is length R and the square’s sides are length 2R.

A circle’s area is defined as π*r², where r is the radius. Notice that this expression includes pi, so if we’re trying to estimate pi then we’re starting to get somewhere. A square’s area is even easier to calculate: [side]*[side].

Now look at what happens if we divide the above circle’s area by the area of its square counterpart:

The R² terms cancel and we’re left with a constant: π/4. We can then simply multiply by 4 to solve for pi.

So where does Python come in? Well, assuming for some reason that we can’t draw the shapes and measure them, the next-best approach might be to randomly generate points within the square and check how many fall within the circle.

If the points are uniformly distributed across the x and y axes, the probability of a point falling within a shape is proportional to the shape’s area. If, as we saw above, [Area of Circle] / [Area of Square] is equal to π/4, then π/4 percent of the points will fall within the circle. That’s about 78.5%. We can run a simulation to randomly generate thousands of points, calculate this ratio, and then multiply by 4 to estimate the value of pi.

Begin by randomly generating 10,000 points within a square of size 2R. This should be enough data to provide a reasonable estimate.

We’ll locate the center of the figure at (0, 0) and set R equal to 1.0. This will keep the code simple, but the final result would be no different if we used a larger square centered at (8, -4.5), for example.

import matplotlib.pyplot as plt
from random import uniform
from math import sqrt, pi

R = 1.0

scatter_x, scatter_y = [], []

for _ in range(10000):
    scatter_x.append(uniform(-R, R))
    scatter_y.append(uniform(-R, R))

With data in hand, check how many points fall within an inscribed circle of radius R. The easiest approach is to test whether a point lies within 1.0 units of the center (0, 0).

inside_circle = 0

for x, y in zip(scatter_x, scatter_y):
    distance = sqrt(x**2 + y**2)
    if distance < R:
        inside_circle += 1

Now we know how many points are within the circle. And we know all 10,000 points fall within the square because that’s how we designed the simulation.

inside_square = len(scatter_x)

These numbers are analogous to the area of each shape, since the probability of a point landing within a shape is equal to its share of the entire figure.

So we have estimates of each shape’s area, and we know their ratio should be equal to π/4. We’re almost there!

ratio = inside_circle / inside_square

pi_estimate = ratio * 4

And we’re there! Just like that.

Before moving on we should calculate how close the simulation came to pi’s true value.

error_percent = (pi_estimate - pi) / pi * 100

Rather than printing the estimate to screen, let’s display everything on a nice plot.

Because the Matplotlib code is so simple I’ll skip the commentary. But as always, the full code will be available at the end of this post.

Without further ado, our estimate of pi:

With 10,000 randomly generated points we came within 0.42% of the actual value (3.14159265…). That’s not bad! Maybe we couldn’t build a particle accelerator with that sort of precision but we could certainly write blog posts all day.

Obviously with more data comes more precision. I put together an animated GIF to show how an estimate emerges from a larger simulation. Error generally trends toward zero but it’s not a straight path!

Download the Matplotlib style.

Full code:

import matplotlib.pyplot as plt
from random import uniform
from math import sqrt, pi


R = 1.0

scatter_x, scatter_y = [], []

for _ in range(10000):
    scatter_x.append(uniform(-R, R))
    scatter_y.append(uniform(-R, R))

inside_circle = 0

for x, y in zip(scatter_x, scatter_y):
    distance = sqrt(x**2 + y**2)
    if distance < R:
        inside_circle += 1

inside_square = len(scatter_x)

ratio = inside_circle / inside_square    # = (pi*R²) / (4*R²) = pi/4

pi_estimate = ratio * 4

error_percent = (pi_estimate - pi) / pi * 100

plt.style.use("wollen_dark.mplstyle")
fig, ax = plt.subplots(figsize=(8, 8))

ax.scatter(scatter_x, scatter_y, alpha=0.3, zorder=1)

ax.add_patch(plt.Circle((0, 0), R, facecolor="None", edgecolor="white", linewidth=3.0))
ax.add_patch(plt.Rectangle((-R, -R), R*2, R*2, facecolor="None", edgecolor="white", linewidth=3.0))

results_message = f"n = {len(scatter_x):,}\n" \
                  f"pi (est.):  {pi_estimate:.4f}\n" \
                  f"error = {error_percent:+.2f}%"
ax.text(0, 0, results_message,
        font="Ubuntu Condensed", size=13, ha="center", va="center",
        bbox={"boxstyle": "round", "facecolor": "#303030", "alpha": 0.9, "pad": 0.35})

plt.savefig("pi_estimate.png", dpi=200)