Random numbers — the modern way

For thirteen lessons, one line has opened almost every example — rng = np.random.default_rng(seed) — and we used what it produced (fake datasets, weight matrices, attention inputs) without once asking what it is. This final NumPy lesson settles the debt. Random numbers are everywhere in data work: synthetic data, train/test splits, bootstrap resampling, weight initialisation, dropout masks. NumPy’s modern API (1.17+) is faster, safely threadable, and much harder to misuse than the old one — and default_rng is the whole reason every example you have read is reproducible.

Use default_rng, not np.random.seed

import numpy as np

# Old API — global state, harder to reason about
np.random.seed(42)
print("old:", np.random.standard_normal(3))

# New API — explicit Generator object, no global state
rng = np.random.default_rng(42)
print("new:", rng.standard_normal(3))

old: [ 0.49671415 -0.1382643   0.64768854]
new: [ 0.30471708 -1.03998411  0.7504512 ]

NumPy’s random numbers are pseudo-random: a deterministic algorithm (a PRNG, Pseudo-Random Number Generator) produces numbers that look random but are fully determined by an initial value, the seed. Same seed → same sequence, every time. Without a seed, NumPy seeds itself from the OS clock, so you get different numbers each run — but that is not true randomness, just an unpredictable (and irreproducible) starting point. (The two APIs print different numbers because they use different underlying generators — that is expected, not a bug.)

The old np.random.seed(...) sets a hidden module-level state, which any library you import can share — and silently reset. The new default_rng(seed) returns a Generator object you pass around explicitly. No spooky action at a distance.

The distributions you’ll actually use

import numpy as np

rng = np.random.default_rng(0)

# Uniform on [0, 1) — for probability draws
print("uniform:        ", rng.random(5).round(3))

# Uniform on [low, high)
print("uniform 0-100:  ", rng.uniform(0, 100, 5).round(1))

# Integers in [low, high)
print("integers:       ", rng.integers(1, 7, 10))   # 10 die rolls

# Standard normal (mean 0, std 1)
print("standard_normal:", rng.standard_normal(5).round(3))

# Normal with custom mean / std
print("normal(170, 8): ", rng.normal(170, 8, 5).round(1))  # heights in cm

# Sampling from a list
fruits = np.array(["apple", "banana", "cherry", "date"])
print("choice:         ", rng.choice(fruits, 3))

# Weighted choice
weights = [0.7, 0.1, 0.1, 0.1]
print("weighted choice:", rng.choice(fruits, 5, p=weights))

uniform:         [0.637 0.27  0.041 0.017 0.813]
uniform 0-100:   [91.3 60.7 72.9 54.4 93.5]
integers:        [2 5 5 1 3 6 4 1 5 5]
standard_normal: [-0.732 -0.544 -0.316  0.412  1.043]
normal(170, 8):  [169.  180.9 164.7 172.8 177.2]
choice:          ['date' 'banana' 'banana']
weighted choice: ['date' 'apple' 'apple' 'apple' 'apple']

For most ML work you will lean on standard_normal (weight init), integers (index sampling), choice (sampling with/without replacement, with optional weights), and random (drawing probabilities). Note integers(1, 7) rolled values 1–6 — the high bound is exclusive — and the weighted choice, with 0.7 mass on “apple”, returned it four times in five.

Reproducible train / test split

Every model evaluation starts here. Shuffle the indices with a seeded generator and slice:

import numpy as np

rng = np.random.default_rng(42)

# Pretend we have 1000 labeled samples
n = 1000
X = rng.standard_normal((n, 4))
y = (X[:, 0] + X[:, 1] > 0).astype(int)

# Reproducible shuffle
idx = rng.permutation(n)

split = int(0.8 * n)
train_idx, test_idx = idx[:split], idx[split:]

X_train, X_test = X[train_idx], X[test_idx]
y_train, y_test = y[train_idx], y[test_idx]

print("train:", X_train.shape, "test:", X_test.shape)
print("class balance, train:", y_train.mean().round(3))
print("class balance, test: ", y_test.mean().round(3))

train: (800, 4) test: (200, 4)
class balance, train: 0.514
class balance, test:  0.455

Because we seeded rng, this 800/200 split is identical every run (and the class balance is close, ~0.51 vs ~0.46, as a random split should be). Change the seed and you get a different — but still reproducible — split. This is exactly what sklearn.model_selection.train_test_split(random_state=42) does under the hood.

