Image Author

Simple random sampling is the first technique most people learn, and the one most people keep using long after they should have moved on. It works — that is not the argument against it. The argument is that it is often wasteful, sometimes dangerously imprecise, and poorly matched to most real-world populations.

To reliably detect an event that occurs in 1% of a population, simple random sampling typically requires thousands of observations just to collect enough rare-event cases. A better-matched technique can get you there with a fraction of the data, at lower cost, and with tighter confidence intervals.

Simple random sampling means every member of the population has an equal probability of being selected, drawn independently. It is the right tool when your population is homogeneous, small, and fully accessible. Most real populations are none of those things.

What follows are five techniques that outperform SRS in specific, common situations. Each solves a real problem, and the Python code for each is runnable.

Technique 1: Stratified Sampling

The Problem It Solves

When your population contains meaningfully different subgroups, simple random sampling treats them as interchangeable. A sample of 100 from a user base with 80% free users, 15% Pro, and 5% Enterprise might return only three or four Enterprise responses by pure chance. Drawing conclusions about your highest-value segment from four data points is not viable.

Stratified sampling fixes this by guaranteeing representation from every subgroup you care about.

How It Works

You divide the population into non-overlapping groups called strata, then sample independently from each one. Proportional allocation gives each stratum a share of the total sample size equal to its share of the population. Optimal allocation goes further, assigning more samples to strata with higher internal variance.

The proportional allocation formula: if stratum h has N_h members out of a total population of N, and you want a total sample of n, then stratum h gets n_h = n × (N_h / N) samples.

When to Use It

Use stratified sampling when your population has known subgroups that behave differently from each other and you need accurate estimates for each group individually. Stratification generally produces more precise estimates than SRS when a heterogeneous population is split into internally homogeneous groups.

Real Scenario

A SaaS company runs a satisfaction survey across a user base that is 80% free-tier, 15% Pro, 5% Enterprise, with a sample of 200. Under SRS: roughly 160 free, 30 Pro, and 10 Enterprise responses — barely enough to say anything reliable about Enterprise. Under proportional stratified sampling, allocations are guaranteed, and Enterprise can be intentionally oversampled and weighted back to population proportions for overall estimates.

Python Code

import pandas as pd
import numpy as np

# Simulate a user population of 10,000
np.random.seed(42)

population = pd.DataFrame({
    "user_id": range(10_000),
    "tier": np.random.choice(
        ["Free", "Pro", "Enterprise"],
        size=10_000,
        # Population proportions: 80% Free, 15% Pro, 5% Enterprise
        p=[0.80, 0.15, 0.05]
    ),
    # Each tier has a different satisfaction distribution
    "satisfaction": np.random.normal(loc=7.0, scale=1.5, size=10_000)
})

def stratified_sample(df, strata_col, total_n, random_state=42):
    """
    Draw a proportional stratified sample from a DataFrame.

    Parameters:
        df           : The full population DataFrame
        strata_col   : Column name that defines the strata
        total_n      : Total desired sample size across all strata
        random_state : Seed for reproducibility

    Returns:
        A DataFrame containing the combined stratified sample
    """
    population_size = len(df)
    sampled_groups = []

    for stratum, group in df.groupby(strata_col):

        # Proportional allocation: each stratum gets a share of total_n
        # equal to its share of the full population
        stratum_proportion = len(group) / population_size
        stratum_n = max(1, round(total_n * stratum_proportion))

        # Sample independently from this stratum
        stratum_sample = group.sample(n=stratum_n, random_state=random_state)
        sampled_groups.append(stratum_sample)

        print(f"  {stratum:12s}: population={len(group):,} | sampled n={stratum_n}")

    return pd.concat(sampled_groups).reset_index(drop=True)


print("Stratified sample allocation (n=200):\n")
sample = stratified_sample(population, strata_col="tier", total_n=200)

print(f"\nTotal sample size: {len(sample)}")

# Compare estimates: SRS vs Stratified
srs = population.sample(n=200, random_state=42)
print(f"\nSRS mean satisfaction:        {srs['satisfaction'].mean():.4f}")
print(f"Stratified mean satisfaction: {sample['satisfaction'].mean():.4f}")
print(f"True population mean:         {population['satisfaction'].mean():.4f}")

Technique 2: Cluster Sampling

The Problem It Solves

