Throwing darts
at the impossible

Chapter 1: An Average Is an Integral in Disguise

In Chapter 0 we estimated π by counting darts inside a region. That was secretly an integral — the area of the quarter-circle. Now let's make the connection explicit, because it is the engine of everything that follows.

Suppose you want the area under the curve f(x) = x² on the interval [0, 1]. Pretend you forgot the antiderivative (it is x³/3, giving the answer 1/3, but pretend). How would Monte Carlo find it?

Sample x uniformly at random on [0, 1], evaluate f(x) = x² at each point, and average those heights. That average is the integral. This is the central identity of Monte Carlo: the expected value of a function over a distribution equals the integral of that function weighted by the distribution.

Read that right to left. The sample mean (1/N)∑f(x_i) approximates the expectation, and for a uniform distribution the expectation is the integral. So averaging function values at random points estimates the integral. An average is an integral in disguise.

The deep property of this estimator is that it is unbiased. Its expected value equals the true integral exactly — for any N, even N = 1. The estimate is correctly centered on the truth; it does not systematically over- or under-shoot. That is a guarantee about the long-run average of many independent estimates, not about any single run.

This distinction trips up almost everyone, so let's be precise. Unbiased does not mean "close." A three-sample estimate can land far from 1/3. Unbiased means that if you repeated the whole three-sample experiment a million times and averaged all those estimates together, you would get 1/3 exactly. How far any single run scatters from the truth is a different question — that is variance, the subject of the next few chapters.

Worked example — three fixed points

Estimate ∫₀¹ x² dx (true value 1/3) using three fixed sample points: u = 0.1, 0.5, 0.9. Evaluate f = x² at each:

The Monte Carlo estimate is the sample mean of these heights:

The true integral is 0.3333. This particular triple reads 0.3567 — a single realization that happens to sit a bit high. Yet the estimator is unbiased: if you averaged over all possible 3-point samples drawn uniformly, the result would be exactly 1/3.

Notice we never integrated anything. We evaluated a function at three points and divided by three. That is the whole method — and it works identically for functions whose integral has no closed form at all.

What if the interval isn't [0,1]? The width factor

The clean identity above hid a convenience: on [0,1] the uniform density is exactly 1, so the average of heights is the integral. On a general interval [a, b] the uniform density is 1/(b−a), and undoing that requires multiplying the average height by the width (b−a):

The logic is "average height × width = area," exactly the rectangle picture. On [0,1] the width is 1 and silently disappears, which is why we never wrote it before. Get this factor wrong and every estimate is off by a constant multiple — a classic Monte Carlo bug.

Worked check on [1, 3]. Estimate ∫₁³ x dx (true value = [x²/2] from 1 to 3 = 4.5) with three fixed points u = 1.5, 2.0, 2.5. The mean height is (1.5 + 2.0 + 2.5)/3 = 2.0; the width is 3 − 1 = 2; so the estimate is 2.0 × 2 = 4.0. The true 4.5 is one realization away — and again, averaging over all triples would land exactly on 4.5.

This width factor generalizes to the deepest form of the identity: for any density q(x) you actually sample from, the integral of f equals the expectation of f(x)/q(x). Uniform sampling is just the special case q = 1/(b−a). Reweighting by a cleverer q is the whole subject of the importance-sampling lesson; here we stay with plain uniform draws.

# From scratch: Monte Carlo estimate of integral of x^2 on [0,1]
import numpy as np
rng = np.random.default_rng(0)

N = 100_000
x = rng.random(N)                # x_i ~ Uniform(0,1)
est = np.mean(x**2)             # average of f(x_i) = average height
print(est)                       # ~0.3333  (true integral of x^2 is 1/3)

est = np.mean(rng.random(N)**2)

And the general-interval version, with the width factor made explicit so it works on any [a, b]:

# From scratch: MC integral on a general interval [a, b]
def mc_integral(f, a, b, N, rng):
    x = rng.uniform(a, b, N)        # x_i ~ Uniform(a, b)
    return (b - a) * np.mean(f(x))  # width * average height = area

print(mc_integral(lambda x: x, 1, 3, 100_000, rng))  # ~4.5 (true integral of x on [1,3])

Chapter 2: Why the Average Settles Down

Chapter 1 promised that the sample mean estimates the integral. But why should a pile of random numbers settle onto a single true value rather than wandering forever? The answer is a theorem so important it has a name: the Law of Large Numbers (LLN).

Picture rolling a fair die over and over, tracking the running average of all rolls so far. Early on it lurches: a 2, then a 5, then a 1 — the average jumps from 2.0 to 3.5 to 2.67. But as the rolls pile up, the average creeps inexorably toward 3.5, the true mean of a fair die.

The LLN is the formal promise that this settling is not luck. It states that the running sample mean converges to the true expectation as the number of samples grows without bound. Any quantity you can write as an expectation can be estimated by averaging, and the average will find the truth.

This is the theorem that makes Monte Carlo legitimate. Without it, averaging random trials would be a hopeful guess. With it, averaging is a guaranteed-convergent procedure: throw enough darts, and you will land arbitrarily close to the answer.

Why does it work? Because the growing denominator dilutes the early noise. After a million rolls, one unlucky streak of low numbers near the start is divided by a million — its influence on the average is crushed. The mean settles not because later rolls "correct" earlier ones, but because each new roll matters less and less to the accumulated average.

That mechanism is also why the LLN does not rescue gamblers. The die is never "due" for a 6. Past rolls are never corrected. The average converges purely because the denominator grows, never because future rolls push back against past imbalance. Independence between rolls is completely untouched.

Worked example — six fixed die rolls

Take a fair six-sided die (true mean 3.5) and a fixed sequence of first rolls: 2, 5, 1, 6, 3, 4. Track the running mean after each roll:

The early values lurch between 2.000 and 3.500 — that is the random scatter of a tiny sample. But by six rolls the running mean has settled onto 3.500, and the mean of these six rolls is (2+5+1+6+3+4)/6 = 21/6 = 3.5, matching the true die mean exactly. The LLN says it stays near 3.5 forever as more rolls accumulate.

The key word is "running." Each new roll updates the average a little, and the updates shrink as N grows. That shrinkage — the funnel narrowing — is the LLN in action.

Why the denominator wins: the count can run away while the fraction settles