Synthetic data for unit tests

When you need data that looks like the real thing, generate it from a known distribution so your tests have a ground truth to check against:

import numpy as np

rng = np.random.default_rng(0)

# Two clusters in 2D — useful for testing clustering / classifiers
n_per = 200
cluster_a = rng.normal(loc=[0, 0], scale=[1, 1], size=(n_per, 2))
cluster_b = rng.normal(loc=[5, 5], scale=[1, 1], size=(n_per, 2))

X = np.vstack([cluster_a, cluster_b])
y = np.concatenate([np.zeros(n_per), np.ones(n_per)]).astype(int)

# Shuffle them together
order = rng.permutation(len(X))
X, y = X[order], y[order]

print("X:", X.shape, "y:", y.shape)
print("cluster means:")
print(" a:", X[y == 0].mean(axis=0).round(2))
print(" b:", X[y == 1].mean(axis=0).round(2))

X: (400, 2) y: (400,)
cluster means:
 a: [-0.1   0.02]
 b: [4.96 5.01]

The recovered cluster means (≈[0,0] and ≈[5,5]) match the centres we asked for — a built-in sanity check that the generator did what you meant.

Bootstrap — confidence intervals from any statistic

Bootstrapping is “what if I had 1000 different datasets?” — you resample with replacement from your one dataset and recompute the statistic on each resample. The spread tells you how uncertain the estimate is.

import numpy as np

rng = np.random.default_rng(7)

# Customer satisfaction scores from a small survey
scores = np.array([7, 8, 9, 6, 7, 8, 5, 9, 10, 7, 6, 8, 7, 9, 8, 7, 6, 8, 9, 7])
print("observed mean:", scores.mean())

# Bootstrap the mean: 5000 resamples
n_boot = 5000
boot_means = np.empty(n_boot)
for i in range(n_boot):
    sample = rng.choice(scores, size=len(scores), replace=True)
    boot_means[i] = sample.mean()

# 95% confidence interval
lo, hi = np.percentile(boot_means, [2.5, 97.5])
print(f"95% CI: [{lo:.2f}, {hi:.2f}]")

observed mean: 7.55
95% CI: [7.00, 8.10]

The observed mean is 7.55, and the bootstrap puts a 95% interval of roughly [7.00, 8.10] around it — no formula required. That loop is the bread and butter of “how confident are you in this number?” and works for any statistic — median, 90th percentile, regression coefficient, anything (it is the same bootstrap the estimation & confidence intervals lesson built).

Bit generators briefly

default_rng uses PCG64 (Permuted Congruential Generator, 64-bit — a modern PRNG that passes rigorous statistical tests). If you need parallel streams (one per worker, with no overlap), call rng.spawn(n) — each child gets a deterministic, independent stream:

rng = np.random.default_rng(42)
children = rng.spawn(4)              # 4 independent generators
# Hand each to a worker process. Each is reproducible from its own seed.

For most work you will never touch this. But the moment you go multi-process and want reproducibility, spawn is the right tool.

In one breath

Use rng = np.random.default_rng(seed) — an explicit Generator (PCG64) — not the old global np.random.seed, which any imported library can silently reset. Numbers are pseudo-random: fully determined by the seed (same seed → same sequence; no seed → an unpredictable-but-still-deterministic OS-clock start, not true randomness). Daily distributions: random/uniform (probabilities), integers (exclusive high), standard_normal/normal (weight init), choice (sampling ± replacement, optional weights). Seed once and your train/test split (rng.permutation), synthetic data, and bootstrap (choice(..., replace=True)) all become reproducible — the contract behind train_test_split(random_state=…). For parallel workers, rng.spawn(n).

Practice

A question to carry forward

That closes the NumPy section — and look at what the array gave us: speed, vectorisation, broadcasting, the whole compiled-C engine under the data stack. But notice what the array withholds. Every array in fourteen lessons was anonymous: position 0, position 1, position 2, one dtype throughout, no labels, and no notion of a missing value beyond a bare NaN. Real data is not like that. A sales table has named columns, a date index you look rows up by, strings beside floats, and holes where a value never arrived.

The next section is the library that wraps a NumPy array in exactly those things — labels, an index, mixed types, first-class missing data — and it is the most-used tool in all of data science: pandas. Its atom is the Series: a 1-D array that has grown an index, a label on every value. What does that one addition actually buy you — and why is it at once pandas’ most magical feature and its most notorious gotcha?