Random numbers and simulation

Why this matters

Every measurement has noise. A balance reads 12.503 g one time and 12.498 g the next; a voltmeter's last digit flickers; a thermocouple drifts. Before you trust a result you need to know how that noise flows through your calculation. Sometimes a formula answers that. Often it is easier, and safer, to simulate: generate thousands of plausible measurements, push each one through the calculation, and look at the spread of results. This is the Monte CarloA method that pushes thousands of simulated inputs through a calculation and studies the spread of the results. method, and NumPy makes it fast.

Core idea: a seeded generator

Modern NumPy uses a Generator object. Create one with a seedA starting number for a random generator; the same seed gives exactly the same numbers every run. so your results are reproducible:

import numpy as np

rng = np.random.default_rng(42)      # same seed -> same numbers every run

rng.random(3)                # 3 uniform numbers in [0, 1)
rng.normal(10.0, 0.2, 5)     # 5 draws, mean 10.0, standard deviation 0.2
rng.uniform(-1, 1, (2, 3))   # a 2x3 array of uniform numbers in [-1, 1)
rng.integers(1, 7, 10)       # 10 dice rolls (upper bound excluded)
rng.poisson(4.0, 1000)       # counts, such as photons per detector window

Why seed? A colleague (or you in six months) should get exactly the same figures from the same code. Without a seed, every run differs a little and nobody can tell a bug from bad luck.

Older code uses np.random.seed(0) and np.random.rand(). It still works, but the Generator API is the recommended one: it has better statistical properties and keeps separate streams independent. Pass the rng object into functions instead of relying on a hidden global state.

Simulating measurement noise

Suppose a length is truly 50.00 mm and the caliper has a random error with standard deviation 0.02 mm. Ten simulated readings:

rng = np.random.default_rng(7)
true_length = 50.00
readings = true_length + rng.normal(0, 0.02, size=10)
print(readings.round(3))
print(readings.mean(), readings.std(ddof=1))

Use ddof=1 for the sample standard deviation of measured data (it divides by n - 1). The standard errorThe uncertainty of a mean, equal to the standard deviation divided by the square root of n. of the mean is std / sqrt(n): averaging 10 readings shrinks the random error by a factor of about 3.

Monte Carlo uncertainty propagation

Worked example: density of a sample

You measure a cylinder's mass m = 25.40 ± 0.05 g and volume V = 9.40 ± 0.10 cm³ (± is one standard deviation). The density is ρ = m / V. What is its uncertainty?

By formula (first-order propagation for a quotient):

σρ / ρ = sqrt( (σm/m)^2 + (σV/V)^2 )
       = sqrt( (0.05/25.40)^2 + (0.10/9.40)^2 ) ≈ 0.0108
ρ = 2.702 g/cm³,   σρ ≈ 0.029 g/cm³

By simulation:

rng = np.random.default_rng(2026)
n = 100_000
m = rng.normal(25.40, 0.05, n)
V = rng.normal(9.40, 0.10, n)
rho = m / V                       # 100,000 densities, vectorised
print(rho.mean(), rho.std())      # about 2.702 and 0.029
lo, hi = np.percentile(rho, [2.5, 97.5])
print(lo, hi)                     # a 95% interval

The two agree. The simulation wins when the formula gets messy: non-linear functions, correlated inputs, or non-normal errors (a uniform resolution error, a skewed drift). You change one line instead of re-deriving a derivative. np.percentile gives you an interval directly from the simulated results without assuming they are normal.

How many samples?

Monte Carlo error shrinks like 1 / sqrt(n). With n = 10,000 the estimated standard deviation is good to about 1%; with n = 1,000,000 to about 0.1%. Vectorised NumPy handles a million samples in a fraction of a second, so start with 100,000 and increase if the answer still moves between seeds.

Estimating a probability

Simulation also answers "how often" questions. If a part's diameter is N(10.00, 0.03) mm and the spec is 9.95 to 10.05 mm, what fraction is out of spec?

d = np.random.default_rng(1).normal(10.00, 0.03, 200_000)
out = (d < 9.95) | (d > 10.05)
print(out.mean())          # about 0.096, the fraction of True values

The mean of a Boolean array is the fraction of True values, which is your estimated probability.

Common mistakes

• No seed. Results change every run; figures cannot be reproduced. Seed once at the top and pass rng around.
• Re-seeding inside a loop. default_rng(0) inside a loop gives the same "random" numbers every iteration.
• std() versus std(ddof=1). For a handful of measured values use ddof=1; for a large simulated population the difference is negligible.
• Confusing standard deviation with standard error. SD describes individual readings; SE (SD / sqrt(n)) describes the uncertainty of the mean.
• Too few samples. With 100 samples, a 95% interval can move noticeably between seeds.

Quick recap

• rng = np.random.default_rng(seed) then rng.normal, rng.uniform, rng.integers, rng.poisson.
• Add simulated noise to a true value to test your analysis before real data arrives.
• Monte Carlo propagation: sample every input, compute the output array, report its mean, SD and percentiles.
• The mean of a Boolean array estimates a probability.

Try it

The Do tab asks you to write a seeded noise simulator and a Monte Carlo density calculation, then compare it with the formula.