How do I generate random samples from a distribution in NumPy?

The short version

Create a generator with np.random.default_rng(seed) and call the distribution method you need: normal, uniform, integers, choice. Always pass a seed so the results are reproducible, and use choice with a p argument when the options have different probabilities.

The common distributions

import numpy as np

rng = np.random.default_rng(0)

rng.normal(loc=0, scale=1, size=3)      # normal, given mean and standard deviation
rng.uniform(low=0, high=10, size=3)     # uniform between two bounds
rng.integers(low=0, high=10, size=3)    # whole numbers
rng.choice(["a", "b", "c"], size=3)     # sampling from a list

Each takes a size argument, so you generate a whole array at once rather than looping, which is both faster and simpler.

Weighted sampling

choice accepts probabilities, which is how you simulate options that are not equally likely.

rng.choice(["pass", "fail"], size=100, p=[0.9, 0.1])

The probabilities must sum to one. Without p, every option is equally likely, which is rarely what real data looks like.

Sampling with and without replacement

rng.choice(range(10), size=3, replace=False)   # each item at most once
rng.choice(range(10), size=3, replace=True)    # repeats allowed (the default)

Sampling without replacement models drawing from a finite pool; with replacement models independent draws, and is what bootstrap resampling uses.

Seed it, always

Without a seed, every run gives different numbers, so a simulation result cannot be reproduced or checked. Seeding the generator fixes the whole stream.

rng = np.random.default_rng(42)     # reproducible from here on

Use a local generator rather than the global np.random.seed, so that other libraries cannot disturb your stream and yours cannot disturb theirs.

Do not confuse the parameters

A frequent error is passing the variance where NumPy expects the standard deviation. normal(scale=...) takes the standard deviation, not the variance. Passing a variance of 4 where you meant a standard deviation of 2 gives you a distribution twice as wide as intended, and nothing warns you.

Similarly, integers excludes the upper bound by default (like range), whereas some other libraries include it. Check the documentation for any distribution you have not used before, because these off-by-one and parameterisation differences are silent.

Simulation needs many runs

A single random draw tells you almost nothing. The point of simulation is to repeat it many times and look at the distribution of outcomes.

results = [simulate(rng) for _ in range(10_000)]
np.percentile(results, [5, 50, 95])

Reporting a median with a plausible range says far more than a single run, and it is the whole reason for simulating rather than calculating a point estimate.