Lab 03: Random numbers, seeding, and reproducible simulation

Lab 03: Random numbers, seeding, and reproducible simulation#

  • Generate numbers. Then generate them again.

  • Keep track of what you ran.

Open your own copy so you can change the code. From the lab page, use the rocket menu to open DataHub or Colab. In Colab, choose File → Save a copy in Drive, then click Connect. You don’t need a GPU.

We’ll use NumPy to generate random numbers and Matplotlib to see what they do. Run the cells from the top. When we repeat a calculation, pay attention to which cells we rerun.

Get the packages ready#

  • Import before using a package.

In Colab or your own local notebook, if an import below gives ModuleNotFoundError, remove the # from this install command and run it. If you’re asked to restart the kernel or runtime, do that, then run the imports again. Put the # back afterward so Run All doesn’t repeat the installation.

On DataHub, these packages should be in the course environment. If they’re missing, check which environment you selected and ask me before installing over it.

# %pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

Make a random number generator#

  • The seed sets the starting point.

  • The generator keeps its place.

These are pseudorandom numbers: an algorithm produces a sequence that we use as random samples. A seed initializes that algorithm. It isn’t a measure of how random or accurate a result is.

default_rng creates a generator. We’ll name it rng and ask it for numbers. Keep the seed fixed for now.

seed = 20261006
rng = np.random.default_rng(seed)

  • Five numbers between zero and one.

random(5) returns an array of five floats in \([0, 1)\): zero is included, one is excluded. This is uniform sampling; equal-width intervals have equal probability.

first = rng.random(5)
first

  • Same generator. Next five numbers.

Predict what happens before running this cell. We haven’t created a new generator.

second = rng.random(5)
second

Start again from the seed#

  • Drawing advances the state.

  • Recreating the generator starts it over.

Run only the previous cell again. You get another draw because rng remembers where it is in the sequence. Now create a new generator with the original seed.

replay_rng = np.random.default_rng(seed)
replayed = replay_rng.random(5)
replayed

  • Compare every entry.

With the same environment and calls, the same starting state reproduces the same values. Here exact equality is the check we want.

np.array_equal(first, replayed)

  • Rerun the draw cell. Then rerun the setup too.

If that comparison is False, rerun the cell that creates rng, then the cell that defines first, then the two replay cells. Running first = rng.random(5) twice without resetting rng replaces first with a later part of the sequence.

Change seed and repeat that sequence of cells. You should reproduce the new sequence. Omitting the seed with np.random.default_rng() asks NumPy to initialize from fresh system entropy; it doesn’t give you a seed you’ve written down.

Choose the distribution#

  • Choose the values your model needs.

  • Draw samples. Then look at their shape.

A random number isn’t always a number between zero and one. We might need a position, a measurement error, a waiting time, or a count. The distribution describes which values are possible and how likely they are.

We’ll use a separate generator for these examples. To reproduce this whole section, rerun its setup cell and work down in order. To change only a plot, rerun the plotting cell; we don’t need to draw new numbers.

distribution_rng = np.random.default_rng(31415)

Uniform: sample evenly across an interval#

  • Set the lower and upper limits.

  • Equal-width intervals have equal probability.

Use a uniform distribution to sample positions evenly along a line or choose an angle in a plane with no preferred direction. Here we’ll draw 5000 values between \(-1\) and \(1\). size is the number of samples, not the number of histogram bars.

Uniform is a modeling choice. Knowing only a minimum and maximum doesn’t tell us that all parts of the interval are equally likely.

uniform_samples = distribution_rng.uniform(low=-1.0, high=1.0, size=5000)
uniform_samples[:8]

Build a histogram#

  • Divide the horizontal axis into bins.

  • Count how many values land in each bin.

A histogram shows where the samples collect. The horizontal axis is the sampled value; each bar covers an interval. With the default settings, the bar height is the number of samples in that interval. It doesn’t show the order we drew them.

First choose the bin edges. Eleven edges make ten bins. linspace includes both endpoints, so these edges cover the whole interval from \(-1\) to \(1\).