The subtlest part of the LLN is that it does not say the raw imbalance shrinks. It says the fraction does. These are different claims, and conflating them is the gambler's fallacy. Let's make the gap concrete with coin flips, where heads = 1, tails = 0, true mean p = 0.5.

The typical size of the head-count imbalance (heads minus N/2) grows like √N — it gets bigger, not smaller. But the fraction's deviation from 0.5 is that imbalance divided by N, so it shrinks like √N / N = 1/√N. Watch both columns move in opposite directions:

From 100 to a million flips the typical raw imbalance grows 100-fold (5 → 500), yet the fraction's error shrinks 100-fold (0.05 → 0.0005). The coin is never "due"; the deficit can widen forever. The fraction settles only because the denominator N outruns the √N imbalance. That is the engine of the LLN — and it is the same 1/√N rate we will name the "iron law" two chapters from now.

# From scratch: running mean of die rolls converges to 3.5
import numpy as np
rng = np.random.default_rng(0)

N = 100_000
rolls = rng.integers(1, 7, N)            # N rolls in {1..6}
running = np.cumsum(rolls) / np.arange(1, N+1)  # running mean after each roll
print(running[-1])                         # ~3.5

running_mean = np.cumsum(rng.integers(1, 7, N)) / np.arange(1, N+1)

And the experiment that shows the count running away while the fraction settles — the LLN's most counterintuitive face:

# From scratch: head COUNT imbalance grows, head FRACTION converges
import numpy as np
rng = np.random.default_rng(0)

N = 1_000_000
flips = (rng.random(N) < 0.5).astype(int)   # 1=head, 0=tail
heads = flips.cumsum()
imbalance = heads - np.arange(1, N+1) / 2   # heads minus expected n/2
frac_err  = heads / np.arange(1, N+1) - 0.5  # fraction minus 0.5
print(abs(imbalance[-1]), abs(frac_err[-1]))  # imbalance ~100s, frac_err ~0.0005

Chapter 3: Measuring Your Own Uncertainty

The LLN promises the average converges — eventually. But you never have infinite samples. You have 100, or 10,000. So the practical question is: how wrong might my finite estimate be? Remarkably, Monte Carlo grades its own homework.

Say you flip a coin 100 times and report a heads-probability of 0.50. How sure are you? Could the true value be 0.43? The beautiful answer is that you can compute your own uncertainty from the samples alone, without ever knowing the true value.

The quantity that captures this is the standard error — the spread of the average, written σ/√N. Here σ is the standard deviation of a single sample, and N is the number of samples. As N grows, the standard error shrinks like one over the square root of N.

Crucially, the standard error is not the same as the standard deviation. The standard deviation σ measures how spread out your raw samples are — and it does not shrink with more samples. The standard error measures how spread out the average is, and it does shrink. Confusing the two is the single most common Monte Carlo mistake.

Think of it physically. Each coin flip is wildly variable — it is 0 or 1, never 0.5. That per-flip spread σ is fixed no matter how many times you flip. But the average of many flips clusters tightly around the true probability, and that cluster tightens as you add flips. The standard error is the width of that tight cluster.

Because you can estimate σ from the same samples you already drew, you get the standard error for free. This is what lets a path tracer report "this pixel is converged to within 1%," or a trader report "this price is good to two cents" — without an oracle.

Worked example — a fair coin, N = 100

Estimate the probability a fair coin lands heads (true p = 0.5) from N = 100 flips. For a Bernoulli outcome with p = 0.5, the per-sample standard deviation is:

The standard error at N = 100 is therefore:

So your 100-flip estimate of p is good to about ±0.05. Now push to N = 10,000:

One hundred times the samples gave ten times the precision (0.05 → 0.005, a factor of exactly 10). That ratio — 100× the work for 10× the accuracy — is the exact signature of the √N law, and it is so important it gets its own chapter next.

The point of this chapter: at N = 100 you already know your answer is good to about ±0.05, computed entirely from p and N. Monte Carlo always tells you how wrong it might be.

Where the √N comes from — in three lines

The formula σ/√N is not a definition to memorize; it falls out of one fact: variances of independent things add. Let each sample X_i have variance σ². The estimate is the average, μ̂ = (1/N)∑X_i. Pull the constant 1/N out of the variance — constants come out squared:

The middle step used independence: the variance of the sum is the sum of the N individual variances, Nσ². Now take the square root to get the standard deviation of the estimate — the standard error:

That square root is the whole story. The variance of the average shrinks like 1/N (fast), but the error we actually feel is its square root, which shrinks only like 1/√N (slow). Chapter 4 is entirely about living with that square root.

Worked example — a non-binary sample (estimating σ from data)

The coin had a known σ = 0.5. Usually you don't — you estimate σ from the very samples you drew. Take five fixed function values from a Monte Carlo integration: f = [0.2, 0.9, 0.4, 0.7, 0.3]. Estimate the mean, the sample standard deviation, and the standard error.

Step 1 — sample mean. (0.2 + 0.9 + 0.4 + 0.7 + 0.3)/5 = 2.5/5 = 0.50.

Step 2 — squared deviations from the mean. (−0.3)², (0.4)², (−0.1)², (0.2)², (−0.2)² = 0.09, 0.16, 0.01, 0.04, 0.04; their sum is 0.34.

Step 3 — sample variance (divide by N−1 = 4, the unbiased "ddof=1" rule).

Step 4 — standard error.

So this five-sample run reports 0.50 ± 0.13 — and every number came from the data alone, no oracle. Dividing by N−1 instead of N is the small "Bessel" correction that keeps s² unbiased on tiny samples; with large N it makes no practical difference, but on N = 5 it matters.

# From scratch: estimate p and its standard error from samples alone
import numpy as np
rng = np.random.default_rng(0)

N = 100
samples = (rng.random(N) < 0.5).astype(float)   # 1 = heads, 0 = tails
est = samples.mean()                          # point estimate of p
se  = samples.std(ddof=1) / np.sqrt(N)        # standard error from the SAME data
print(est, se)                            # ~0.5, ~0.05

se = samples.std(ddof=1) / np.sqrt(len(samples))

Chapter 4: The Iron Law — 100× Work for 10× Accuracy

Chapter 3 showed the standard error is σ/√N. Stare at that √N and you arrive at the single most important — and most painful — fact about Monte Carlo. The error does not shrink in proportion to N. It shrinks in proportion to √N. This is the iron law.

