Sampling Algorithms: A Practical Python Guide

Sampling algorithms either select records from existing data or generate values governed by a target probability distribution. Python has direct tools for common population samples, while streams and custom distributions require algorithms such as reservoir sampling, inverse transform, rejection sampling, or Markov chain Monte Carlo.

A practical sampling workflow first identifies whether the task concerns observed records or a target distribution, then chooses by replacement, weights, stream size, CDF availability, and independence requirements.

What Is a Sampling Algorithm?

The word sampling covers two related but distinct jobs.

Population sampling selects items that already exist. The population might be a list of customers, rows in a dataset, test cases, or events arriving from an iterator. Questions about replacement, category balance, weights, and storage determine the right method.

Distribution sampling generates numeric values governed by a probability law. The law might have a cumulative distribution function, known as a CDF, that can be inverted. In harder cases, you may know only a density up to a constant. That changes which algorithms are available.

Keep these categories stable:

Task	Input	Result	Typical methods
Population sampling	Existing records	Selected records	random.sample, random.choices, stratification, bootstrap, reservoir sampling
Distribution sampling	A probability law	Generated values	Inverse transform, rejection sampling, Metropolis-Hastings, Gibbs sampling
Expectation estimation	Target and proposal densities	A weighted estimate	Importance sampling

Before choosing an algorithm, answer these questions:

• Are you selecting records or generating values?
• Can the same item appear more than once?
• Do items have equal probabilities?
• Is the full population available in memory?
• Can you evaluate and invert the target CDF?
• Can you sample from a convenient proposal distribution?
• Must the returned draws be independent?
• Is the result used in a security-sensitive setting?

That decision path prevents a common category error: using a method that produces valid random output but answers the wrong statistical question.

Sampling from Lists and Arrays

For an in-memory Python sequence, replacement is the first decision.

Without replacement means a population position can be selected at most once. Python’s random.sample(population, k) selects k distinct positions, so equal-valued entries can both appear.

With replacement means each draw starts from the complete population, so an item may appear several times. random.choices(population, weights=None, k=1) uses this rule and can apply relative weights.

This standalone example demonstrates both:

import random

rng = random.Random(17)
records = ["Ada", "Grace", "Linus", "Guido", "Margaret"]

without_replacement = rng.sample(records, k=3)
with_replacement = rng.choices(records, k=8)
weighted = rng.choices(
    records,
    weights=[1, 1, 1, 4, 1],
    k=8,
)

assert len(without_replacement) == 3
assert len(set(without_replacement)) == 3
assert len(with_replacement) == 8
assert len(weighted) == 8

A seed makes the pseudorandom sequence repeatable for testing. It does not make one particular sample more representative.

NumPy’s Generator.choice is useful when the population or result already belongs in a NumPy workflow. It samples with replacement by default. Pass replace=False for sampling without replacement and p for per-entry probabilities.

import numpy as np

rng = np.random.default_rng(17)
records = np.array(["Ada", "Grace", "Linus", "Guido", "Margaret"])

uniform = rng.choice(records, size=3, replace=False)
weighted = rng.choice(
    records,
    size=8,
    replace=True,
    p=[0.1, 0.1, 0.1, 0.6, 0.1],
)

assert uniform.shape == (3,)
assert len(set(uniform.tolist())) == 3
assert weighted.shape == (8,)

Weights must match the intended selection probabilities. Validate their length, signs, finiteness, and total before treating the output as meaningful. A successful function call does not prove that the weights encode the right model.

How does stratified sampling preserve groups?

Uniform sampling can underrepresent a small but important group. Stratified sampling divides observed records into defined groups, called strata, and samples within each group.

The allocation can follow population proportions, or it can deliberately take more records from a small group. The second design is useful for analysis, but its raw sample proportions no longer match the original population.

import random
from collections import defaultdict

def stratified_sample(records, group_key, sizes, rng):
    groups = defaultdict(list)
    for record in records:
        groups[group_key(record)].append(record)

    selected = []
    for group, size in sizes.items():
        selected.extend(rng.sample(groups[group], k=size))
    return selected

employees = [
    {"name": "A", "team": "platform"},
    {"name": "B", "team": "platform"},
    {"name": "C", "team": "platform"},
    {"name": "D", "team": "product"},
    {"name": "E", "team": "product"},
]

sample = stratified_sample(
    employees,
    group_key=lambda row: row["team"],
    sizes={"platform": 2, "product": 1},
    rng=random.Random(17),
)