uniform_edges = np.linspace(-1, 1, 11)
uniform_edges

  • Give hist the samples and the bin edges.

  • Label what each axis counts.

As in Lab 02, plt.subplots gives us a figure and an axes object. ax.hist counts the samples and draws the bars. edgecolor outlines them so we can see the bins. The returned values are the counts, edges, and bar objects; we’ll inspect the counts next.

fig_uniform, ax_uniform = plt.subplots(figsize=(6, 4), constrained_layout=True)
counts, edges, bars = ax_uniform.hist(uniform_samples, bins=uniform_edges, edgecolor="white")
ax_uniform.set(xlabel="Sample value", ylabel="Number of samples", title="Uniform: ten bins")
plt.show()

  • Ten bins. About 500 samples in each.

  • Add the counts: did we include every sample?

We drew 5000 values. Each of these ten equal-width bins has probability \(1/10\), so we expect about 500 values per bin. The counts fluctuate; they don’t have to match. Bins include their left edge and exclude their right edge, except that the last bin includes both edges.

print("Counts:", counts)
print("Total counted:", int(counts.sum()))
print("Samples drawn:", len(uniform_samples))

Change the bins, keep the samples#

  • Same numbers. Different picture.

  • More bins. Fewer samples per bin.

An integer such as bins=8 asks Matplotlib for eight equal-width bins. range sets their endpoints; samples outside that range would not be counted. Our uniform samples all lie within this range.

Here we plot the same array twice. Narrower bins contain fewer samples on average, so the counts and their vertical scales differ. Changing the bins doesn’t change the data or make the generator better.

fig_bins, axes_bins = plt.subplots(1, 2, figsize=(10, 4), constrained_layout=True)
for ax, bin_count in zip(axes_bins, [8, 40]):
    ax.hist(uniform_samples, bins=bin_count, range=(-1, 1), edgecolor="white")
    ax.set(xlabel="Sample value", ylabel="Number of samples", title=f"{bin_count} bins")
plt.show()

Discrete uniform: roll a fair die#

  • The upper limit is excluded.

  • Put one bin around each integer.

Use integers when equally likely outcomes are whole numbers: a die face or an index chosen from an array. For a six-sided die, use low=1 and high=7. Using high=6 would leave out the six.

die_rolls = distribution_rng.integers(low=1, high=7, size=5000)
die_rolls[:12]

  • Half-integer edges center the bars on the outcomes.

Edges at 0.5, 1.5, …, 6.5 give one bin for each face. Arbitrary continuous bins can split these integer outcomes awkwardly or leave empty bars between them. The dashed line is the expected count per face, not a requirement on this run.

die_edges = np.arange(0.5, 7, 1)
fig_die, ax_die = plt.subplots(figsize=(6, 4), constrained_layout=True)
ax_die.hist(die_rolls, bins=die_edges, edgecolor="white", rwidth=0.85)
ax_die.axhline(len(die_rolls) / 6, color="black", linestyle="--", label="Expected per face")
ax_die.set(xlabel="Die face", ylabel="Number of rolls", xticks=np.arange(1, 7))
ax_die.legend()
plt.show()

Normal: model fluctuations around a mean#

  • loc is the mean.

  • scale is the standard deviation.

The normal distribution is the familiar bell curve. It is often a useful approximation for measurement noise when many small, roughly independent effects add together. It can produce negative values, so it isn’t a good default for a nonnegative count or waiting time.

scale is the standard deviation, not the variance. We’ll start with mean zero and standard deviation one. A finite sample won’t have exactly those values.

normal_mean = 0.0
normal_std = 1.0
normal_samples = distribution_rng.normal(loc=normal_mean, scale=normal_std, size=5000)
normal_samples[:8]

  • Measure the sample we actually drew.

ddof=1 uses the sample-variance denominator \(N-1\) when estimating the spread. This standard deviation describes the spread of individual draws.

print("Sample mean:", np.mean(normal_samples))
print("Sample standard deviation:", np.std(normal_samples, ddof=1))