Spell out the consequences. To halve the error, you need √N to double, which means N must quadruple. To gain a single decimal digit of accuracy — a 10× reduction in error — you need √N to grow tenfold, which means N must grow a hundredfold.

Concretely: your π estimate is accurate to one decimal after a thousand darts. You want a second decimal. The math demands a hundred thousand darts. Brute force alone can never give high precision cheaply — the returns diminish too steeply.

This is not a defect of any particular sampler or random number generator. It is baked into the variance of an average. No amount of cleverness in how you sample can change the √N rate itself. The only lever is the constant σ sitting in front — and shrinking that constant is exactly what the variance-reduction chapter is about.

Why √N and not N? Because variances of independent things add, but standard deviations do not. The variance of a sum of N independent samples grows like N; the variance of their average (dividing by N) shrinks like 1/N; and the standard deviation — the square root of variance — shrinks like 1/√N. The square root is unavoidable.

Once you internalize the iron law you stop being surprised by Monte Carlo's costs. Grainy renders, jittery option prices, slow-converging Bayesian integrals — all are the same 4×-per-halving tax, paid in samples. Every later chapter in this lesson pays it.

Worked example — the price of decimals

A Monte Carlo estimate has standard error 0.02 at N = 1000. We want to (a) halve the error to 0.01 and (b) cut it tenfold to 0.002. Because SE = σ/√N, the required N scales as 1/SE², so the factor by which N must grow is the square of the factor by which error must shrink.

(a) Halve the error. Error must shrink by a factor of 2, so N must grow by 2² = 4:

(b) Cut the error tenfold. Error must shrink by a factor of 10, so N must grow by 10² = 100:

And the often-misjudged case: doubling N from 1000 to 2000 improves the error by only:

So one more decimal digit of accuracy costs a hundredfold more work, and doubling the samples buys a mere 41% improvement. That is the defining pain of Monte Carlo — and the reason variance reduction (Chapter 10) earns its keep.

The cost ladder, rung by rung

Tabulate the same σ = 0.5 estimator at a sequence of target errors. Each row uses N = (σ/target)². Read down the table and watch N explode while the error merely halves step to step:

Every halving of the error costs exactly 4× the samples; cutting the error to a quarter (0.025 → 0.00625 would be) costs 16×. The "vs previous" column never lies: the work factor is always the square of the precision factor. There is no row where you get precision cheaply.

# From scratch: how many samples to hit a target error?
import numpy as np

sigma = 0.5            # per-sample std (estimated from a pilot run)
target_error = 0.01    # desired standard error
N_needed = int(np.ceil((sigma / target_error)**2))  # invert SE = sigma/sqrt(N)
print(N_needed)        # 2500 — note: tighten target 2x -> N grows 4x

N_needed = int(np.ceil((sigma / target_err)**2))

And the cost ladder printed rung by rung, so the 4×-per-halving tax is impossible to miss:

# From scratch: print the cost ladder for sigma=0.5
import numpy as np
sigma = 0.5
for target in [0.05, 0.025, 0.0125, 0.00625]:
    N = (sigma / target)**2
    print(f"target {target:.5f} -> N = {N:,.0f}")
# target 0.05000 -> N = 100   (each halving of target multiplies N by 4)
# target 0.02500 -> N = 400
# target 0.01250 -> N = 1,600
# target 0.00625 -> N = 6,400

Chapter 5: Error Bars from the CLT

You have a point estimate and a standard error. But a bare number with a spread is still not quite an honest report. A trader asks: "Your option price is 10.45 — but how confident are you?" To answer, you need a confidence interval, and that requires one more theorem.

The Central Limit Theorem (CLT) says that the Monte Carlo average, whatever the messy shape of the underlying samples, is approximately a bell curve (a normal distribution) centered on the true value, with width equal to the standard error. This holds even when individual samples are wildly non-normal — binary coin flips, lumpy option payoffs, anything.

Because the average is approximately normal, you can attach a calibrated interval. A 95% confidence interval is simply the estimate plus or minus 1.96 standard errors, where 1.96 is the normal multiplier that captures 95% of the bell's mass.

Different confidence levels use different multipliers: 1.645 for 90%, 1.96 for 95%, 2.576 for 99%. Higher confidence means a wider interval — you trade precision for certainty. The CLT supplies the bell shape that justifies these exact numbers.

Now the subtle part, the part that trips up even working scientists. A 95% confidence interval does not mean "there is a 95% probability the true value lies in this interval." The true value is fixed; it is not random. It is the interval that is random — it shifts every time you rerun the experiment.

The correct reading: if you repeated the whole experiment many times, building a fresh interval each time, about 95% of those intervals would contain the true value. This particular interval either does or does not — you just do not know which. "95% confidence" is a property of the procedure, not of one fixed interval.

Worked example — a 95% interval by hand

A Monte Carlo run estimates a probability as μ̂ = 0.30, with sample standard deviation s = 0.46, from N = 400 samples. Build a 95% confidence interval (z = 1.96).

Step 1 — standard error.

Step 2 — the 95% half-width.

Step 3 — the interval.

So we report: "0.30, 95% CI [0.255, 0.345]." And by the iron law of Chapter 4, halving that band to ±0.0225 would require 4× the samples — 1600 instead of 400.

Where the multipliers come from

Why 1.96 for 95% and not some other number? It is a pure property of the standard bell curve. The CLT says the estimate is approximately normal, and for a normal distribution a fixed fraction of the mass lies within a fixed number of standard deviations of the center. The multiplier z is exactly the number of standard errors that brackets the desired fraction:

Higher confidence pushes z outward to capture more of the bell's tails, so the interval widens: 90% gives ±0.038, but 99% costs ±0.059 for the same data. You are not getting more accurate — you are casting a wider net to be more often right. That is the precision-for-certainty trade, and the CLT's bell shape is what fixes these exact z-values.

One honest caveat: the CLT is an approximation that improves with N. For small N or wildly skewed integrands the bell is imperfect and these intervals are only roughly calibrated — coverage might be 92% when you asked for 95%. With the N in the thousands that real Monte Carlo runs use, the approximation is excellent.

# From scratch: 95% confidence interval for a Monte Carlo mean
import numpy as np

mu = 0.30; s = 0.46; N = 400
se = s / np.sqrt(N)            # 0.023
lo, hi = mu - 1.96*se, mu + 1.96*se   # [0.2549, 0.3451]
print(lo, hi)

