NumPy Random Generator

In one line: use np.random.default_rng(seed) and pass the generator around — the legacy np.random.seed global API is deprecated in practice and unsafe in parallel code.

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).

The old API still works and will not be removed — NumPy guarantees its stream for backwards compatibility, which is precisely why it cannot be improved. Do not use it in new code.

The naming differences

Legacy	Modern	Note
np.random.rand(3)	rng.random(3)
np.random.randn(3)	rng.standard_normal(3)
np.random.randint(0, 10)	rng.integers(0, 10)	endpoint=True for inclusive
np.random.choice	rng.choice	modern version is faster without replacement
np.random.permutation	rng.permutation

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

rng = np.random.default_rng(42)
a = rng.random(5)

rng = np.random.default_rng(42)
b = rng.random(5)
np.array_equal(a, b)          # True

Seed once, at the top, and pass rng down. Do not re-seed inside functions — a function that re-seeds destroys the caller's stream and makes independent draws identical.

def bootstrap(data, n, rng):          # ← take the rng as a parameter
    idx = rng.integers(0, len(data), size=(n, len(data)))
    return data[idx].mean(axis=1)

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

# read counts: negative binomial (overdispersed Poisson) — the RNA-seq model
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 — the canonical use

def permutation_test(group_a, group_b, n_perm=10_000, rng=None):
    rng = rng or np.random.default_rng()
    observed = group_a.mean() - group_b.mean()
    pooled = np.concatenate([group_a, group_b])
    n_a = len(group_a)

    null = np.empty(n_perm)
    for i in range(n_perm):
        p = rng.permutation(pooled)
        null[i] = p[:n_a].mean() - p[n_a:].mean()

    # +1 in numerator and denominator: never report p = 0
    return (np.sum(np.abs(null) >= np.abs(observed)) + 1) / (n_perm + 1)

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