Counts and probability density#

  • Counts add to the sample size.

  • Density times bin width gives a fraction.

So far, the vertical axis has counted samples. density=True instead divides each count by the total number counted and by the bin width. The histogram’s total area becomes one. A bar’s height alone isn’t its probability; multiply by its width to get the fraction of samples in that bin.

This lets us compare the histogram with a probability density curve. We’ll calculate the normal density using the same mean and standard deviation as the sampling call.

x_grid = np.linspace(normal_mean - 4 * normal_std, normal_mean + 4 * normal_std, 301)
normal_density = np.exp(-0.5 * ((x_grid - normal_mean) / normal_std)**2)
normal_density /= normal_std * np.sqrt(2 * np.pi)

fig_normal, ax_normal = plt.subplots(figsize=(6, 4), constrained_layout=True)
density, normal_edges, normal_bars = ax_normal.hist(
    normal_samples, bins=35, density=True, alpha=0.6, edgecolor="white", label="Samples"
)
ax_normal.plot(x_grid, normal_density, color="black", label="Normal density")
ax_normal.set(xlabel="Sample value", ylabel="Probability density", title="Normal")
ax_normal.legend()
plt.show()

  • Check the area, not the sum of the heights.

np.diff subtracts adjacent edges to give each bin’s width. Multiply by the heights and add the areas. The result should be one, up to rounding.

bin_widths = np.diff(normal_edges)
print("Total histogram area:", np.sum(density * bin_widths))

Fewer samples, same distribution#

  • Keep the bins fixed for this comparison.

  • Normalize before comparing different sample sizes.

Compare the first 100 values with all 5000. The smaller sample is part of the larger one, so these aren’t independent samples. Both histograms use the same edges and density scale; otherwise the larger sample would simply have more counts. We use all the samples, including those in the tails.

fig_sizes, axes_sizes = plt.subplots(1, 2, figsize=(10, 4), sharex=True, sharey=True,
                                     constrained_layout=True)
for ax, sample in zip(axes_sizes, [normal_samples[:100], normal_samples]):
    ax.hist(sample, bins=normal_edges, density=True, alpha=0.6, edgecolor="white")
    ax.plot(x_grid, normal_density, color="black")
    ax.set(xlabel="Sample value", ylabel="Probability density", title=f"N = {len(sample)}")
plt.show()

  • Change the mean. Then change the width.

Try normal_mean = 2.0, then normal_std = 0.5. Rerun the normal sampling, statistics, density, and plotting cells. Which change shifts the histogram? Which squeezes it? Keep the standard deviation positive for this density formula.

A rough histogram doesn’t mean the generator is broken. A smooth histogram alone doesn’t prove it is good either; it doesn’t tell us whether successive draws are related.

Exponential: time until the next event#

  • Waits are nonnegative.

  • scale is the mean waiting time.

Use an exponential distribution for waiting times between events in a process with independent arrivals at a constant average rate, such as an idealized background counter. Here the average wait is two seconds, corresponding to half an event per second. NumPy takes the mean wait, the inverse of the rate.

Most waits are short, but occasional long waits make a tail to the right. The mean isn’t where the tallest bar has to be.

mean_wait_s = 2.0
waiting_times = distribution_rng.exponential(scale=mean_wait_s, size=5000)
waiting_times[:8]

fig_wait, ax_wait = plt.subplots(figsize=(6, 4), constrained_layout=True)
ax_wait.hist(waiting_times, bins=35, color="tab:orange", edgecolor="white")
ax_wait.axvline(mean_wait_s, color="black", linestyle="--", label="Model mean")
ax_wait.set(xlabel="Time between events (s)", ylabel="Number of waits", title="Exponential")
ax_wait.legend()
plt.show()

Binomial: successes in a fixed number of trials#

  • Fix the number of attempts and the success probability.

  • One sample is a count from one batch.

Suppose 20 photons reach a detector and each has a 30% chance of being detected. The binomial distribution models the number detected if these trials are independent and share the same probability. It is also useful for counting heads in a fixed number of coin flips.