from scipy import stats
lo, hi = stats.norm.interval(0.95, loc=mu, scale=s/np.sqrt(N))

The interval is a procedure, so its 95% claim is testable: build many intervals from fresh data and count how many cover the truth. It should be about 95%:

# From scratch: verify that ~95% of CIs cover the true p
import numpy as np
rng = np.random.default_rng(0)

true_p, N, trials = 0.30, 400, 10_000
covered = 0
for _ in range(trials):
    x  = (rng.random(N) < true_p).astype(float)
    se = x.std(ddof=1) / np.sqrt(N)
    lo, hi = x.mean() - 1.96*se, x.mean() + 1.96*se
    covered += (lo <= true_p <= hi)
print(covered / trials)        # ~0.95 — the procedure's coverage

Chapter 6: The Curse Grid Quadrature Can't Escape

So far Monte Carlo looks slow — that brutal 1/√N rate. Why would anyone choose it over the obvious alternative, a regular grid of evaluation points? The answer is the single deepest reason Monte Carlo exists, and it only reveals itself in high dimensions.

Here is the astonishing fact: the 1/√N rate is dimension-independent. Monte Carlo error does not care whether you integrate over a 1-dimensional interval or a 1000-dimensional cube. The dimension d simply does not appear in the error formula σ/√N.

Now contrast a grid. To keep m points along each axis in d dimensions, you need m^d points total. With just 10 points per axis in 10 dimensions, that is 10¹⁰ = ten billion evaluations — and even then the accuracy per point is poor.

The grid's accuracy scales like N^−1/d: in high d, even a vast number of grid points barely helps, because the exponent 1/d is tiny. Monte Carlo's exponent is 1/2 for every dimension. Since 1/2 > 1/d whenever d > 2, Monte Carlo's rate beats the grid for any dimension above 2 — and the cost gap is astronomical.

This is why high-dimensional integration runs on Monte Carlo and never on grids. Option pricing over 50 correlated assets, path tracing over the effectively infinite-dimensional space of light paths, Bayesian posterior integrals over thousands of parameters — all are Monte Carlo, because no grid could ever cover those spaces.

One careful caveat. Monte Carlo does not get more accurate in higher dimensions — its rate stays put at 1/√N. What changes is the competition: the grid degrades catastrophically as d grows while Monte Carlo holds steady, so Monte Carlo wins by comparison. The constant σ can still grow with d, but the rate never does.

Worked example — cost and rate as d grows

Compare grid quadrature against Monte Carlo as the dimension d rises, requiring 10 grid points per axis (m = 10). The grid cost is m^d = 10^d:

At d = 6 the grid already needs 10⁶ = 1,000,000 points; at d = 10 it needs 10¹⁰ = 10,000,000,000 — ten billion. Meanwhile the Monte Carlo rate exponent stays fixed at 0.5 for every d.

The crossover is sharp. The grid rate exponent at d = 10 is 1/10 = 0.1, far worse than MC's 0.5. Since 1/2 > 1/d for any d > 2, Monte Carlo's rate beats the grid above two dimensions — and the grid's astronomical point count makes it hopeless long before that.

Where the grid's N^−1/d rate comes from

The grid's collapse isn't asserted — it follows from counting. Lay m points along each of d axes; the total is N = m^d. Invert that to find the per-axis resolution as a function of the budget N:

A grid quadrature rule has error proportional to a power of its spacing h (for the midpoint rule, error ~ h² per dimension; even in the best case the error is tied to h, hence to N^−1/d). The exponent 1/d is the killer: as d grows, 1/d → 0, so the error barely responds to more points. Monte Carlo's exponent is frozen at 1/2 forever.

Worked crossover. Set the two rates equal to find where MC overtakes the grid. MC error ~ N^−1/2, grid error ~ N^−1/d. They cross when the exponents match: 1/2 = 1/d, i.e. d = 2. For every d > 2 the grid's exponent 1/d is smaller (worse), so MC wins on rate — and on raw point count it wins catastrophically earlier:

A million points buys a 1000-per-axis grid in 2-D but only a 2-per-axis grid in 20-D — barely two slabs per dimension. The same million Monte Carlo points give SE ~ 1/√10⁶ = 0.001 in every dimension. That is the gap that makes high-dimensional integration synonymous with Monte Carlo.

# From scratch: the SAME code and SAME 1/sqrt(N) error for ANY dimension d
import numpy as np
rng = np.random.default_rng(0)

def f(x):                      # e.g. integrate ||x||^2 over the unit cube
    return (x**2).sum(axis=1)

N, d = 100_000, 10          # d=1 or d=1000 — identical convergence rate
x = rng.random((N, d))         # N points in the d-dimensional cube
est = f(x).mean()             # a grid here would need 10**10 points!
print(est)

est = f(rng.random((N, d))).mean()   # d=1 or d=1000, identical rate

And a head-to-head: give the grid and Monte Carlo the same budget in 5-D and compare their errors against the known answer (the integral of ||x||² over the unit cube is d/3):

# From scratch: grid vs MC on the SAME budget in d=5
import numpy as np, itertools
rng = np.random.default_rng(0)
d = 5; truth = d / 3          # exact integral of sum(x_i^2)

m = 7; N = m**d                  # grid: 7 per axis -> 16,807 points
axis = (np.arange(m) + 0.5) / m       # midpoints of m bins
grid = np.array(list(itertools.product(axis, repeat=d)))
grid_est = (grid**2).sum(1).mean()

x = rng.random((N, d))             # MC: SAME N points, but random
mc_est = (x**2).sum(1).mean()
print(abs(grid_est-truth), abs(mc_est-truth))  # MC error is competitive even here

Chapter 7: Pricing an Option by Averaging Futures

Time for the payoff. You sold a European call option: the buyer has the right to purchase a stock at a fixed strike price K = 100 in one year. What is that right worth today? This is a real, money-on-the-table question — and it is an expectation, so Monte Carlo answers it.

The fair price is the average, over all possible futures, of the discounted payoff. The payoff at expiry is max(S_T − K, 0): if the stock ends above the strike, you profit by the difference; if below, the option is worthless. Average that payoff over many simulated futures, discount back to today, and you have the price.

