Monte Carlo methods

Why it matters

Many physics problems are too messy for a formula: radiation passing through a detector, neutrons in a reactor, magnets near a phase transition, the path integrals of quantum field theory. Monte Carlo methods replace the exact calculation with random sampling. They are easy to write, they work in any number of dimensions, and they come with an honest error bar.

The core idea

An integral is an average times a size. If U is uniform on [a, b], then

integral from a to b of f(x) dx = (b - a) * mean of f(U)

So draw n random points, average f at those points, and multiply by (b - a). The standard error of the estimate is

error = (b - a) * s / sqrt(n)

where s is the sample standard deviation of the f values. The key fact: the error shrinks as 1 / sqrt(n). Four times as many samples halves the error. That is slow in one dimension, where grid methods do better, but the 1 / sqrt(n) rate does not get worse with dimension. A grid with 10 points per axis in 10 dimensions needs 10^10 evaluations; Monte Carlo does not care.

Worked example 1: estimating pi

Throw points uniformly in the unit square. The fraction inside the quarter circle x^2 + y^2 <= 1 approaches pi / 4.

import random

rng = random.Random(2026)
n = 100_000
hits = sum(1 for _ in range(n) if rng.random() ** 2 + rng.random() ** 2 <= 1)
print(4 * hits / n)     # about 3.14, typically within 0.01

Worked example 2: an integral with no elementary antiderivative

integral from 0 to 1 of exp(-x^2) dx = 0.746824... (it is related to the error function).

import math, random

def mc_integrate(f, a, b, n, seed=0):
    rng = random.Random(seed)
    total = total_sq = 0.0
    for _ in range(n):
        y = f(a + (b - a) * rng.random())
        total += y
        total_sq += y * y
    mean = total / n
    var = (total_sq - n * mean * mean) / (n - 1)
    return (b - a) * mean, (b - a) * math.sqrt(var / n)

est, err = mc_integrate(lambda x: math.exp(-x * x), 0, 1, 100_000, seed=7)
print(f"{est:.4f} +/- {err:.4f}")   # about 0.7468 +/- 0.0006

Always report the error with the estimate. If two runs with different seeds disagree by much more than the error bar, something is wrong.

Worked example 3: a random walk

A particle takes N steps of length 1, each in a random direction left or right. Its average position stays 0, but its mean squared distance grows linearly: <x^2> = N. That is the microscopic picture of diffusion. The typical distance grows only as sqrt(N), which is why a drop of ink spreads quickly at first and slowly later.

rng = random.Random(1)
walkers, steps = 2000, 100
msd = sum(sum(rng.choice((-1, 1)) for _ in range(steps)) ** 2 for _ in range(walkers)) / walkers
print(msd)   # close to 100

Beyond plain sampling

• Importance sampling: draw more points where f is large to reduce the variance.
• Markov chain Monte Carlo (Metropolis): sample from a probability distribution such as the Boltzmann distribution exp(-E / kT) by accepting or rejecting random moves. It is the workhorse of statistical physics.
• Random number quality: Python's random uses the Mersenne Twister, fine for these exercises. NumPy's np.random.default_rng(seed) is the modern choice for large simulations.

Common mistakes

• Not seeding. Use a fixed seed for anything you report, and change it on purpose to check stability.
• Reporting no error bar, or the standard deviation of f instead of the standard error of the mean.
• Expecting quick precision. One more correct digit costs 100 times more samples.
• Reusing the same random numbers across "independent" runs, which makes results look more certain than they are.

Recap

• Integral = (b - a) times the mean of f at uniform random points.
• Error = (b - a) s / sqrt(n); it shrinks as 1 / sqrt(n) in any dimension.
• Random walks give <x^2> = N, the root of diffusion.

Try it

Open the Do tab to write a Monte Carlo integrator that returns an estimate and its standard error.