Here n is the number of trials in each batch. size=5000 means repeat the whole batch 5000 times. Each result must be an integer from zero through n. The expected count per batch is n * p.

n_trials = 20
detection_probability = 0.3
detected_counts = distribution_rng.binomial(n=n_trials, p=detection_probability, size=5000)
detected_counts[:12]

binomial_edges = np.arange(-0.5, n_trials + 1.5, 1)
fig_detected, ax_detected = plt.subplots(figsize=(6, 4), constrained_layout=True)
ax_detected.hist(detected_counts, bins=binomial_edges, color="tab:green", rwidth=0.85)
ax_detected.set(xlabel=f"Photons detected out of {n_trials}", ylabel="Number of batches",
                title="Binomial", xticks=np.arange(0, n_trials + 1, 2))
plt.show()

Poisson: events in a fixed interval#

  • Fix the observation time, not the number of attempts.

  • lam is the expected count in that interval.

For independent arrivals at a constant average rate, use a Poisson distribution to model how many events a counter records in a fixed time. The exponential example asked how long to wait. This one asks how many events arrive.

An average wait of two seconds gives five expected events in ten seconds: interval divided by mean wait. lam is this expected count, not the rate per second. Each sample below is the count in one interval. Unlike the binomial model, there is no fixed upper limit on the count.

interval_s = 10.0
expected_events = interval_s / mean_wait_s
event_counts = distribution_rng.poisson(lam=expected_events, size=5000)
event_counts[:12]

  • Use integer bins again. Include the largest observed count.

The last edge is half a unit beyond the largest sampled count. This shows every sample; it doesn’t put an upper limit on the Poisson distribution itself.

poisson_edges = np.arange(-0.5, event_counts.max() + 1.5, 1)
fig_events, ax_events = plt.subplots(figsize=(6, 4), constrained_layout=True)
ax_events.hist(event_counts, bins=poisson_edges, color="tab:purple", rwidth=0.85)
ax_events.set(xlabel=f"Events in {interval_s:g} seconds", ylabel="Number of intervals",
              title=f"Poisson: expected count = {expected_events:g}",
              xticks=np.arange(0, event_counts.max() + 1, 2))
plt.show()

Pick a model before writing the sampling call#

  • What does one draw represent?

  • Which values are possible?

If one draw represents…

Start with…

A position sampled evenly along an interval

Uniform

A face of a fair die

Discrete uniform (integers)

An approximate measurement error from many small effects

Normal

A wait between independent arrivals at a constant average rate

Exponential

Successes in a fixed number of independent trials with the same probability

Binomial

An event count over a fixed interval in that constant-rate arrival model

Poisson

Try doubling interval_s and rerun the Poisson draw and histogram. Predict how the center changes before running it. Why doesn’t fixing the observation time give us a fixed maximum count?

Now change the binomial detection probability to 0.8. Rerun its draw and plot. What moves, and what remains the largest possible count? Changing histogram bins alone would not change either model.

Find the reseeding mistake#

  • Where does the generator get created?

  • Predict the list before running it.

This code runs, but it doesn’t produce the five successive draws we intended.

repeated = []
for index in range(5):
    bad_rng = np.random.default_rng(1234)
    repeated.append(bad_rng.random())
repeated

  • Create it once, outside the loop.

Each iteration above goes back to the same starting state and takes the first number again. Move the setup outside the loop so the generator can advance.

good_rng = np.random.default_rng(1234)
draws = []
for index in range(5):
    draws.append(good_rng.random())
draws

Use the numbers in a simulation#

  • Sample points in a square.

  • Count the points inside a quarter circle.

Take a square with \(0 \leq x < 1\) and \(0 \leq y < 1\). The quarter circle of radius one has area \(\pi/4\), and the square has area one. For points sampled uniformly in the square, the probability of landing inside the quarter circle is \(\pi/4\).

That gives us an estimate: four times the fraction inside. We already know \(\pi\), which gives us something to check. This is a Monte Carlo calculation: we use random samples to estimate a quantity.