There is one deep subtlety that makes or breaks the answer. You must simulate the stock under the risk-neutral measure, where the drift is the risk-free rate r, not the asset's real-world expected return. Plugging in the return you would actually forecast gives the wrong price — this is the heart of arbitrage-free pricing.

Under risk-neutral dynamics, the year-end stock price follows geometric Brownian motion:

Here S₀ is today's price, r the risk-free rate, σ the volatility, T the time to expiry, and Z a standard normal draw. The (r − σ²/2) term is the drift; the σ√T·Z term is the random diffusion. Each draw of Z gives one possible future.

Run thousands of these, take each payoff, discount each by e^−rT, and average. The result converges to the famous Black-Scholes analytic price — to within the standard error. This is the showcase: watch the simulated stock-price fan spread, the in-the-money paths light up, and the running average lock onto the analytic line.

Worked example — one path by hand

Price a one-year European call with S₀ = 100, K = 100, r = 0.05, σ = 0.20, T = 1. The Black-Scholes analytic anchor (computed once, exactly) is 10.4506. Now trace ONE simulated path using a fixed standard-normal draw Z = 0.5.

Step 1 — the drift term. (r − σ²/2)T = (0.05 − 0.20²/2)(1) = (0.05 − 0.02)(1) = 0.03.

Step 2 — the diffusion term. σ√T · Z = 0.20 × 1 × 0.5 = 0.10.

Step 3 — the year-end price. Add the exponents: 0.03 + 0.10 = 0.13, then

Step 4 — the payoff. max(S_T − K, 0) = max(113.8828 − 100, 0) = 13.8828.

Step 5 — discount to today. e^−rT × payoff = e^−0.05 × 13.8828 = 0.951229 × 13.8828 = 13.2058.

This single path's discounted payoff is 13.21 — well above the average. That makes sense: Z = 0.5 is an above-average draw, pushing the stock high. Average many such discounted payoffs (most with smaller or even zero payoff) and the mean converges to the Black-Scholes 10.4506.

A second path — where the option expires worthless

One path is not Monte Carlo; it is a single sample. The whole method lives in the averaging, and the average only behaves because the estimator mixes large payoffs with many zeros. Trace a second path with a below-average draw Z = −1.0 to see the other half of the distribution.

Step 1 — the exponent. drift + diffusion = (r − σ²/2)T + σ√T·Z = 0.03 + (0.20)(1)(−1.0) = 0.03 − 0.20 = −0.17.

Step 2 — the payoff. max(84.3665 − 100, 0) = max(−15.63, 0) = 0. The stock finished below the strike, so the call is worthless. Discounting zero is still zero.

Step 3 — average the two hand paths. ( 13.2058 + 0 ) / 2 = 6.6029. With N = 2 our estimate is 6.60 — far from the true 10.45, because two samples carry enormous standard error. This is the iron law from Chapter 4 in action: only as N climbs into the thousands does the running mean settle onto the analytic line.

The analytic anchor, dissected

Where does the magic number 10.4506 come from? The Black-Scholes formula is itself just the risk-neutral expectation done in closed form. For a call it reads C = S₀·N(d₁) − K·e^−rT·N(d₂), where N is the standard-normal CDF and

Plug in our numbers: the log-moneyness ln(100/100) = 0, so d₁ = (0 + (0.05 + 0.02)(1)) / 0.20 = 0.07/0.20 = 0.35, and d₂ = 0.35 − 0.20 = 0.15. Reading a normal table, N(0.35) = 0.6368 and N(0.15) = 0.5596. Then

The Monte Carlo run does not need this formula — that is the whole point — but having it lets us measure whether the simulation is right. The two factors N(d₁) and N(d₂) are roughly the chance the option finishes in the money under two related probability measures; you will meet that measure-change idea again as reweighting in importance-sampling.

How wrong is the estimate? Read the standard error

Chapter 3 taught that every Monte Carlo answer ships with its own error bar: SE = s/√N, where s is the sample standard deviation of the discounted payoffs. Run N = 200,000 paths and you measure a sample mean of about 10.4645 with s ≈ 14.7, giving

So the 95% confidence half-width is 1.96 × 0.033 ≈ 0.064, and the estimate 10.46 ± 0.06 comfortably brackets the true 10.4506. The widget reports exactly this band live; watch it shrink as √N grows. Notice the deep lesson: to halve that ±0.06 you must quadruple the paths to 800,000 — the iron law again, now with money attached.

# From scratch: Monte Carlo European call price
import numpy as np
rng = np.random.default_rng(0)

S0, K, r, sig, T, N = 100, 100, 0.05, 0.20, 1.0, 200_000
Z  = rng.standard_normal(N)                              # one draw per future
ST = S0 * np.exp((r - 0.5*sig**2)*T + sig*np.sqrt(T)*Z) # risk-neutral GBM
payoff = np.maximum(ST - K, 0)                       # call payoff
price  = np.exp(-r*T) * payoff.mean()                   # discount + average -> ~10.45
print(price)

price = np.exp(-r*T)*np.maximum(S0*np.exp((r-sig**2/2)*T+sig*np.sqrt(T)*rng.standard_normal(N))-K,0).mean()

The generality the callout promised is one line of code. Reuse the same simulated ST array and only swap the payoff function — the put for the parity check, an Asian option averaging the path, a digital that pays a fixed amount. The simulate-and-discount skeleton never changes; the payoff is the only knob.

# Reuse ST, change only the payoff. Standard error ships for free.
put     = np.exp(-r*T) * np.maximum(K - ST, 0)   # put payoff -> ~5.57 (parity check)
se      = put.std(ddof=1) / np.sqrt(N)               # its own error bar
print(f"put = {put.mean():.4f} +/- {1.96*se:.4f}")  # 5.5735 +/- 0.04

Chapter 8: Los Alamos and the Neutron's Random Walk

Where did this whole method come from? In the 1940s at Los Alamos, Stanislaw Ulam was recovering from illness, playing solitaire, and wondering about his odds of winning. Rather than solve the combinatorics, he realized he could just deal many games and count. The same thought, applied to neutrons, became Monte Carlo.

The nuclear problem was this: a neutron entering a slab of shielding material either scatters off a nucleus, gets absorbed, or leaks out a boundary. The transport equations describing this are brutal to solve analytically. So Ulam, with John von Neumann and Nicholas Metropolis, proposed: don't solve them — simulate the neutrons and count what happens.