Stratified sampling is precise but requires reaching individuals spread across the entire population. If you need to survey patients across 500 hospitals or students across 2,000 schools, travelling to every location is not realistic.

Cluster sampling trades some statistical precision for a major reduction in logistical cost. Instead of sampling individuals across the whole population, you sample entire natural groups and only visit the selected ones.

How It Works

You divide the population into clusters — natural groupings such as schools, hospitals, or zip codes — then randomly select some clusters and collect data from everyone inside them. This is single-stage cluster sampling. In two-stage cluster sampling, you draw a random sample within each selected cluster rather than surveying all members, adding a second layer of efficiency.

The key distinction from stratified sampling: stratified sampling works best when strata are internally homogeneous and different from each other. Cluster sampling works best when clusters are internally heterogeneous and similar to each other — each cluster should be a microcosm of the full population. You sample from all strata but only from some clusters.

When to Use It

Use cluster sampling when your population is geographically or logistically spread out, clusters are naturally defined and roughly equivalent to each other, and the cost of reaching individuals across the whole population is the binding constraint.

Real Scenario

A public health agency wants to measure vaccination rates across a country of 50,000 villages. They randomly select 200 villages and survey every household in those villages. If the villages are demographically similar, the estimates will be valid for the whole country at a fraction of the cost.

Python Code

import pandas as pd
import numpy as np

np.random.seed(42)

# Simulate a population of students spread across 50 schools (clusters)
schools = []
for school_id in range(1, 51):
    school_size = np.random.randint(80, 121)  # Each school has 80-120 students
    school_data = pd.DataFrame({
        "student_id": range(school_size),
        "school_id": school_id,
        # Students within a school have similar scores, but schools differ
        "test_score": np.random.normal(
            loc=np.random.uniform(60, 90),  # Each school has its own mean
            scale=10,
            size=school_size
        )
    })
    schools.append(school_data)

population = pd.concat(schools, ignore_index=True)
print(f"Total population: {len(population):,} students across 50 schools\n")


def two_stage_cluster_sample(df, cluster_col, n_clusters, sample_per_cluster, random_state=42):
    """
    Two-stage cluster sampling.

    Stage 1: Randomly select n_clusters from all available clusters.
    Stage 2: Within each selected cluster, randomly sample sample_per_cluster individuals.

    Parameters:
        df                 : Full population DataFrame
        cluster_col        : Column that identifies cluster membership
        n_clusters         : How many clusters to select in Stage 1
        sample_per_cluster : How many individuals to sample per cluster in Stage 2
        random_state       : Seed for reproducibility
    """
    rng = np.random.default_rng(random_state)

    # Stage 1: randomly select clusters from the full cluster list
    all_clusters = df[cluster_col].unique()
    selected = rng.choice(all_clusters, size=n_clusters, replace=False)
    print(f"Stage 1 -- Selected {n_clusters} clusters from {len(all_clusters)} total")

    # Stage 2: within each selected cluster, sample individuals
    sampled = []
    for cluster_id in selected:
        cluster_data = df[df[cluster_col] == cluster_id]

        # If a cluster is smaller than sample_per_cluster, take the whole cluster
        n = min(sample_per_cluster, len(cluster_data))
        cluster_sample = cluster_data.sample(n=n, random_state=random_state)
        sampled.append(cluster_sample)

    return pd.concat(sampled).reset_index(drop=True)


# Select 10 schools, draw 15 students from each = 150 total
sample = two_stage_cluster_sample(
    population,
    cluster_col="school_id",
    n_clusters=10,
    sample_per_cluster=15
)

print(f"Stage 2 -- Sampled 15 students from each selected school")
print(f"\nTotal sample size:               {len(sample):,}")
print(f"Cluster sample mean test score:  {sample['test_score'].mean():.2f}")
print(f"True population mean:            {population['test_score'].mean():.2f}")

What this code does: Stage 1 uses np.random.default_rng().choice() to select 10 schools at random from 50. Stage 2 loops through each selected school and draws 15 students. The min(sample_per_cluster, len(cluster_data)) guard handles small clusters that cannot satisfy the per-cluster quota. The final sample of 150 students from 10 schools never requires contact with the other 40 schools.

Technique 3: Systematic Sampling

The Problem It Solves