simulation_seed = 271828
n_points = 10000
simulation_rng = np.random.default_rng(simulation_seed)

  • One row per point. Two coordinates per row.

The tuple (n_points, 2) sets the shape of the array. Keep this sampling call unchanged when we reproduce the run later.

points = simulation_rng.random((n_points, 2))
points[:5]

  • Separate the two columns.

x = points[:, 0]
y = points[:, 1]

  • One True or False per point.

A point is inside when \(x^2 + y^2 \leq 1\).

inside = x**2 + y**2 <= 1
inside[:10]

  • The mean of a mask is the fraction that’s True.

NumPy counts True as one and False as zero in this average.

hit_fraction = np.mean(inside)
pi_estimate = 4 * hit_fraction

  • Look at a subset; calculate with all the points.

We’ll plot the first 600 points to keep the picture readable. The estimate uses all n_points points. ~ flips a Boolean mask.

shown = 600
shown_inside = inside[:shown]
fig_points, ax_points = plt.subplots(figsize=(5, 5), constrained_layout=True)
ax_points.scatter(x[:shown][shown_inside], y[:shown][shown_inside], s=12, label="Inside")
ax_points.scatter(x[:shown][~shown_inside], y[:shown][~shown_inside], s=12, label="Outside")
circle_x = np.linspace(0, 1, 201)
ax_points.plot(circle_x, np.sqrt(1 - circle_x**2), color="black")
ax_points.set(xlabel="x", ylabel="y", xlim=(0, 1), ylim=(0, 1), aspect="equal")
ax_points.legend()
plt.show()

Put an uncertainty on it#

  • A repeatable result can still fluctuate between seeds.

  • Report the estimate and its uncertainty together.

Each point gives a zero or one with probability \(p\) of getting one. For \(N\) independent points, the variance of their average is \(p(1-p)/N\). We don’t know \(p\) in a calculation like this, so we estimate it with our hit fraction \(\hat p\). Multiplying the fraction by four also multiplies its standard error by four:

\[ \widehat{\mathrm{SE}}(\hat\pi) = 4\sqrt{\frac{\hat p(1-\hat p)}{N}}. \]

This is an approximate standard error of the estimate, not the spread of individual points and not a guaranteed bound. It is useful here because we have independent draws and plenty of points both inside and outside the circle.

pi_standard_error = 4 * np.sqrt(hit_fraction * (1 - hit_fraction) / n_points)
print(f"seed={simulation_seed}, N={n_points}")
print(f"pi estimate: {pi_estimate:.4f} +/- {pi_standard_error:.4f} (one standard error)")

  • Compare the discrepancy with the uncertainty.

The difference from np.pi is the discrepancy we observe in this run. The standard error estimates the typical sampling variation. They aren’t the same thing. Even a correct simulation won’t always land within one standard error of the known value.

difference = pi_estimate - np.pi
print(f"Estimate minus pi: {difference:.4f}")
print(f"Difference / standard error: {difference / pi_standard_error:.2f}")

Make the run easy to repeat#

  • Pass the generator into the function.

  • Keep the same sampling call.

Put the calculation into a function so we can change the sample count without copying the code. The function receives a generator and advances it; it doesn’t choose or reset a seed. It returns the estimate and its standard error.

def estimate_pi(n_points, rng):
    points = rng.random((n_points, 2))
    inside = points[:, 0]**2 + points[:, 1]**2 <= 1
    fraction = np.mean(inside)
    estimate = 4 * fraction
    standard_error = 4 * np.sqrt(fraction * (1 - fraction) / n_points)
    return estimate, standard_error

  • Recreate the generator. Repeat the calculation.

Start from the same seed and use the same sample count. If you changed earlier cells, rerun the simulation from its settings cell through the uncertainty cell before comparing.

check_rng = np.random.default_rng(simulation_seed)
check_estimate, check_error = estimate_pi(n_points, check_rng)
print("Same estimate:", check_estimate == pi_estimate)
print("Same uncertainty:", check_error == pi_standard_error)

  • The seed isn’t the whole record.

