Advanced randomness

Independent streams with `SeedSequence`

When multiple parts of the program (or multiple processes) need random numbers, does not use "consecutive" seeds (1, 2, 3...). The documented way is to start from one a root seed and derive independent streams with spawn. The same root reproduces the same streams, and they do not overlap:

import numpy as np

raiz = np.random.SeedSequence(2026)
geradores = [np.random.default_rng(filho) for filho in raiz.spawn(3)]
amostras = [g.random(2000) for g in geradores]
print(len(amostras), abs(np.corrcoef(amostras)[0, 1]) < 0.1)

a = np.random.default_rng(np.random.SeedSequence(2026).spawn(3)[1]).random(2)
b = np.random.default_rng(np.random.SeedSequence(2026).spawn(3)[1]).random(2)
print(np.array_equal(a, b))

3 True
True

The generator itself also knows how to split, with the method spawn, which returns new independent generators:

filhos = np.random.default_rng(1).spawn(2)
print(type(filhos[0]).__name__)

Generator

In a parallel program, you create the children in the main process and hand one to each worker. This way each has its own stream, and the total result is reproducible, no matter the order in which the workers finish.

Saving and restoring state

The state of the generator can be saved and restored, allowing a simulation to resume from where it left off:

rng = np.random.default_rng(1)
estado = rng.bit_generator.state
primeiro = rng.random(2)
rng.bit_generator.state = estado
print(np.array_equal(primeiro, rng.random(2)))

True

Monte Carlo simulation

Many problems without a closed formula are solved by sampling many times and counting. The classic is estimating π: sample points in a square of side 1, and the fraction that falls inside the quarter circle is worth π/4. The whole calculation is vectorized:

import math


def estimar_pi(n, rng):
    xy = rng.random((n, 2))
    dentro = (xy ** 2).sum(axis=1) <= 1
    return 4 * dentro.mean()


estimativa = estimar_pi(1_000_000, np.random.default_rng(0))
print(abs(estimativa - math.pi) < 0.01)

True

The error decreases with the square root of the number of samples: to gain one decimal place of precision, you need a hundred times more samples. You can see this by repeating the experiment and measuring how much the estimates vary:

def dispersao(n, repeticoes=30):
    rng = np.random.default_rng(3)
    return np.std([estimar_pi(n, rng) for _ in range(repeticoes)])


print("mais amostras, menos variação:", dispersao(100_000) < dispersao(1_000))

mais amostras, menos variação: True

Bootstrap: uncertainty without a formula

How to know how much to trust an average calculated from a sample? The bootstrap resamples the sample itself, with replacement, many times, calculates the average of each resample, and uses their dispersion as a measure of uncertainty. The matrix of resamples is generated at once, and the average is taken along an axis:

rng = np.random.default_rng(5)
amostra = rng.normal(100, 15, 200)
medias = rng.choice(amostra, size=(2000, amostra.size), replace=True).mean(axis=1)
baixo, alto = np.percentile(medias, [2.5, 97.5])
print(baixo < amostra.mean() < alto, alto - baixo < 10)

True True

The interval between the 2.5% and 97.5% percentiles is a 95% confidence interval for the mean. Here the width was below 10, close to what was expected by theory (about 4, since the standard error is 15 / √200).

Sampling from any distribution

If you know the cumulative distribution function, the inverse transform transforms uniform draws into draws from the distribution you want. For the exponential, -ln(1 − u) / λ:

lam = 2.0
x = -np.log(1 - rng.random(100_000)) / lam
print(abs(x.mean() - 1 / lam) < 0.01)

True

Register the generator and the version

For an experiment that needs to be repeated years later, keep the root seed, the type of generator (the default is the PCG64) and the version of NumPy. The same seed reproduces the same numbers in the same version, but the project reserves the right to change the algorithms between versions of the Generator.