You have a large ordered list — server logs, production line outputs, transaction records — and you want an evenly distributed sample without running a full random number pass across millions of records. Simple random sampling requires assigning a random number to every element first. Systematic sampling skips that entirely.

How It Works

Pick a random starting point within the first k positions, then select every kth element after that. The sampling interval k is the population size divided by the desired sample size. For 100,000 records and a desired sample of 1,000, k = 100. You pick a random start between 0 and 99, then collect elements at positions start, start+100, start+200, and so on.

When to Use It

Use systematic sampling when data arrives sequentially or is already ordered, you want even coverage across the full list, and speed matters. Quality control, log auditing, and production line sampling are the natural homes for this technique.

One limitation: if the ordered list has a hidden periodicity that aligns with your interval k, systematic sampling introduces bias. A classic example is shift-based manufacturing data where every 8th record is a shift-change entry and k happens to be 8. Before applying systematic sampling, check whether the ordering has any rhythmic structure.

Real Scenario

A backend team stores 1 million API request logs per day and wants a 1% sample of 10,000 records without loading the full dataset into memory. Systematic sampling with k = 100 reads the log file in a single linear pass, selecting every 100th entry. The sample is evenly distributed across the day’s activity, and memory usage stays flat throughout.

Python Code

import random

def systematic_sample(population, sample_size):
    """
    Draw a systematic sample from any ordered sequence.

    Parameters:
        population  : A list or sequence of elements to sample from
        sample_size : The desired number of elements in the output sample

    Returns:
        A list of sample_size elements drawn at regular intervals

    How it works:
        1. Calculate k = len(population) // sample_size  (the interval)
        2. Pick a random start position between 0 and k-1
        3. Select elements at positions: start, start+k, start+2k, ...
    """
    N = len(population)

    if sample_size > N:
        raise ValueError(f"sample_size ({sample_size}) cannot exceed population size ({N})")

    # Sampling interval: step size between selected elements
    k = N // sample_size

    # Random start within the first interval -- this is what makes it random
    # Without this, you always pick from the same fixed positions
    start = random.randint(0, k - 1)

    # Collect every kth element starting from the random start
    sample = [population[i] for i in range(start, N, k)]

    # Trim to exact size -- integer division can produce one extra element
    return sample[:sample_size]


# Example: sample 1,000 records from 100,000 API logs in one linear pass
api_logs = [f"request_{i:06d}" for i in range(100_000)]

sample = systematic_sample(api_logs, sample_size=1_000)

k = len(api_logs) // 1_000
print(f"Population size:        {len(api_logs):,}")
print(f"Sample size:            {len(sample):,}")
print(f"Sampling interval k:    every {k}th record")
print(f"First 5 records:  {sample[:5]}")
print(f"Last 5 records:   {sample[-5:]}")

Technique 4: Reservoir Sampling

The Problem It Solves

The previous three techniques all assume you know the population size upfront. Reservoir sampling is for the situation where you do not — and where the data may be too large to hold in memory even if you did. Streaming data arrives continuously, log files grow without a declared end, and database exports can contain tens of millions of rows that cannot be fully loaded at once. Simple random sampling is structurally impossible here because you cannot assign probabilities without knowing n.

How It Works

The algorithm (Algorithm R, introduced by Jeffrey Vitter in 1985) maintains a reservoir of exactly k items. Fill it with the first k elements you see. Then, for each subsequent element at position i, generate a random integer j between 0 and i inclusive. If j is less than k, replace the jth element in the reservoir with the new element; otherwise, discard it.

At any point during the stream, the reservoir contains a perfectly uniform random sample of everything seen so far — every element has a k/n probability of being in the reservoir, where n is the number of elements processed. No knowledge of the total stream size is ever required.

Worth noting: Python’s built-in random.sample() uses reservoir sampling internally when called on an iterator of unknown size.

When to Use It

Use reservoir sampling for data streams, very large files that cannot be fully loaded into memory, or any dataset where the total size is unknown or unbounded. It enables representative sampling for model training, anomaly detection, and statistical analysis without processing the entire dataset.

Real Scenario

A data platform processes 50 million user events per day and wants a live 10,000-record representative sample for real-time dashboarding. Loading 50 million rows into memory is not viable. Reservoir sampling maintains the sample with a single pass through the stream and constant memory usage regardless of how many events have arrived.

Python Code