assert len(sample) == 3
assert sum(row["team"] == "platform" for row in sample) == 2
assert sum(row["team"] == "product" for row in sample) == 1

How does bootstrap resampling work?

A bootstrap sample has the same size as the observed dataset and is drawn with replacement. Some records appear multiple times, while others do not appear at all. Repeating that process lets you examine how a statistic changes across resampled datasets.

import numpy as np

rng = np.random.default_rng(17)
observed = np.array([12.0, 14.0, 15.0, 18.0, 21.0])

bootstrap_samples = rng.choice(
    observed,
    size=(2_000, observed.size),
    replace=True,
)
bootstrap_means = bootstrap_samples.mean(axis=1)

assert bootstrap_samples.shape == (2_000, 5)
assert bootstrap_means.shape == (2_000,)

These methods operate on observed values. They do not replace a model for generating new values outside the data you supplied. For a refresher on the containers feeding these APIs, see Python lists, tuples, sets, and dictionaries.

Reservoir Sampling for Unknown-Length Streams

random.sample needs a finite population that is available when the call begins. A stream can be different: it may arrive through an iterator, and its total length may be unknown.

Reservoir sampling selects up to k items without replacement while storing only the reservoir. Algorithm R starts by keeping the first k items. When item number i arrives after that, it chooses an integer from 0 through i - 1. If that integer names a reservoir position, the new item replaces the value in that position.

Here is a complete standalone implementation:

import random
from collections.abc import Iterable
from typing import TypeVar

T = TypeVar("T")

def reservoir_sample(
    stream: Iterable[T],
    k: int,
    rng: random.Random,
) -> list[T]:
    if k < 0:
        raise ValueError("k must be non-negative")

    reservoir: list[T] = []

    for i, item in enumerate(stream, start=1):
        if i <= k:
            reservoir.append(item)
            continue

        j = rng.randrange(i)
        if j < k:
            reservoir[j] = item

    return reservoir

sample = reservoir_sample(
    (value * value for value in range(100)),
    k=5,
    rng=random.Random(17),
)

assert len(sample) == 5
assert len(set(sample)) == 5
assert all(value in {x * x for x in range(100)} for value in sample)

After processing i items, every one of those items has inclusion probability k / i, provided i is at least k. Existing items survive the next step with exactly the probability needed to preserve that invariant. By the end, every possible group of k input items has equal probability.

The algorithm stores O(k) items. It still reads the complete stream once, so it is memory-efficient rather than magically able to skip unseen input.

A repeated-trial test can catch biased replacement logic:

import random
from collections import Counter

def reservoir_sample(stream, k, rng):
    if k < 0:
        raise ValueError("k must be non-negative")

    reservoir = []
    for i, item in enumerate(stream, start=1):
        if i <= k:
            reservoir.append(item)
        else:
            j = rng.randrange(i)
            if j < k:
                reservoir[j] = item
    return reservoir

rng = random.Random(17)
counts = Counter()
trials = 5_000
population_size = 10
k = 2

for _ in range(trials):
    counts.update(reservoir_sample(range(population_size), k, rng))

expected = trials * k / population_size

assert sum(counts.values()) == trials * k
assert all(abs(counts[item] - expected) < 200 for item in range(population_size))

This tolerance is a test guard, not a proof of uniformity. The equal-inclusion argument explains why the algorithm works; the trial catches implementation mistakes such as choosing from the wrong integer range.

For more context on deterministic pseudorandom sequences, see the pseudo random code generator guide.

Sampling Directly from a Distribution

Distribution sampling starts with a probability law rather than a collection of records. Two direct methods are especially useful because their accepted draws can be independent when their underlying uniform values and proposals are independent.

How does inverse-transform sampling work?

If U is uniform on [0, 1] and F is a target CDF, then applying the inverse CDF gives a value with the target distribution:

X = F⁻¹(U)

For an exponential distribution with positive rate rate, the inverse calculation is:

X = -log(1 - U) / rate

This complete implementation uses log1p(-u), which directly computes log(1 - u):

import numpy as np

def sample_exponential(
    rng: np.random.Generator,
    rate: float,
    size: int,
) -> np.ndarray:
    if rate <= 0:
        raise ValueError("rate must be positive")
    if size < 0:
        raise ValueError("size must be non-negative")

    u = rng.random(size)
    return -np.log1p(-u) / rate

rng = np.random.default_rng(17)
draws = sample_exponential(rng, rate=2.0, size=50_000)