Each neutron's history is a random walk. At every collision, a random draw decides its fate: absorbed (the history ends), leaked out (it escaped), or scattered (it continues in a new random direction). Run thousands of these histories and tally the escapes, and you have estimated the shielding probability — the very quantity no formula could deliver.

The codename "Monte Carlo," after the famous casino, came from this work — a nod to Ulam's gambling uncle and the method's reliance on chance. It is literally where the field got its name.

The crucial conceptual point: a single neutron history is cheap and crude — a few random collisions, ending in one binary verdict (escaped or not). One neutron tells you nothing about the escape probability. The probability emerges only from the fraction over thousands of independent histories.

This is the same sample-and-count idea as the dartboard, dressed in physics. The dart's "inside vs outside" becomes the neutron's "escaped vs absorbed." The answer is a probability, and a probability is estimated by a fraction, and a fraction converges by the Law of Large Numbers at the iron 1/√N rate.

Worked example — one history by hand

Use a simplified per-collision model: at each collision a neutron is absorbed with probability 0.3, leaks out with probability 0.2, or scatters (continues) with probability 0.5. We want P(eventually LEAKS) rather than being absorbed.

The analytic answer. Scattering merely postpones the verdict — it does not resolve the neutron's fate, just sends it onward. So the eventual outcome is a race between absorb (0.3) and leak (0.2). Conditioning on a resolving event:

One fixed history. Use thresholds [0, 0.3) = absorb, [0.3, 0.5) = leak, [0.5, 1) = scatter, with fixed uniform draws U = 0.6 then U = 0.1.

First draw U = 0.6 falls in [0.5, 1), so the neutron scatters and continues. Second draw U = 0.1 falls in [0, 0.3), so the neutron is absorbed. This single history ends "absorbed" — it contributes a 0 to the escape tally.

One history gives one outcome and tells you nothing about the 0.4 escape probability. Only averaging thousands of such histories — counting the fraction that leak — recovers the 0.4. The physics lives in the count.

Why "race between absorb and leak" is exactly right

The clean formula P(leak) = 0.2/(0.2+0.3) deserves a derivation, because it shows why scattering drops out. A history leaks if it survives some number of scatters and then resolves to "leak." Summing over every possible number of scatters k = 0, 1, 2, ... gives a geometric series, where each scatter carries probability 0.5:

The bracket is a geometric series with ratio 0.5, summing to 1/(1−0.5) = 2. So P(leak) = 0.2 × 2 = 0.4 — matching the race formula exactly. The same series gives P(absorb) = 0.3 × 2 = 0.6, and the two sum to 1 as they must (every history eventually resolves). Scattering only stretches the series; it never changes the ratio of the two resolving outcomes.

How many collisions does a typical history take? Each collision resolves with probability 0.2 + 0.3 = 0.5, so the number of collisions until resolution is geometric with mean 1/0.5 = 2. Histories are short and cheap — usually one or two collisions — which is exactly why you can afford to run millions of them. Cheap-but-crude per-history, precise-in-aggregate: the Monte Carlo bargain in physical clothing.

# From scratch: estimate the neutron escape probability
import numpy as np
rng = np.random.default_rng(0)

def history(rng):
    while True:                 # walk until resolved
        u = rng.random()
        if u < 0.3: return 0     # absorbed
        if u < 0.5: return 1     # leaked / escaped
        # else u in [0.5,1): scatter, loop again

esc = np.mean([history(rng) for _ in range(100_000)])
print(esc)                      # ~0.4  (analytic 0.2/(0.2+0.3))

escape_prob = np.mean([history(rng) for _ in range(N)])

A vectorized version makes the count-not-the-particle point vivid: simulate all histories' collision counts and outcomes at once, then report the escape fraction and its standard error — the same SE = σ/√N from Chapter 3, now grading a reactor-shielding estimate:

# From scratch: escape fraction + its standard error
import numpy as np
rng = np.random.default_rng(0)

N = 200_000
esc = np.zeros(N, dtype=bool)
alive = np.ones(N, dtype=bool)
while alive.any():                      # resolve every live history
    u = rng.random(N)
    absorbed = alive & (u < 0.3)
    leaked   = alive & (u >= 0.3) & (u < 0.5)
    esc[leaked] = True
    alive = alive & ~absorbed & ~leaked   # scatterers stay alive
p_hat = esc.mean()
se = esc.std(ddof=1) / np.sqrt(N)
print(p_hat, se)                       # ~0.400, SE ~0.0011

Chapter 9: Why Renders Are Grainy

Every frame in a Pixar film, and increasingly every frame in a real-time path-traced game, asks one question per pixel: of all the light bouncing around the scene, how much reaches this point on the sensor? That question is an integral — and a hard one.

It is called the rendering equation, and it integrates incoming light over the entire hemisphere of directions above a surface, accounting for light that bounced off other surfaces first. The space of possible light paths is enormous — effectively infinite-dimensional. No closed form exists.

So renderers do exactly what we have done all lesson: they throw random rays into the scene, gather the light each one carries back, and average. Each ray is a random sample of a light path; the pixel's color is the Monte Carlo estimate of that path integral.

And here is the punchline that ties the whole lesson together: when you use too few rays, the image is grainy — salt-and-pepper noise scattered across the picture. That grain is not a display bug or a rendering glitch. It is Monte Carlo estimation variance, the same σ/√N spread we have studied since Chapter 3, made visible as pixel noise.

Because it is genuine standard error, it obeys the iron law of Chapter 4. To halve the visible grain you must quadruple the rays. Going from a noisy 64 rays-per-pixel to half the noise means 256 rays-per-pixel. Clean path-traced images are expensive for exactly the reason clean π estimates are expensive.

Modern denoisers (often neural) can clean up the grain — but they do it by trading variance for bias, smoothing the image based on learned priors. The principled, unbiased fix is always the same: more rays, paid at the 4×-per-halving rate. There is no free lunch in the rate.

Worked example — four rays into one pixel

Estimate one pixel's radiance from four fixed ray samples returning normalized radiance values L = [0.2, 0.8, 0.5, 0.1]. Then ask how many rays halve the noise.

Step 1 — the pixel estimate. The mean of the four returns:

Step 2 — the sample standard deviation. The spread of the four ray returns (using the ddof=1 sample formula) is 0.3162.