Adding an extra draw before the simulation advances the state. Changing the order or shape of the sampling calls can change which values get used. Keep the code and inputs as well as the seed. NumPy also doesn’t guarantee identical Generator results across all versions and environments.

The NumPy compatibility policy explains those limits. Matching a previous run checks reproducibility; it doesn’t prove the calculation itself is correct. That’s why we also checked the geometry and compared with a known value.

Change the sample count#

  • Four times the points: about half the standard error.

  • More points don’t guarantee a closer answer in every run.

Run the calculation at three sample counts. We start a fresh generator with the same seed for each count, so the smaller samples are contained in the larger ones here. These are related estimates, not three independent trials.

for count in [1000, 4000, 16000]:
    count_rng = np.random.default_rng(simulation_seed)
    estimate, standard_error = estimate_pi(count, count_rng)
    print(f"seed={simulation_seed}, N={count:5d}: {estimate:.4f} +/- {standard_error:.4f}")

  • Change the seed, keep the sample count.

Change simulation_seed in the settings cell and rerun from there. Does the estimate change? Does its uncertainty change much? Don’t search for a seed that gives the closest answer and report only that run.

The uncertainty should fall roughly as \(1/\sqrt{N}\) as we increase the sample count. The discrepancy from the known value can move either way in one run.

Save what you ran#

  • Seed, inputs, code, and environment.

  • Save the result with its uncertainty.

We’ll write a small JSON file next to the notebook. JSON stores named values in text, so you can open it later and see what produced the result. Python’s json and platform modules come with Python; they don’t need a separate install.

import json
import platform
from pathlib import Path

  • Record the generator, not just the seed.

This records the main simulation run. The generator name identifies the algorithm used underneath NumPy’s sampling methods. Save your edited notebook too: the file below records inputs and results, but it doesn’t contain the code.

run_record = {
    "seed": simulation_seed,
    "n_points": n_points,
    "sampling_call": "rng.random((n_points, 2))",
    "bit_generator": type(simulation_rng.bit_generator).__name__,
    "python_version": platform.python_version(),
    "numpy_version": np.__version__,
    "pi_estimate": float(pi_estimate),
    "standard_error": float(pi_standard_error),
}
print(json.dumps(run_record, indent=2))

  • Write it. Read it back.

Reusing this filename replaces the previous record. Use a different name if you want to keep several runs. In Colab, download the file from the Files panel before leaving the runtime.

record_path = Path("pi_run.json")
record_path.write_text(json.dumps(run_record, indent=2))

saved_record = json.loads(record_path.read_text())
saved_record == run_record

Try it with a die#

  • Use the same pieces for a new calculation.

  • Explain what changes when you repeat it.

Simulate 6000 rolls of a fair six-sided die with a seed you choose.

  1. Use rng.integers and find the fraction of rolls that are six.

  2. Plot a histogram with one bin per face. Label the axes and compare the counts with the expected count per face, 6000 / 6.

  3. Estimate the standard error of the fraction of sixes with \(\sqrt{\hat p(1-\hat p)/N}\). There is no factor of four for this calculation.

  4. Report the seed, number of rolls, fraction, and standard error. Compare the result with \(1/6\).

  5. Recreate the generator with the same seed and confirm you get the same rolls.

  6. Try a new seed, then four times as many rolls. Explain which differences you expect and which you can’t predict for a single run.

Use new names for your variables so this doesn’t replace the circle calculation.

# Create your generator, roll the die, plot a histogram, and analyze the results here.

  • Restart. Run from the top.

  • Check that the saved record matches the run.

Save your edits, restart the kernel or runtime, and run all cells in order. Do the replay comparisons still pass? Can someone else tell which seed, sample count, and code produced your result without asking you?

For other distributions and their arguments, use NumPy’s random generator reference.

For the plotting options used here, see Matplotlib’s histogram reference. NumPy documents the parameters for exponential, binomial, and Poisson sampling.

Contents