π, normally

The classic Monte Carlo party trick for π goes: throw darts uniformly at a square, inscribe a circle, count. A dart lands inside the circle with probability π/4, so four times the inside-fraction converges to π. A question on Math Stack Exchange asks for the same trick with a twist: the darts aren’t uniform, but normally distributed around the centre with σ = 0.6. Work out the probability of landing inside the circle — 0.905, straight off the Z-table — and the trick stops dead. Uniform darts measure the circle’s area, which is where the π lives; normal darts pile their mass in the middle, and the inside-share stops being an area ratio altogether. There is no π in 0.905. Nothing to solve for.

The orthodox exit is to push the samples through their own empirical CDF: they flatten back into uniforms, and the usual recipe applies. Perfectly correct — and it answers the question by undoing it. Handed bell-shaped randomness, the first move is to iron the bell flat. It seemed more interesting to leave the bell standing.

The circle never left

Here is what the dartboard framing hides: you don’t need to draw a circle, because the normal distribution already has one inside. Every one of those darts is drawn from the density

One part is the shape — the bare exponential bell, and no π anywhere in it. The other part is the normalising constant, 1/√(2π), which is nothing but π wearing a disguise. The reason it is there at all is one of the small scandals of mathematics: e^−x² has no elementary antiderivative, and the only way anyone integrates it over the whole line is by squaring the integral and switching to polar coordinates — at which point a circle appears, does its work, and leaves a √π on the table. Every bell curve is a circle that has been integrated out. The dartboard construction smuggles a circle in by hand; the darts were already carrying one.

So flip the problem. Treat the sampler as a black box that emits bell-shaped randomness. Somewhere in its output, by construction, sits 1/√(2π). If you can pin down the same quantity twice — once by calculus, which refuses to produce the π, and once by counting darts, which have it baked in — then the ratio of the two measurements is √(2π), and π falls out.

The one crack in the bell

“Pin it down by calculus” is where the Gaussian is uncooperative: the density won’t integrate in closed form, which is precisely why its probabilities live in look-up tables. But there is exactly one crack. The density itself won’t integrate; its first central moment will. Multiply the bell by (x−μ) and substitution suddenly works:

Try it by parts or by substitution — the intuition comes easiest in the clean case μ = 0, σ = 1. The right-hand side is elementary: an exponential, a σ, a constant of integration, and no π anywhere except the 1/√(2π) we dragged along on both sides for bookkeeping. That asymmetry is the entire method.

Three integrals, two of them countable

Now fix a slice [a, b] of the axis and name three integrals over it:

A is the moment of the shape — the π-free half of (1) — and thanks to (2) it collapses to a closed form you can evaluate on a pocket calculator. B and C are defined through the true density p(x), so both of them wear the 1/√(2π), and calculus can’t touch either (that is the scandal again). But darts can:

• B is a frequency. Throw n darts, count the share that lands in the slice.
• C/B is an average. It is the mean of the darts that made the cut, less μ — the sample average.

And since A and C are the same integral up to the missing constant,

and rearranging for π:

If the B in (4) looks redundant — surely it cancels? — remember that C cannot be computed directly: its analytical form is (2), and (2) needs the very π we are hunting. B and C/B are measured separately, a frequency and an average, and only their product stands in for C.

One genuine booby-trap survives. Choose the slice wisely: centre it on μ — any a, b with a−μ = −(b−μ) — and both A and C vanish by symmetry, (4) reads 0/0, and the darts teach you nothing. Any lopsided slice will do.

Two functions

The whole of pinorm is equation (2) and a dart-counter, verbatim:

from numpy import exp, inf
from scipy.stats import norm


# equation (2), evaluated between a and b: this is A — no π in sight
def firstMom(loc = 0, scale = 1, a = 0, b = 1):
    return scale * (-exp(-0.5*((b-loc)/scale)**2)--exp(-0.5*((a-loc)/scale)**2))


# the darts: returns the average C/B and the frequency B — π baked in
def monteCarloIntegration(dist, loc = 0, scale = 1, a = 0, b = 1, n = 10000):
    r = dist(loc, scale).rvs(size = n)
    r = r[(r >= a) & (r <= b)]
    return r.sum() / r.size - loc , r.size / n


if __name__ == "__main__":

    a = 0
    b = 10
    loc = -2
    scale = 5
    n = 1000000

    mean, mass = monteCarloIntegration(norm, loc, scale, a, b, n)
    moment = firstMom(loc, scale, a, b)

    res = (moment / (mass*mean))**2 / 2   # equation (4)

    print(res)   # 3.1529, that day

That is the whole repository — pinorm, as in pi-by-normal. firstMom is A in its closed form; monteCarloIntegration throws the darts and hands back the pair (C/B, B); the last line is (4). Note what is not here: no empirical CDF to estimate, no transformation back to uniformity — one frequency, one average, one squaring.

On the original numbers

Run it with the numbers from the question — μ = 0, σ = 0.6, and take the slice [0, ∞) — and the pieces come out clean. A is σ itself, 0.6, because the shape evaluates to 1 at the mean and to 0 at infinity. B is exactly ½, no calculus required: half the darts land right of the mean. The only number the darts must actually earn is the average, C/B, which a million of them put at 0.47873. Then (4) says

The bracketed ratio is √(2π) to four digits; the squaring and halving just unwraps it.

The bill

Monte Carlo error shrinks like 1/√n, and squaring the ratio at the end doesn’t do the noise any favours. A hundred darts will tell you π is somewhere between 2.3 and 6.7. A million darts buy you 3.14, give or take a hundredth. At ten million the runs are still arguing about the third decimal. The uniform dartboard pays the same tax — the price of asking randomness to do arithmetic, not a failing of the bell.

But efficiency was never the point. You can measure π by drawing a circle and laying a string along it. Getting the same number out of a million bell-shaped random numbers — asked nothing but how many landed in a slice, and where their average sits — is a different kind of measurement, and the one this script exists to make.

One coincidence, noticed only while writing this down: with the script’s defaults — μ = −2, σ = 5 and the slice [0, 10] — the average dart inside the slice, the marked dot in the first figure, sits at 3.1411, within half a thousandth of π. That is closer than the method itself typically gets with a million darts, by accident, for free. I’ve left it in.