import random
def reservoir_sample(stream, k):
    """
    Algorithm R: reservoir sampling for streams of unknown or very large size.

    Guarantees a uniformly random sample of exactly k elements
    from any iterable, using O(k) memory regardless of stream length.

    Parameters:
        stream : Any iterable (list, generator, file object, API stream...)
        k      : Desired sample size

    Returns:
        A list of k elements sampled uniformly at random from the stream
    """
    reservoir = []

    for i, element in enumerate(stream):

        if i < k:
            # Phase 1: fill the reservoir with the first k elements
            # No decision needed -- just append everything
            reservoir.append(element)

        else:
            # Phase 2: for each new element, randomly decide whether
            # it replaces something already in the reservoir.
            #
            # Probability of replacement = k / (i + 1)
            # This ensures every element seen so far has equal probability
            # of being in the reservoir at any given moment.
            j = random.randint(0, i)   # Random index from 0 to i inclusive

            if j < k: # This element replaces the one at position j reservoir[j] = element # If j >= k, discard the new element -- no action needed

    return reservoir


# Example 1: sample from a large known list
print("=== Example 1: Large list ===")
data = list(range(1_000_000))
sample = reservoir_sample(data, k=1_000)

print(f"Stream size:  {len(data):,}")
print(f"Sample size:  {len(sample):,}")
print(f"Sample mean:  {sum(sample)/len(sample):,.1f}  (true mean: 499,999.5)")


# Example 2: sample from a generator where the total size is never known
def event_stream():
    """
    Simulates a live data stream of unknown length.
    Could be Kafka events, file lines, API responses, sensor readings, etc.
    The reservoir sampler never knows how many elements will arrive.
    """
    i = 0
    while True:  # Stream runs indefinitely
        yield {"event_id": i, "value": random.gauss(100, 15)}
        i += 1

print("\n=== Example 2: Unknown-length stream (stops at 1M events) ===")

stream = event_stream()
events_to_process = 1_000_000

# Process 1M events from the infinite stream
bounded_stream = (next(stream) for _ in range(events_to_process))
stream_sample = reservoir_sample(bounded_stream, k=10_000)

print(f"Events processed:  {events_to_process:,}")
print(f"Reservoir size:    {len(stream_sample):,}")
print(f"Sample mean value: {sum(e['value'] for e in stream_sample)/len(stream_sample):.2f}")
print(f"Expected mean:     100.00")

What this code does: Phase 1 fills the reservoir without any decision logic. Phase 2 begins once the reservoir is full: for each new element at index i, random.randint(0, i) generates a number between 0 and i inclusive, and the probability it falls below k is exactly k/(i+1) — the probability needed to keep every seen element equally likely to remain in the reservoir. The generator example is the key demonstration: the sampler processes events from an infinite stream, never knowing the total count, while holding exactly 10,000 records in memory throughout.

Technique 5: Importance Sampling

The Problem It Solves

Some of the most consequential questions in data science involve rare events — fraud, equipment failure, extreme medical outcomes, financial tail risk. If fraud occurs in 0.1% of transactions and you draw a random sample of 10,000, you get roughly 10 fraud cases — far too few to estimate tail risk reliably. Importance sampling solves this by deliberately breaking the uniform sampling assumption and then correcting for it mathematically.

How It Works

Instead of sampling uniformly, you sample from a proposal distribution that over-represents the region you care about, then assign each observation an importance weight — the ratio of its likelihood under the true distribution versus the proposal. These weights correct for the biased sampling and restore the validity of your estimates.

In plain terms, importance sampling is deliberately over-fishing in the rare-event pond, then applying a discount to your catch to account for having done so.

When to Use It

Use importance sampling in risk modelling, fraud detection, rare disease research, Monte Carlo integration over heavy-tailed distributions, and anywhere the events you most need to understand are the ones SRS will least often find.

Real Scenario

A fintech company wants to estimate the probability that a transaction loss exceeds $50,000 — rare but catastrophic. Under their historical distribution, such losses occur roughly 0.05% of the time. A simple random sample of 100,000 transactions yields around 50 such events — a noisy estimate. Importance sampling over-samples the high-loss region and reweights, producing a far tighter estimate with the same sample size.

Python Code

import numpy as np

np.random.seed(42)