Step 3 — the per-pixel standard error — this is the visible noise:

Step 4 — rays to halve the noise. By the iron law, halving SE needs 4× the rays: from 4 to 16. From a typical 64 rays you would need:

Same 1/√N iron law from Chapter 4, now wearing the clothes of pixel noise. The dartboard, the option price, and the rendered frame all pay the identical 4×-per-halving tax.

Why some pixels are grainier than others

Look at any noisy render and the grain is uneven: smooth walls converge fast, but caustics, glossy highlights, and edges of bright light sources stay speckled long after. The reason is in the numerator of SE = σ/√N. The N (rays per pixel) is the same everywhere, but σ — the spread of radiance values a pixel's rays return — varies enormously.

A pixel on a flat matte wall gets nearly the same radiance from every ray: σ is tiny, so even a handful of rays converge. A pixel that might or might not catch a tight specular glint gets mostly dark returns with the occasional blazing hit: σ is huge, so it needs vastly more rays for the same SE. The visible noise is σ/√N pixel by pixel, and high-σ pixels are the last to clean up.

Worked contrast. Two pixels, both at N = 64 rays. Pixel A (matte) has ray-return spread σ = 0.05; pixel B (glint) has σ = 0.6. Their standard errors:

Pixel B is 12× noisier than A on identical rays. To bring B down to A's noise level you would need to scale N by (0.0750/0.00625)² = 12² = 144× the rays — the iron law's square again. This is why renderers spend their ray budget adaptively, dumping rays into the high-σ pixels and sparing the easy ones.

# From scratch: Monte Carlo estimate of a pixel's radiance
import numpy as np

L = np.array([0.2, 0.8, 0.5, 0.1])   # radiance from 4 sampled light paths
pixel = L.mean()                        # 0.40 — the rendering-equation MC estimate
se    = L.std(ddof=1) / np.sqrt(len(L))  # 0.1581 — the visible grain
print(pixel, se)

pixel_color = np.mean([trace_ray(pixel, rng) for _ in range(spp)])

And the adaptive idea in code: estimate each pixel's σ from a pilot batch, then assign rays where the noise is — the same fixed budget spent smarter:

# From scratch: pixel SE drives an adaptive ray budget
import numpy as np
rng = np.random.default_rng(0)

def pixel_se(radiances):              # visible noise of one pixel
    return radiances.std(ddof=1) / np.sqrt(len(radiances))

matte = 0.4 + 0.05*rng.standard_normal(64)   # low-variance pixel
glint = 0.4 + 0.60*rng.standard_normal(64)   # high-variance pixel
se_m, se_g = pixel_se(matte), pixel_se(glint)
extra_rays = (se_g / se_m)**2          # how much more N the glint needs
print(se_m, se_g, extra_rays)        # glint ~12x noisier -> ~144x rays to match

Chapter 10: Same Answer, a Tenth of the Samples

The iron law says you cannot beat the 1/√N rate. But look again at the standard error: σ/√N. There are two factors. You cannot touch the √N. But you absolutely can shrink the constant σ sitting in front — and shrinking it is exactly as good as buying more samples, only far cheaper.

This is the art of variance reduction. Four classic techniques each exploit some structure in the problem to lower the variance of the estimator without changing what it estimates. Done right, all four are unbiased — same answer on average, tighter spread for the same N.

Antithetic variates. Pair each random draw u with its mirror 1−u. For a monotone integrand, when f(u) is high, f(1−u) is low — they are negatively correlated, so their average cancels variance. This is the technique we will compute by hand below.

Control variates. Subtract off a related quantity whose exact answer you already know, then add that known answer back. The estimator now only has to account for the small residual, which has far less variance than the original. A path tracer, for instance, can use a cheap analytic approximation of the lighting as a control and Monte-Carlo only the correction.

Stratified sampling. Partition the domain into bins and force one sample into each, guaranteeing coverage of every region instead of leaving it to chance. No region is accidentally over- or under-sampled, which tightens the estimate. A renderer stratifies the pixel and the hemisphere of directions for exactly this reason.

Common random numbers. When comparing two configurations (this reward function vs that one in an RL evaluation), reuse the same random seed across both. The shared randomness cancels in the difference, so the comparison is far less noisy than two independent runs would allow.

Worked example — antithetic variance, analytically

Estimate ∫₀¹ e^x dx (true value e−1 = 1.7183). Compare the variance of an antithetic pair-average against a plain two-sample average — no simulation, pure analysis, so it is exactly reproducible.

Step 1 — variance of a single sample. For f(U) = e^U with U uniform on [0,1]:

Step 2 — plain two-sample average. Two independent samples, averaged, have variance Var(f)/2:

Step 3 — the antithetic covariance. Pair u with 1−u. Because e^x is increasing, f(U) and f(1−U) move oppositely — the covariance is negative:

Step 4 — antithetic pair-average variance. The average of the pair has variance ( Var(f) + Cov ) / 2:

Step 5 — the reduction factor. Compare antithetic to plain, for the same two function evaluations:

The antithetic variance is 0.0323 of the plain variance — a ~31× reduction (1 / 0.0323 = 30.9). Same two evaluations, same unbiased answer, roughly one thirty-first the variance. The negative covariance did all the work.

Worked example — control variates, and the optimal coefficient

Antithetic is a clever trick, but control variates are the workhorse, so let us derive one in full. Suppose alongside f(U) = e^U we have a cheap helper g(U) = U whose exact mean we know: E[U] = 1/2. The control-variate estimator subtracts a scaled copy of the helper's error:

Because we subtract a quantity with mean zero, θ̂ is unbiased for every choice of c — its expectation is still e−1. We get to pick c purely to minimize variance. Differentiating Var(θ̂) and setting it to zero gives the optimal coefficient:

Step 1 — the helper's variance. For U uniform on [0,1], Var(U) = 1/12 = 0.0833.

Step 2 — the covariance. Cov(e^U, U) = E[U·e^U] − E[U]·E[e^U]. The first moment integrates by parts: ∫₀¹ u·e^u du = [u·e^u − e^u]₀¹ = (e − e) − (0 − 1) = 1. So Cov = 1 − (½)(e−1) = 1 − 0.8591 = 0.1409.

Step 3 — the optimal coefficient. c* = 0.1409 / 0.0833 = 1.6903.

