NumPy Random Generator

Core Principle: Pseudo-random number generation in NumPy should utilize explicit np.random.default_rng(seed) generator instances rather than deprecated global state functions.

The modern API

rng = np.random.default_rng(42)

rng.random(10)                       # uniform [0, 1)
rng.normal(loc=0, scale=1, size=(3, 4))
rng.integers(0, 100, size=10)        # note: high is EXCLUSIVE by default
rng.choice(genes, size=100, replace=False)
rng.permutation(labels)              # returns a shuffled copy
rng.shuffle(arr)                     # shuffles in place
rng.binomial(n=100, p=0.5, size=10)
rng.poisson(lam=5, size=10)
rng.negative_binomial(n=5, p=0.3, size=10)
rng.multinomial(100, [0.25]*4)

Why not np.random.seed

The legacy interface (np.random.seed, np.random.rand, np.random.randn, np.random.randint) uses a single hidden global state. That causes three real problems:

1. Non-local coupling. Any library you import can consume from the same stream, so adding an unrelated dependency changes your results.
2. No parallel safety. Two threads or processes drawing from one global state either race or produce correlated streams.
3. No reproducible independence. You cannot say "give me 10 independent streams" cleanly.

The Generator API fixes all three, and is also faster and statistically better (PCG64 rather than Mersenne Twister).

While the legacy np.random.seed API is retained for backwards compatibility, new workflows should exclusively instantiate localized Generator objects.

The naming differences

rng.integers excludes the high value by default, matching range(). The legacy randint did too, but random_integers did not, which is part of why the old API was a mess.

Reproducibility

Initialize the generator once with an explicit seed and pass the rng instance to subordinate functions, avoiding internal re-seeding that compromises stochastic independence.

Record the seed in your output. A figure that cannot be regenerated is not reproducible, and "I used a random seed somewhere" is not a record.

Independent parallel streams

ss = np.random.SeedSequence(42)
children = ss.spawn(8)                              # 8 guaranteed-independent seeds
rngs = [np.random.default_rng(c) for c in children]

# in each worker:
result = worker(data, rngs[i])

spawn is the correct way to seed parallel workers. Seeding worker i with seed + i is a common hack that can produce correlated streams and is not safe.

Distributions you will actually use in bioinformatics

# Overdispersed count model (Negative Binomial distribution for RNA-seq read counts)
counts = rng.negative_binomial(n=10, p=0.3, size=(n_genes, n_samples))

# simulate coverage
rng.poisson(lam=30, size=chrom_len)

# allele counts at a het site
rng.binomial(n=depth, p=0.5)

# simulate a random sequence with given base composition
rng.choice(list("ACGT"), size=1000, p=[0.3, 0.2, 0.2, 0.3])

# random genomic positions
rng.integers(0, chrom_len, size=1000)

Negative binomial rather than Poisson for RNA-seq counts is the whole basis of DESeq2/edgeR: biological replicates are overdispersed relative to Poisson counting noise.

Permutation Testing Implementation

The +1 correction matters. With 10,000 permutations the smallest reportable p-value is 1/10,001, not 0.

A fully vectorised version avoids the loop:

idx = np.argsort(rng.random((n_perm, len(pooled))), axis=1)
perms = pooled[idx]
null = perms[:, :n_a].mean(axis=1) - perms[:, n_a:].mean(axis=1)

Bootstrap resampling

n = len(data)
boot_idx = rng.integers(0, n, size=(n_boot, n))       # WITH replacement
boot_means = data[boot_idx].mean(axis=1)
ci = np.percentile(boot_means, [2.5, 97.5])

This is what Statistical Estimation in seaborn does internally for its error bars.

Common mistakes

• Using np.random.seed in new code.
• Re-seeding inside a function, collapsing independent draws to identical ones.
• seed + i for parallel workers. Use SeedSequence.spawn.
• Not recording the seed in the output or the filename.
• rng.choice with replace=False on a huge population: it can be slow; use rng.permutation(n)[:k] for large n.
• Reporting p = 0 from a permutation test. Use the (count + 1) / (n + 1) estimator.
• Assuming rng.integers(0, 10) includes 10. It does not, unless endpoint=True.

See also

Array Creation · Fancy Indexing · Vectorization · Statistical Estimation in seaborn · Applied - Reproducible Environment

Exercise

Write draw_report(seed, n, k), which exercises a seeded Generator and checks the properties a shuffle and a sample are supposed to have. Build every generator with np.random.default_rng(seed), taking the seed from the input, so the run is reproducible. Return ["bad k"] when k is below one or above n.

Otherwise return six lines: valid_permutation=<bool>, that rng.permutation(n) is a rearrangement of 0..n-1; reproducible=<bool>, that a freshly seeded generator repeats it exactly; other_seed_differs=<bool>, that seed + 1 does not; unique_draws=<bool> in_range=<bool> for rng.choice(n, size=k, replace=False); multiset_preserved=<bool>, that an in-place shuffle keeps the same elements; and stream_advances=<bool>, that a second permutation drawn from the same generator differs from the first. The report asserts properties rather than specific draws, because a Generator's bit stream is not guaranteed identical across NumPy versions. For n = 1 the last two comparisons have only one possible outcome, so do not assume they are true.