# Setup: transaction losses follow a log-normal distribution
# Most losses are small. We want to estimate P(loss > THRESHOLD).
TRUE_MEAN = 3.0       # Log-normal parameters (log scale)
TRUE_STD  = 1.0
THRESHOLD = 50.0      # Rare high-loss event we want to estimate probability for
N_SAMPLES = 100_000


# --- Method 1: Simple Random Sampling ---
srs_samples = np.random.lognormal(mean=TRUE_MEAN, sigma=TRUE_STD, size=N_SAMPLES)

# Count how many samples exceed the threshold
srs_estimate = np.mean(srs_samples > THRESHOLD)
srs_variance = np.var(srs_samples > THRESHOLD) / N_SAMPLES

print("=== Simple Random Sampling ===")
print(f"Total samples:        {N_SAMPLES:,}")
print(f"Samples > {THRESHOLD}:     {int(np.sum(srs_samples > THRESHOLD)):,}")
print(f"P(loss > {THRESHOLD}):     {srs_estimate:.6f}")
print(f"Estimate variance:    {srs_variance:.2e}")


# --- Method 2: Importance Sampling ---
# Proposal distribution: shift the mean upward so the tail region
# (where losses > THRESHOLD) is sampled much more frequently
PROPOSAL_MEAN = 5.0   # Higher mean = more samples in the tail
PROPOSAL_STD  = 1.0

# Draw samples from the proposal (biased) distribution
proposal_samples = np.random.lognormal(
    mean=PROPOSAL_MEAN, sigma=PROPOSAL_STD, size=N_SAMPLES
)

def lognormal_log_pdf(x, mean, std):
    """
    Compute the log of the log-normal PDF at values x.
    Using log-space arithmetic avoids numerical underflow
    when dealing with very small probabilities.
    """
    log_x = np.log(x)
    return (
        -np.log(x * std * np.sqrt(2 * np.pi))
        - ((log_x - mean) ** 2) / (2 * std ** 2)
    )

# Importance weights = true_pdf(x) / proposal_pdf(x)
# Computed in log space to avoid underflow, then exponentiated
log_weights = (
    lognormal_log_pdf(proposal_samples, TRUE_MEAN, TRUE_STD)
    - lognormal_log_pdf(proposal_samples, PROPOSAL_MEAN, PROPOSAL_STD)
)
weights = np.exp(log_weights)

# Normalize weights so they sum to 1
normalized_weights = weights / weights.sum()

# IS estimate: weighted average of (sample > threshold)
is_estimate = np.sum(normalized_weights * (proposal_samples > THRESHOLD))

# Effective sample size: measures how many "equivalent" SRS samples
# our weighted IS sample is worth. Higher = better.
eff_sample_size = 1.0 / np.sum(normalized_weights ** 2)
is_variance = is_estimate * (1 - is_estimate) / eff_sample_size

print("\n=== Importance Sampling ===")
print(f"Total samples:          {N_SAMPLES:,}")
print(f"Samples > {THRESHOLD}:       {int(np.sum(proposal_samples > THRESHOLD)):,}")
print(f"Effective sample size:  {eff_sample_size:,.0f}")
print(f"P(loss > {THRESHOLD}):       {is_estimate:.6f}")
print(f"Estimate variance:      {is_variance:.2e}")

print(f"\nVariance reduction: {srs_variance / is_variance:.1f}x tighter than SRS")
print(f"(achieved with the same {N_SAMPLES:,} samples)")

What this code does: The SRS block draws from the true loss distribution and counts samples that clear the threshold — very few, producing a noisy estimate. The importance sampling block shifts the proposal distribution upward so rare high-loss events appear frequently, then computes log-ratio weights to mathematically undo that shift. Log-space arithmetic is important here, as directly computing probability ratios for very small probabilities can produce floating-point underflow that silently corrupts results. The effective sample size at the end quantifies how many SRS-equivalent samples the weighted IS sample is worth.

Conclusion

Simple random sampling is not broken. For a homogeneous population where every subgroup is equally interesting, it is the right call. The problem is that real-world populations have structure, datasets have size constraints, and events have unequal importance.

The five techniques in this article are the standard tools statisticians and data scientists reach for when SRS underperforms. Stratified sampling appears in every serious survey workflow. Cluster sampling underpins every national study. Systematic sampling runs inside every quality control system. Reservoir sampling sits at the core of every streaming data pipeline. Importance sampling drives every serious risk model. The Python implementations above are production-ready starting points.