Estimating π with random darts

Here is a wonderfully sneaky way to compute π using nothing but randomness.

Throw darts uniformly at the unit square [0,1] × [0,1]. A dart lands inside the quarter circle of radius 1 when x² + y² ≤ 1. The square has area 1 and the quarter circle has area π/4, so the fraction of darts that land inside should be about π/4. Multiply by 4 and you have an estimate of π.

This is Monte Carlo estimation: answer a hard question by sampling at random. The catch is accuracy — the error only shrinks like 1/√N, so squeezing out more digits costs a lot more darts. Drag the slider and watch the estimate settle, all running in your browser as WebAssembly.

begin
    using PlutoUI, WasmMakie
end

number of darts = 4000

random seed = 7

let
    # throw `ndarts` darts; after each one, record the running estimate of pi so we
    # can watch it converge. The randomness is an inline Park-Miller generator (pure
    # integer arithmetic), so it compiles straight to WebAssembly.
    s = (seed * 2654435761 + 12345) % 2147483647
    if s == 0
        s = 1
    end
    xs = Vector{Float64}(undef, ndarts)
    ys = Vector{Float64}(undef, ndarts)
    inside = 0
    for k in 1:ndarts
        s = (s * 16807) % 2147483647
        x = Float64(s) / 2147483647.0
        s = (s * 16807) % 2147483647
        y = Float64(s) / 2147483647.0
        if x * x + y * y <= 1.0
            inside = inside + 1
        end
        xs[k] = Float64(k)
        ys[k] = 4.0 * Float64(inside) / Float64(k)   # running estimate of pi
    end
    fig = Figure(size = (600, 340))
    ax = Axis(fig[1, 1])
    lines!(ax, [1.0, Float64(ndarts)], [3.14159265, 3.14159265])  # the true value of pi
    lines!(ax, xs, ys)                                            # our running estimate
    fig
end

mc_stats = let
    # the final estimate, computed in a SEPARATE bond-dependent cell so the markdown
    # below interpolates it live (values inside a markdown cell's own `let` bake to
    # the slider defaults).
    s = (seed * 2654435761 + 12345) % 2147483647
    if s == 0
        s = 1
    end
    inside = 0
    for k in 1:ndarts
        s = (s * 16807) % 2147483647
        x = Float64(s) / 2147483647.0
        s = (s * 16807) % 2147483647
        y = Float64(s) / 2147483647.0
        if x * x + y * y <= 1.0
            inside = inside + 1
        end
    end
    est = 4.0 * Float64(inside) / Float64(ndarts)
    err = est - 3.14159265
    if err < 0.0
        err = -err
    end
    (floor(est * 100000.0) / 100000.0, floor(err * 100000.0) / 100000.0)
end;

Estimate from 4000 darts: 3.14 (true value 3.14159...). Off by about 0.00159. Quadruple the darts and the error only roughly halves – that is the 1/sqrt(N) law.

Why the error shrinks so slowly

Each dart is an independent yes/no trial (inside or outside), so the count inside behaves like a sum of coin flips. The standard error of such an average falls like 1/√N: to get one more correct decimal digit (10× less error) you need 100× more darts. That is why Monte Carlo is fantastic for rough answers to otherwise impossible problems — high-dimensional integrals, physics simulations — but a poor way to compute π to many digits.

The same principle powers Monte Carlo methods across data science: estimate something hard by averaging many easy random samples, and know that your uncertainty shrinks like one over the square root of the work you do.

Appendix

The original lecture uses Julia's rand and Plots.jl. WebAssembly can't run those in the browser, so the darts come from an inline Park-Miller generator and the picture from WasmMakie. The idea — area by sampling, error like 1/√N — is exactly the same.