assert draws.shape == (50_000,)
assert np.all(draws >= 0)
assert abs(draws.mean() - 0.5) < 0.03

Inverse transform is attractive when the inverse CDF exists and can be evaluated conveniently. It becomes awkward when inversion is unavailable or costly.

How does rejection sampling work?

Rejection sampling uses four pieces:

• A target density q(x).
• A proposal density p(x) that is easy to sample.
• An envelope constant M such that q(x) <= M p(x) wherever the target is positive.
• An acceptance test U < q(X) / (M p(X)).

The proposal must cover the full target support. If p(x) is zero where q(x) is positive, no finite envelope can fix the gap.

This standalone implementation samples the triangular density q(x) = 2x on [0, 1]. The proposal is uniform on [0, 1], so p(x) = 1, M = 2, and the acceptance probability for proposal x simplifies to x.

import numpy as np

def sample_triangular(
    rng: np.random.Generator,
    size: int,
) -> tuple[np.ndarray, float]:
    if size < 0:
        raise ValueError("size must be non-negative")
    if size == 0:
        return np.empty(0), 0.0

    accepted = []
    proposals = 0

    while len(accepted) < size:
        x = rng.random()
        u = rng.random()
        proposals += 1

        if u < x:
            accepted.append(x)

    return np.asarray(accepted), size / proposals

rng = np.random.default_rng(17)
draws, acceptance_rate = sample_triangular(rng, size=40_000)

assert draws.shape == (40_000,)
assert np.all((0 <= draws) & (draws <= 1))
assert abs(draws.mean() - (2 / 3)) < 0.02
assert 0 < acceptance_rate < 1

The expected number of proposals per accepted draw is M. A loose envelope therefore wastes proposals. In higher dimensions, finding a proposal that covers the target without leaving large low-density regions can become difficult.

Accepted draws from this procedure are independent because each attempt uses a fresh independent proposal and uniform value. That property separates rejection sampling from the correlated sequences produced by MCMC.

When Direct Sampling Is Not Practical

Direct methods are not always available. The inverse CDF may be inaccessible, or rejection sampling may require an impractical envelope. Two different responses are easy to confuse: importance sampling estimates an expectation, while MCMC generates a dependent sequence associated with the target.

What does importance sampling estimate?

Suppose the goal is an expectation under target density q, but samples come from proposal density p. Importance sampling corrects proposal draws with weights based on q(x) / p(x).

The raw draws still follow p. The weighted average estimates a target expectation.

This standalone example estimates the second moment of a standard normal target using draws from a wider normal proposal:

import numpy as np

def normal_density(x, standard_deviation):
    scale = standard_deviation
    return (
        np.exp(-0.5 * (x / scale) ** 2)
        / (np.sqrt(2 * np.pi) * scale)
    )

rng = np.random.default_rng(17)
proposal_sd = 2.0
proposal_draws = rng.normal(0.0, proposal_sd, size=100_000)

target_density = normal_density(proposal_draws, 1.0)
proposal_density = normal_density(proposal_draws, proposal_sd)
weights = target_density / proposal_density

estimated_second_moment = np.sum(
    weights * proposal_draws**2
) / np.sum(weights)

assert abs(estimated_second_moment - 1.0) < 0.08

A proposal that rarely visits regions important to the target can produce highly uneven weights. More draws do not repair a proposal that misses target support.

How does Metropolis-Hastings sampling work?

Metropolis-Hastings builds a Markov chain. Each new state is proposed using the current state, then accepted or rejected using a ratio involving the target and proposal.

With a symmetric random-walk proposal, proposal terms cancel. For an unnormalized standard normal target, the following complete implementation is enough:

import math
import numpy as np

def metropolis_hastings_normal(
    rng: np.random.Generator,
    steps: int,
    proposal_scale: float,
    start: float = 0.0,
) -> tuple[np.ndarray, float]:
    if steps <= 0:
        raise ValueError("steps must be positive")
    if proposal_scale <= 0:
        raise ValueError("proposal_scale must be positive")

    chain = np.empty(steps)
    current = start
    current_log_target = -0.5 * current * current
    accepted = 0

    for i in range(steps):
        proposal = current + rng.normal(0.0, proposal_scale)
        proposal_log_target = -0.5 * proposal * proposal

        log_acceptance_ratio = proposal_log_target - current_log_target

        if log_acceptance_ratio >= 0.0 or rng.random() < math.exp(log_acceptance_ratio):
            current = proposal
            current_log_target = proposal_log_target
            accepted += 1

        chain[i] = current

    return chain, accepted / steps