Step 4 — the reduced variance. Plugging c* back in, the control-variate variance collapses to Var(f)·(1 − ρ²), where ρ = Corr(f, g):

Step 5 — the reduction factor. 0.2420 / 0.0039 = 61.4×. Here the correlation between e^U and U is ρ = 0.9918, so ρ² = 0.9837 — the control explains 98.4% of f's variance, and only the 1.6% residual remains to be Monte-Carlo'd. The better your free helper correlates with the target, the more variance vanishes.

Worked example — stratified sampling

The fourth technique forces coverage. Split [0,1] into two equal strata, [0, ½] and [½, 1], and place exactly one uniform sample in each, then average. Because each draw is confined to half the domain, neither can wander into the other's territory, and the chunk of variance caused by where the samples land is eliminated.

Step 1 — within-stratum variances. For f = e^x restricted to a stratum [a,b], Var = E[f²] − E[f]². On [0, ½] this is 0.0349; on [½, 1] it is 0.0949 (the function is steeper there, so more variance).

Step 2 — the stratified estimator's variance. The estimate is ½(f₁ + f₂), one sample per stratum, so

Step 3 — compare to plain. A plain two-sample average had variance 0.1210. So stratification gives 0.1210 / 0.0325 = 3.7× less variance — more modest than antithetic here, but it never hurts and it shines in higher dimensions where coverage gaps are the dominant error. The widget's bar chart shows all four side by side so you can see antithetic and control crush the noise while stratified takes a steady, reliable bite.

# From scratch: antithetic variates for integral of exp(x) on [0,1]
import numpy as np
rng = np.random.default_rng(0)

N = 10_000
u = rng.random(N // 2)              # half the draws...
anti = (np.exp(u) + np.exp(1 - u)) / 2  # ...each paired with its mirror
est = anti.mean()                   # ~1.7183, far less variance than plain MC
print(est)

est = ((np.exp(u:=rng.random(N//2)) + np.exp(1-u))/2).mean()

The other two estimators are equally short. The control variate estimates c from the same draws (using the known E[g] = ½) — in practice you never know c* exactly, so you plug in the sample covariance, which costs a negligible bias and recovers nearly the full 61×:

# Control variate: helper g(u)=u with known mean 1/2
u  = rng.random(N)
f  = np.exp(u)
g  = u
c  = np.cov(f, g)[0, 1] / np.var(g)     # estimate c* ~ 1.69
cv = f - c * (g - 0.5)               # subtract the known-mean error
print(cv.mean())                          # ~1.7183, variance ~61x smaller

# Stratified: one uniform draw forced into each of K equal bins
K     = N
edges = (np.arange(K) + rng.random(K)) / K   # one point per 1/K-wide bin
print(np.exp(edges).mean())               # ~1.7183, no bin left to chance

Chapter 11: Monte Carlo Is Eating Modern AI

The dartboard estimator you built in Chapter 0 is not a quaint numerical trick from the 1940s. It is the silent engine under a huge swath of 2024–2026 machine learning. Once you can see the sample-and-average pattern, you start spotting it everywhere.

When a probabilistic programming language like NumPyro, Pyro, or Stan computes a Bayesian posterior expectation, it is doing Monte Carlo integration — averaging over samples of the parameters. When a reinforcement-learning agent estimates a state's value by averaging episode returns, the "MC" in "Monte Carlo methods" is literally this estimator.

When a neural network reports calibrated uncertainty via MC-dropout, it is averaging many stochastic forward passes — sample and average. When a diffusion model denoises pure noise into an image, each step is a Monte Carlo move along a sampling trajectory. The image is, in a real sense, a Monte Carlo sample.

And path-traced rendering, which we met in Chapter 9, now runs in real time in game engines on consumer GPUs — billions of Monte Carlo light-path samples per frame. Far from legacy, the method is more central to computing now than at any point since Los Alamos.

From here the family branches, and each branch lifts one limitation of plain Monte Carlo. Importance sampling reweights samples drawn from a smarter proposal. MCMC lets you sample from distributions you cannot draw from directly. Quasi-Monte Carlo replaces random points with cleverly structured ones to beat the 1/√N rate. Monte Carlo Tree Search turns sampling into planning and powers game-playing and LLM reasoning.

One diffusion sampler even reframes denoising as annealed Langevin Monte Carlo — see diffusion and hamiltonian-mc for that gradient-guided lineage. The same dartboard, scaled all the way up.

Worked example — the "MC" in MC policy evaluation

Estimate a reinforcement-learning state's value as the average discounted return over episodes. One episode yields rewards r = [1, 0, 2] at times t = 0, 1, 2, with discount γ = 0.9.

Step 1 — the discounted return of episode 1. Weight each reward by γ raised to its time index:

The γ² weight on the t=2 reward is 0.9² = 0.81, so the reward of 2 contributes 1.62.

Step 2 — a second episode. Suppose another episode gives return G₂ = 1.0.

Step 3 — the Monte Carlo value estimate. Average the returns — the very same sample-and-average estimator from Chapter 1, now estimating a value function:

As more episodes accumulate, this average converges to the true state value by the Law of Large Numbers — the identical machinery as the dartboard, with episodes playing the role of darts. See policy-gradients and particle-filter for where this leads.

The same iron law, in RL clothing

This RL estimator inherits everything from this lesson, including the cost. The two-episode estimate V̂ = 1.81 has a standard error you can read off exactly as in Chapter 3. With episode returns G = [2.62, 1.0], the sample spread and standard error are:

An SE of 0.81 on an estimate of 1.81 is enormous — two episodes barely constrain the value at all. To trust the value to ±0.08 you would need (0.81/0.08)² ≈ 100× the episodes, the iron law of Chapter 4 charging its usual tax. This is precisely why RL is sample-hungry, and why the variance-reduction tricks of Chapter 10 — baselines, control variates, common random numbers across policy comparisons — are everywhere in modern RL.

# From scratch: Monte Carlo policy evaluation (average discounted returns)
import numpy as np
gamma = 0.9
episodes = [[1, 0, 2], [1]]      # two episodes of rewards

returns = [sum(gamma**t * r for t, r in enumerate(ep)) for ep in episodes]
V_est = np.mean(returns)            # 1.81 — the MC value estimate
print(returns, V_est)            # [2.62, 1.0] 1.81

V_est = np.mean([np.polyval(ep[::-1], gamma) for ep in episodes])