rng = np.random.default_rng(17)
chain, acceptance_rate = metropolis_hastings_normal(
    rng,
    steps=30_000,
    proposal_scale=1.0,
)

demo_tail = chain[2_000:]

assert abs(demo_tail.mean()) < 0.12
assert abs(demo_tail.var() - 1.0) < 0.15
assert 0 < acceptance_rate < 1

The discarded prefix is a choice for this demonstration, not a universal burn-in rule. Starting point, target shape, proposal, and convergence evidence all matter.

A tiny proposal can accept often while moving slowly. A large proposal can travel farther but reject often. No universal acceptance rate settles that tradeoff for every target.

How does Gibbs sampling differ?

Gibbs sampling updates one variable at a time from its conditional distribution given the others. It is useful when those conditional distributions are available even though drawing the full joint distribution directly is difficult.

For a two-variable normal target with correlation rho, the conditional updates are normal:

import numpy as np

def gibbs_bivariate_normal(
    rng: np.random.Generator,
    steps: int,
    rho: float,
) -> np.ndarray:
    if steps <= 0:
        raise ValueError("steps must be positive")
    if not -1 < rho < 1:
        raise ValueError("rho must be between -1 and 1")

    draws = np.empty((steps, 2))
    x = 0.0
    y = 0.0
    conditional_sd = np.sqrt(1.0 - rho * rho)

    for i in range(steps):
        x = rng.normal(rho * y, conditional_sd)
        y = rng.normal(rho * x, conditional_sd)
        draws[i] = (x, y)

    return draws

rng = np.random.default_rng(17)
draws = gibbs_bivariate_normal(rng, steps=30_000, rho=0.7)
demo_tail = draws[2_000:]

assert np.all(np.abs(demo_tail.mean(axis=0)) < 0.1)
assert abs(np.corrcoef(demo_tail.T)[0, 1] - 0.7) < 0.05

Successive Gibbs and Metropolis-Hastings states are generally correlated. Ten thousand chain states therefore should not be treated as equivalent to ten thousand independent direct draws.

Method	Primary objective	Main requirement	Independent draws?	Main checks
Importance sampling	Estimate expectations	Proposal density and target-to-proposal weights	Proposal draws can be independent, but they follow the proposal	Weight concentration and estimate stability
Metropolis-Hastings	Explore a target distribution	Target ratios and a proposal transition	No	Trace behavior, repeated chains, autocorrelation
Gibbs sampling	Explore a joint target	Sampleable conditional distributions	No	Trace behavior, repeated chains, autocorrelation

How to Choose the Right Sampling Algorithm

Use one decision rule from start to finish.

1. Start with the object being sampled. Existing records belong to population sampling. A mathematical target belongs to distribution sampling. A target expectation may call for importance sampling.
2. For existing records, decide replacement and weighting. Use random.sample for a small uniform sample without replacement. Use Generator.choice when arrays, replacement, or explicit probabilities fit the task.
3. Preserve important groups deliberately. Use stratified sampling when category representation is part of the sampling design.
4. Use reservoir sampling for an unknown-length stream. Algorithm R keeps a uniform sample without loading the stream. For unequal weights without replacement, use a weighted reservoir method such as priority-key sampling.
5. Use inverse transform when the inverse CDF is practical. This gives direct independent draws.
6. Use rejection sampling when a good envelope is available. The proposal must cover the target, and a smaller valid M means fewer rejected proposals.
7. Use importance sampling when the deliverable is an expectation. Do not present its unweighted proposal draws as target samples.
8. Consider Metropolis-Hastings when direct independent sampling is impractical. Choose a proposal that explores the target, then inspect dependence and convergence evidence.
9. Consider Gibbs sampling when full conditional distributions are available. Expect correlated draws and validate the resulting chains.

Constraint	First method to consider
Distinct positions from an in-memory sequence	random.sample
Weighted or repeated array selections	Generator.choice
Required representation from defined groups	Stratified sampling
Equal-weight iterator with unknown final length	Reservoir sampling (Algorithm R)
Unequally weighted iterator without replacement	Weighted reservoir sampling, such as priority-key sampling
Invertible target CDF	Inverse transform
Easy proposal with a tight target envelope	Rejection sampling
Weighted expectation under a target	Importance sampling
Target ratios available, direct draws difficult	Metropolis-Hastings
Conditional distributions are sampleable	Gibbs sampling

What mental model keeps these methods separate?

Picture two workbenches.

On the population bench, the objects already exist. You choose which cards to remove from a box, whether to put a card back, and whether some cards deserve a higher chance.

On the distribution bench, there is no box of finished values. You have a rule describing which values should appear more often. Inverse transform builds values directly from uniform inputs. Rejection sampling filters proposals. MCMC walks through the target and keeps a dependent history.

Importance sampling sits beside the second bench with a calculator. Its usual product is a weighted estimate, not a box of unweighted target draws.

Validate Before You Trust the Sample

A seed makes a failed test repeatable, but validation must examine the property each algorithm promises.

• Finite populations: confirm the result size, uniqueness when sampling without replacement from a population whose values are distinct, membership in the source, and long-run inclusion frequencies.
• Weighted categories: compare empirical frequencies with the supplied probabilities across many trials.
• Continuous distributions: inspect a histogram and compare sample moments or quantiles with values implied by the target.
• Rejection sampling: check support, distribution shape, and the observed acceptance rate. A high rate alone does not prove correctness.
• Importance sampling: inspect the weights and repeat the estimate with independent seeds.
• MCMC: inspect traces, compare repeated chains started from different points, and calculate autocorrelation at several lags.

This helper measures lagged autocorrelation for a one-dimensional sequence:

import numpy as np

def autocorrelation(values, lag):
    values = np.asarray(values, dtype=float)

    if lag <= 0 or lag >= values.size:
        raise ValueError("lag must be between 1 and len(values) - 1")

    centered = values - values.mean()
    denominator = np.dot(centered, centered)

    if denominator == 0:
        raise ValueError("autocorrelation is undefined for constant data")

    return np.dot(centered[:-lag], centered[lag:]) / denominator

chain_autocorrelations = {
    lag: autocorrelation(demo_tail, lag)
    for lag in (1, 5, 20)
}
independent = np.random.default_rng(17).normal(size=20_000)
independent_lag_one = autocorrelation(independent, lag=1)

assert all(-1 <= value <= 1 for value in chain_autocorrelations.values())
assert chain_autocorrelations[1] > 0.5
assert abs(independent_lag_one) < 0.05

The positive chain autocorrelation means the 28,000 retained states contain less information than 28,000 independent draws. Estimating an effective sample size requires autocorrelations across enough lags, not just lag one.

For a quick histogram check without relying on a plot, compare observed bin proportions with probabilities implied by the target:

import numpy as np

rng = np.random.default_rng(17)
draws = rng.random(40_000)

counts, _ = np.histogram(
    draws,
    bins=[0.0, 0.25, 0.5, 0.75, 1.0],
)
proportions = counts / counts.sum()

assert np.all(np.abs(proportions - 0.25) < 0.02)

Tests like these detect obvious coding errors. They do not establish that a complex model is appropriate for the real problem.

For experiments, create a local random.Random(seed) or np.random.default_rng(seed) and pass it into the code that needs it. This avoids hidden dependence on module-level state and makes failures reproducible.

For passwords, reset links, hard-to-guess identifiers, or adversarial selection, do not use random or a predictable seed. Use secrets:

import secrets

alphabet = "ABCDEFGHJKLMNPQRSTUVWXYZ23456789"
identifier = "".join(secrets.choice(alphabet) for _ in range(16))

assert len(identifier) == 16
assert set(identifier) <= set(alphabet)

Common mistakes

• Forgetting that Generator.choice samples with replacement by default.
• Asking for a sample larger than the population when sampling without replacement.
• Supplying weights that describe convenience rather than the intended selection design.
• Calling bootstrap output new observations from the real process.
• Loading a complete iterator before applying reservoir sampling, which defeats its storage advantage.
• Using rejection sampling with a proposal that misses part of the target support.
• Treating importance proposal draws as though they came from the target.
• Reading an MCMC histogram while ignoring autocorrelation and poor movement through the target.
• Choosing a seed after inspecting results, which turns reproducibility into result shopping.
• Using ordinary pseudorandom functions for security-sensitive values.

These mistakes often look reasonable in short scripts. The common Python mistakes guide covers the surrounding issues with mutable data, iterator consumption, and hidden state.

Where the Books Fit

Sampling combines Python collections, iterators, numerical code, probability, and careful testing. The Python learning roadmap places those skills in a practical order. For a structured path from core syntax through generators, testing, data structures, and interview-level algorithms, Python in Three Months provides the longer job-ready route. Keep the decision table nearby, implement each method once, and validate the property it claims before using the sample.