Bootstrap and permutation tests: confidence intervals and p-values without a formula

You have two versions of an API endpoint. The new version appears faster, but its latency distribution is skewed, a handful of requests take much longer than the rest, and the median matters more to your users than the average.

You calculate the difference between the sample medians. The result looks promising. But two important questions remain:

• How uncertain is the observed difference? You need a confidence interval.

• Could a difference this large arise from random variation if the two versions were equivalent? You need a p-value.

For a difference of means, classical statistical methods often provide convenient formulas. For a difference of medians, an arbitrary percentile, or a custom business metric, getting reliable uncertainty estimates can be considerably harder.

This is where bootstrap confidence intervals and permutation tests become useful.

Both methods use repeated computation to approximate distributions that would otherwise be difficult to derive analytically. Instead of requiring a convenient mathematical formula for every statistic, they repeatedly recalculate the statistic on carefully constructed versions of the observed data.

The crucial word is carefully. Bootstrap resampling and permutation testing answer different questions, rely on different assumptions, and are not interchangeable.

In this tutorial, you’ll implement both from scratch, reproduce them with NumPy and SciPy, interpret their results, and learn when the methods can fail.

1. Two questions, two different resampling methods

Imagine that you collect request latency measurements from two independent groups:

• Baseline: requests served by the existing endpoint.

• Candidate: requests served by a redesigned endpoint.

Your statistic is:

Observed effect = median(candidate) − median(baseline)

A negative value means the candidate has lower median latency.

The same approach applies to many other statistics:

• Difference in median order values between two checkout experiences.

• Difference in 90th-percentile response times.

• Difference in average revenue per transaction.

• A trimmed-mean difference that reduces sensitivity to outliers.

• A custom reliability or performance score.

The bootstrap and permutation test use the statistic differently.

Bootstrap: approximate repeated sampling

Suppose you could repeatedly run your experiment, collecting a new baseline sample and a new candidate sample each time.

Every experiment would produce a slightly different difference of medians.

The distribution of those estimates is called the sampling distribution. It tells you how much your estimate varies because of random sampling.

In practice, you typically have only one observed dataset.

The bootstrap approximates repeated sampling by treating the empirical distribution of your sample as a stand-in for the unknown population distribution.

It repeatedly draws samples with replacement, preserving the original sample size, and recomputes the statistic.

Permutation test: construct a null world

A permutation test asks what your statistic would look like if the group labels carried no information under the null hypothesis.

For two independent groups, the conventional permutation procedure pools all observations and redistributes them between groups without replacement, preserving the original group sizes.

Each redistribution produces a new value of the statistic.

Comparing your observed statistic with this null distribution gives a p-value.

Importantly, the ordinary independent-sample permutation test assumes exchangeability under the null. Equal population medians alone do not guarantee this assumption.

2. Set up a reproducible Python environment

We’ll use NumPy for numerical operations and SciPy for the production-ready statistical implementations.

Create a virtual environment and install the dependencies:

python -m venv .venv
source .venv/bin/activate
python -m pip install "numpy>=2.0" "scipy>=1.15"

The examples use NumPy’s modern random number generator and SciPy’s current random-generator interface.

We’ll set explicit random seeds so results are reproducible, while keeping independent random streams for data generation, bootstrapping, and permutation testing.

The official SciPy bootstrap documentation and permutation-test documentation describe these interfaces, including vectorization, resampling counts, and memory controls.

3. Build a realistic synthetic latency experiment

Rather than downloading a dataset, we’ll generate request latencies from a gamma distribution.

This is a reasonable illustrative model for positive, right-skewed measurements: most requests complete within a moderate range, while a smaller number experience unusually high latency.

We’ll create:

• 80 baseline observations.

• 95 candidate observations.

• A common gamma-distributed source of random latency variation.

• A 50 ms baseline offset and a 15 ms candidate offset.

The population distributions therefore differ by a constant shift of −35 ms. Because adding a constant shifts the median by the same amount, the true population difference in medians is exactly −35 ms.

That known ground truth will help us check whether our estimate and confidence interval are sensible.

The observed difference will not necessarily equal −35 ms because our samples are finite.

4. A complete bootstrap CI and permutation p-value with SciPy

Here is a self-contained example you can save as a Python script and run directly.

import numpy as np
from scipy.stats import bootstrap, permutation_test


# Generate two independent latency samples.
data_rng = np.random.default_rng(2026)

baseline = 50.0 + data_rng.gamma(
    shape=5.0,
    scale=24.0,
    size=80,
)

candidate = 15.0 + data_rng.gamma(
    shape=5.0,
    scale=24.0,
    size=95,
)


# Define our statistic.
# A negative result means the candidate is faster.
def median_difference(a, b, *, axis=0):
    return (
        np.median(b, axis=axis)
        - np.median(a, axis=axis)
    )


observed = median_difference(baseline, candidate)


# 95% bootstrap confidence interval.
ci_result = bootstrap(
    data=(baseline, candidate),
    statistic=median_difference,
    paired=False,
    vectorized=True,
    n_resamples=20_000,
    confidence_level=0.95,
    method="percentile",
    batch=1_000,
    rng=np.random.default_rng(101),
)


# Two-sided independent-sample permutation test.
permutation_result = permutation_test(
    data=(baseline, candidate),
    statistic=median_difference,
    permutation_type="independent",
    vectorized=True,
    n_resamples=19_999,
    alternative="two-sided",
    batch=1_000,
    rng=np.random.default_rng(202),
)


print(f"Baseline median: {np.median(baseline):.2f} ms")
print(f"Candidate median: {np.median(candidate):.2f} ms")
print(f"Observed difference: {observed:.2f} ms")

print(
    "95% bootstrap CI: "
    f"[{ci_result.confidence_interval.low:.2f}, "
    f"{ci_result.confidence_interval.high:.2f}] ms"
)

print(
    "Bootstrap standard error: "
    f"{ci_result.standard_error:.2f} ms"
)

print(
    "Permutation p-value: "
    f"{permutation_result.pvalue:.4f}"
)

Representative output with the specified random seeds:

Baseline median: 154.21 ms
Candidate median: 128.66 ms
Observed difference: -25.55 ms
95% bootstrap CI: [-55.47, -8.53] ms
Bootstrap standard error: 12.05 ms
Permutation p-value: 0.0026

The precise resampling results can vary with implementation details and library versions.

Interpret the result

The observed candidate median is approximately 25.55 ms lower than the baseline median.

The percentile bootstrap confidence interval runs from approximately −55.47 ms to −8.53 ms.

Under the sampling assumptions and with adequate bootstrap calibration, the interval provides a range of plausible values for the population median difference.

Since its upper endpoint is below zero, it supports the conclusion that the candidate has a lower population median latency.

The permutation test returns a p-value of approximately 0.0026.

Under the test’s exchangeability null hypothesis, a result at least this extreme would be uncommon.

These are complementary findings:

• The confidence interval expresses the magnitude and uncertainty of the effect.

• The permutation p-value expresses how incompatible the observed statistic is with the null model.

Neither establishes that the latency improvement is necessarily important to users. Statistical evidence and practical significance are different things.

And neither fixes bias caused by an invalid experiment, nonrepresentative sampling, or incorrect assumptions about the independence of observations.

5. Implement the bootstrap from scratch

SciPy handles the machinery, but understanding what it does is important when you need to debug or adapt the method.

Let’s implement the percentile bootstrap ourselves.

Step 1: Resample each group independently

For every bootstrap iteration:

1. Draw 80 observations from the baseline sample, with replacement.

2. Draw 95 observations from the candidate sample, with replacement.

3. Calculate the difference between the two resampled medians.

4. Save the result.

With replacement means that some observations can appear multiple times in a resample while others do not appear at all.

This is intentional. It approximates drawing another random sample from the empirical distribution.

We resample within each group, rather than combining the groups, because we want to estimate uncertainty about their difference without imposing the null hypothesis that their distributions are identical.

Step 2: Construct a confidence interval

After 20,000 bootstrap iterations, we have 20,000 estimates.

For a two-sided 95% percentile interval, take:

• The 2.5th percentile of the estimates as the lower bound.

• The 97.5th percentile as the upper bound.

This is the percentile bootstrap confidence interval.

Notice that we have not needed to derive a standard-error formula for the median difference.

6. Implement the permutation test from scratch

The permutation algorithm is similar computationally but fundamentally different statistically.

Instead of sampling each group with replacement, we combine the observations and randomly reassign them into groups of their original sizes.

For each random reassignment, we calculate the difference of medians.

If the observed value is unusually far into the tails of this null distribution, it provides evidence against the exchangeability null hypothesis.

Complete manual implementation

This example independently regenerates the dataset, so it can run without executing the previous script.

import numpy as np


# Recreate the dataset.
data_rng = np.random.default_rng(2026)

baseline = 50.0 + data_rng.gamma(
    shape=5.0,
    scale=24.0,
    size=80,
)

candidate = 15.0 + data_rng.gamma(
    shape=5.0,
    scale=24.0,
    size=95,
)


def median_difference(a, b):
    return np.median(b) - np.median(a)


observed = median_difference(baseline, candidate)


# PART 1: Bootstrap confidence interval.
bootstrap_rng = np.random.default_rng(101)

n_bootstrap = 20_000
bootstrap_statistics = np.empty(n_bootstrap)

for i in range(n_bootstrap):
    baseline_sample = bootstrap_rng.choice(
        baseline,
        size=len(baseline),
        replace=True,
    )

    candidate_sample = bootstrap_rng.choice(
        candidate,
        size=len(candidate),
        replace=True,
    )

    bootstrap_statistics[i] = median_difference(
        baseline_sample,
        candidate_sample,
    )

ci_low, ci_high = np.quantile(
    bootstrap_statistics,
    [0.025, 0.975],
)


# PART 2: Permutation hypothesis test.
permutation_rng = np.random.default_rng(202)

pooled = np.concatenate([baseline, candidate])
n_baseline = len(baseline)

n_permutations = 19_999
extreme_count = 0

for i in range(n_permutations):
    shuffled = permutation_rng.permutation(pooled)

    permuted_baseline = shuffled[:n_baseline]
    permuted_candidate = shuffled[n_baseline:]

    permuted_statistic = median_difference(
        permuted_baseline,
        permuted_candidate,
    )

    if abs(permuted_statistic) >= abs(observed):
        extreme_count += 1


# Plus-one correction for a randomized test.
p_value = (
    extreme_count + 1
) / (
    n_permutations + 1
)


print(f"Observed difference: {observed:.2f} ms")
print(f"Bootstrap CI: [{ci_low:.2f}, {ci_high:.2f}] ms")
print(f"Permutation p-value: {p_value:.4f}")
print(f"Extreme permutations: {extreme_count}")

Representative output:

Observed difference: -25.55 ms
Bootstrap CI: [-54.99, -8.78] ms
Permutation p-value: 0.0021
Extreme permutations: 41

The results closely resemble the SciPy implementation. The confidence intervals differ slightly because the two implementations consume random numbers differently.

The p-values differ for an additional reason: the two implementations use different definitions of a two-sided test.

Why do the permutation p-values differ?

Our manual test counts how often a permutation produces a statistic whose absolute magnitude is at least as large as the observed statistic.

In plain notation:

p ≈ (number of extreme permutations + 1) / (number of permutations + 1)

Here, extreme means the absolute permuted statistic is at least the absolute observed statistic.

SciPy’s two-sided permutation test uses another common convention: it doubles the smaller of the two one-sided tail probabilities, with the result capped at 1.

These approaches coincide for many symmetric null distributions but can differ for asymmetric distributions.

The distinction matters when your sample sizes differ or the statistic’s null distribution is asymmetric.

There is no universal definition of two-sided extremeness for every statistic. The correct approach depends on how the alternative hypothesis is formulated and how the test is calibrated.

For reproducible analysis, document the definition you choose.

Why add one to the numerator and denominator?

Suppose you generate 999 random permutations and none produces a statistic as extreme as the observed value.

The naive empirical estimate would be zero.

But zero is not a sensible Monte Carlo p-value for a finite randomized permutation test. It implies more precision than the simulation provides.

The plus-one adjustment counts the observed assignment along with the randomly generated assignments, giving a minimum p-value of 1 / 1000 in this example.

For Monte Carlo permutation tests under the appropriate null assumptions, this convention produces a valid conservative p-value.

It also prevents reporting a misleading zero.

7. Understand what a confidence interval actually means

A confidence interval is frequently interpreted incorrectly.

Suppose your 95% bootstrap interval is:

• Lower bound: −55 ms

• Upper bound: −9 ms

It is tempting to conclude that there is a 95% probability that the true improvement falls inside those bounds.

That is not the standard frequentist interpretation.

The true population parameter is treated as fixed. The interval is random because it is constructed from random data.

A properly calibrated 95% confidence-interval procedure should contain the population parameter approximately 95% of the time across repeated independent studies under its assumptions.

For the particular interval you computed, the parameter either lies inside or outside it.

Confidence intervals depend on the sampling model

Our independent bootstrap treats observations within each group as independent draws from their respective populations.

If that assumption is inappropriate, the resulting interval may have incorrect coverage.

For example, imagine that 80 latency measurements actually come from only four independent servers, each contributing 20 strongly correlated observations.

Treating all 80 measurements as independent can substantially underestimate uncertainty.

A confidence interval cannot correct an incorrectly modeled sampling process simply by increasing the number of resamples.

Why the sample statistic matters

The bootstrap is flexible because the target statistic is programmable.

However, not every statistic has equally good bootstrap behavior.

The median is more complicated than the mean because it is determined by the ordering of observations. Small data changes can move the sample median abruptly.

For large samples from sufficiently regular continuous distributions, bootstrap inference for medians is generally well motivated.

For very small samples, highly discrete observations, or distributions with substantial probability mass at the median, percentile intervals can have poor coverage.

This is an important limitation of the promise of inference without analytical formulas: avoiding a formula does not eliminate statistical regularity conditions.

8. Percentile, basic, and BCa intervals

So far, we’ve used the percentile bootstrap because its construction is intuitive.

SciPy also supports two alternatives worth knowing.

Percentile interval

The percentile interval takes quantiles directly from the bootstrap distribution.

Advantages:

• Easy to understand and implement.

• Preserves the natural scale of the statistic.

• Requires relatively little additional computation.

Limitations:

• Can exhibit coverage errors when sampling distributions are biased or skewed.

• Can perform poorly with small or highly discrete datasets.

• Does not explicitly correct for the estimator’s bias or acceleration.

It is a useful starting point, but not a universal best choice.

Basic bootstrap interval

The basic, or reverse-percentile, interval reflects bootstrap quantiles around the observed estimate.

If the bootstrap distribution is asymmetric, the basic and percentile intervals need not agree.

The basic interval can sometimes perform better, but it is not guaranteed to outperform the percentile method.

It also behaves differently under nonlinear transformations of the parameter.

Bias-corrected and accelerated interval

The BCa method adjusts the interval’s quantile levels to account for estimated bias and changes in the estimator’s sampling variability.

It often improves coverage accuracy for sufficiently regular problems.

SciPy currently uses BCa as the default bootstrap interval method.

However, BCa has practical limitations:

• It requires additional computation, including jackknife calculations.

• It can become unstable for certain small or highly discrete samples.

• A degenerate bootstrap distribution can produce undefined interval endpoints.

• Bias correction cannot repair dependence, confounding, or nonrepresentative sampling.

For advanced statistical work, calculate more than one interval method when uncertainty calibration matters.

Large differences between percentile, basic, and BCa intervals may indicate that the estimated uncertainty is sensitive to the method.

In that situation, investigate the sampling distribution and assumptions rather than selecting whichever interval supports the preferred conclusion.

9. Paired data require a different resampling strategy

Consider a different experiment.

You execute the same 30 test workloads against both API versions, recording the latency for each workload.

Now every baseline observation is naturally paired with a corresponding candidate observation.

These are not two independent samples.

A particularly useful statistic is the median of the within-pair latency differences.

This is not generally equal to the difference between the two group medians.

Correct paired bootstrap and permutation test

For a paired bootstrap, resample entire pairs together.

For a paired permutation test, randomly swap the baseline and candidate labels within each pair.

Here is a self-contained example of such a matched-workload experiment:

import numpy as np
from scipy.stats import bootstrap, permutation_test


rng = np.random.default_rng(77)

baseline = 45.0 + rng.lognormal(
    mean=3.0,
    sigma=0.4,
    size=30,
)

candidate = (
    baseline
    - 8.0
    + rng.normal(0.0, 4.0, size=30)
)


# Median of paired differences, not difference of medians.
def paired_median_difference(a, b, *, axis=0):
    return np.median(b - a, axis=axis)


observed = paired_median_difference(
    baseline,
    candidate,
)


# Resample matched pairs together.
ci_result = bootstrap(
    (baseline, candidate),
    paired_median_difference,
    paired=True,
    vectorized=True,
    n_resamples=10_000,
    method="percentile",
    rng=np.random.default_rng(301),
)


# Swap labels within pairs under the null.
test_result = permutation_test(
    (baseline, candidate),
    paired_median_difference,
    permutation_type="samples",
    vectorized=True,
    n_resamples=9_999,
    alternative="two-sided",
    rng=np.random.default_rng(302),
)


print(f"Median paired effect: {observed:.2f} ms")

print(
    "95% paired bootstrap CI: "
    f"[{ci_result.confidence_interval.low:.2f}, "
    f"{ci_result.confidence_interval.high:.2f}] ms"
)

print(f"Paired permutation p-value: {test_result.pvalue:.4f}")

Two details are especially important.

Our run (NumPy 2.5.4, SciPy 1.18.1) printed:

Median paired effect: -8.03 ms
95% paired bootstrap CI: [-9.43, -6.41] ms
Paired permutation p-value: 0.0002

Pairing removes most of the between-request variation, so the interval is far narrower than in the independent-samples analysis.

First, the paired bootstrap uses the same resampled indices for both arrays. This preserves the within-pair relationship.

Second, the permutation test exchanges labels within each pair rather than pooling all observations globally.

Such label swapping is justified when the experimental design or null hypothesis makes the two labels exchangeable within pairs. For a randomized matched-pair experiment, this can follow directly from the assignment mechanism.

It is not automatically justified for an observational before-and-after comparison, where temporal trends or other changes may confound the result.

The correct resampling unit is a property of your experimental design, not merely a software setting.

10. The biggest trap: equal medians do not imply exchangeability

This is arguably the most important subtlety in the entire tutorial.

Suppose two groups have identical population medians but substantially different variances or distribution shapes.

For example:

• Group A has a narrow distribution concentrated around 100 ms.

• Group B has a much wider distribution, also centered around a median of 100 ms.

The population median difference is zero.

But the two distributions are not identical.

An ordinary independent-sample permutation test assumes that observations could exchange group labels under its null model.

That assumption may fail even though the population medians are equal.

Consequently, an unstudentized permutation test using the difference of medians is not generally an exact test of the hypothesis that the two population medians are equal.

What does this mean in practice?

If your scientific null hypothesis is that the complete distributions are identical, the standard independent-sample permutation design is appropriate under the relevant exchangeability assumptions.

If you specifically want to test equality of medians while allowing other distributional features to differ, the standard permutation test may have incorrect type-I error.

Alternatives include:

• Methods developed specifically for quantile comparisons under unequal distributions.

• Appropriately calibrated or studentized resampling statistics.

• Quantile regression when covariates and experimental structure matter.

• A randomization test based on the actual assignment mechanism when a sharp causal null is the target.

Studentization can improve asymptotic validity in some settings, but it must be justified for the statistic and sampling design.

It is not a universal repair.

The broader principle is that a permutation test is only as valid as the exchangeability argument underlying its rearrangements.

11. How many resamples should you use?

You must distinguish two sources of uncertainty:

1. Statistical uncertainty because your original dataset is finite.

2. Monte Carlo uncertainty because you approximate a distribution using finitely many random resamples.

Increasing the resample count reduces the second source, not the first.

Bootstrap resample counts

For a routine 95% interval, several thousand bootstrap replicates can be an adequate starting point.

For more stable tail quantiles, 10,000 to 20,000 replicates are reasonable.

If interval endpoints are important to a high-stakes decision, assess their Monte Carlo stability by repeating the procedure with different random seeds or increasing the number of replicates.

The stability of an endpoint depends on the distribution and quantile being estimated, not just the nominal number of draws.

Permutation counts

For a randomized permutation test, the smallest attainable p-value using the plus-one convention is:

Minimum p-value = 1 / (number of random permutations + 1)

For example:

If you need to investigate extremely small p-values, more random permutations may be necessary.

Even with 19,999 permutations, a p-value near 0.003 has noticeable Monte Carlo variability in its final decimal places.

Avoid reporting more significant digits than the simulation supports.

Exact permutation tests

For small datasets, you may be able to enumerate every distinct reassignment.

With three observations in one group and three in another, there are only 20 distinct ways to choose which three of the six pooled observations go into the first group.

An exact enumeration eliminates Monte Carlo error in the permutation distribution.

It does not eliminate the need for a valid null hypothesis and exchangeability assumptions.

For larger datasets, the number of arrangements grows combinatorially, so random sampling of permutations is usually more practical.

12. Performance engineering for resampling

Resampling can become computationally expensive when the statistic is costly.

A 20,000-replicate bootstrap of a median is typically manageable. A 20,000-replicate bootstrap of an expensive optimization routine or complex machine-learning pipeline can be far more demanding.

Several practices help.

Vectorize the statistic

NumPy supports median calculations along a specified axis.

SciPy can exploit this to evaluate many resamples in batches without repeatedly invoking Python-level statistic functions for each individual resample.

The examples above define an axis-aware statistic for that reason.

Not every statistic can be vectorized conveniently, but it is worth considering when performance matters.

Bound memory usage

Vectorizing an entire resampling computation at once can allocate large intermediate arrays.

SciPy’s batch parameter lets you limit the number of resamples processed together.

A moderate batch size can provide much of the speed benefit of vectorization while avoiding unnecessarily large allocations.

Preserve reproducibility

Record:

• NumPy and SciPy versions.

• The statistic definition and its sign convention.

• Resampling method and relevant parameters.

• The number of resamples.

• Random seeds.

• The experimental unit and any dependence structure.

• The exact hypotheses and confidence level.

A reported p-value without its test definition and resampling design can be surprisingly difficult to reproduce.

Avoid unnecessary recomputation

If you need several confidence levels from the same bootstrap distribution, you do not necessarily need to generate all the resamples again.

The bootstrap result returned by SciPy includes the calculated bootstrap statistics and supports reusing a previous result.

For manual implementations, you can retain the bootstrap estimates and compute different quantiles from that distribution.

13. Common mistakes that invalidate resampling analysis

Resampling is powerful, but simple implementation errors can create convincing yet misleading results.

Mistake 1: Mixing up replacement rules

The ordinary nonparametric bootstrap samples with replacement.

A permutation test reassigns the existing observations without replacement.

Confusing these operations changes the reference distribution and the statistical question being answered.

Mistake 2: Pooling groups for the ordinary bootstrap CI

Pooling two treatment groups before generating a bootstrap confidence interval for their difference imposes a shared empirical distribution.

That is generally inappropriate when estimating uncertainty about the difference between two potentially distinct populations.

Resample independently within the groups for the ordinary independent-sample bootstrap.

Mistake 3: Ignoring repeated observations

Server telemetry, customer transactions, repeated sensor readings, and time-series measurements often contain dependence.

If thousands of requests come from a small number of servers or users, a naive observation-level bootstrap can underestimate uncertainty.

Depending on the design, consider:

• Cluster bootstrap methods that resample independent clusters.

• Block bootstrap methods that preserve important temporal dependence.

• Stratified or design-aware resampling when sampling or assignment uses strata.

These approaches require choosing a resampling unit that reflects how uncertainty arises in the original data.

Mistake 4: Interpreting the p-value as a probability that the null is true

A p-value of 0.0026 does not mean there is a 0.26% probability that the null hypothesis is true.

It measures how extreme the result is relative to the specified null distribution.

It also does not imply that there is only a 0.26% chance the result occurred because of noise.

Those interpretations reverse the conditional reasoning of the test.

Mistake 5: Treating the observed sample as universally representative

Bootstrap methods cannot discover missing customer segments, correct selection bias, or repair broken measurement instrumentation.

They reuse the information already present in your sample.

If the sample does not adequately represent the target population, the bootstrap distribution can provide precise-looking estimates for the wrong target.

Mistake 6: Ignoring multiplicity

Suppose an engineering team tests 40 endpoints and publishes only the one with the smallest p-value.

A small unadjusted p-value from that selected comparison no longer carries the same evidential interpretation as a single preregistered test.

When screening many statistics, account for multiple comparisons or use an appropriate error-control procedure.

Resampling can help with multiple-testing procedures, but ordinary per-metric bootstrap intervals do not automatically provide simultaneous coverage.

Mistake 7: Treating statistical significance as sufficient for deployment

Suppose a new service architecture reduces median latency by 2 ms with an extremely small p-value.

If it also substantially increases operating costs and complicates reliability, the improvement may not justify deployment.

Conversely, an interval may be too wide to reject a null hypothesis while still containing improvements large enough to be commercially interesting.

Always report effect sizes, uncertainty, and practical decision thresholds together.

14. Cherry on the cake: the bootstrap was born partly because the jackknife could fail for the median

One of the most interesting facts about the bootstrap is that the median—our central example—played a role in its original development.

In 1979, statistician Bradley Efron published Bootstrap Methods: Another Look at the Jackknife in The Annals of Statistics.

At the time, the jackknife was already an established resampling technique.

Its basic idea was to remove one observation at a time, recompute the statistic, and use the resulting variation to estimate uncertainty.

For many sufficiently smooth statistics, this works well.

But Efron highlighted a surprising result: the ordinary delete-one jackknife variance estimator can be inconsistent for the sample median, even under regular conditions.

In other words, simply increasing the sample size does not guarantee that this jackknife variance estimate converges to the correct target.

The ordinary nonparametric bootstrap, meanwhile, can consistently estimate the sampling variance of the median under suitable regularity conditions.

Why?

The sample median is not a smooth linear function of individual observations. Removing one observation often changes it very little, while other deletions can move it abruptly between adjacent order statistics.

The leave-one-out changes do not necessarily capture the correct sampling variability.

Bootstrap resampling replaces multiple observations in each replicate, producing a different approximation to sampling variability.

The result was more than a historical curiosity. It showed that two procedures commonly grouped under the broad heading of resampling could have radically different statistical behavior.

The lesson remains relevant today: repeated computation is not, by itself, evidence of a valid uncertainty estimate. The resampling scheme has to match the statistic and the sampling model.

15. A practical workflow for production analysis

When analyzing an unconventional statistic, follow a disciplined sequence.

First, define the estimand.

Be explicit about what population quantity you want to estimate.

The difference between population medians is not the same as the population median of paired differences. Choosing a different statistic changes the question.

Second, understand the data-generating process.

Identify whether observations are independent, paired, clustered, stratified, or serially dependent.

The resampling unit should preserve the structure that matters for statistical uncertainty.

Third, calculate and inspect the observed effect.

Check sample sizes, missing observations, extreme values, distribution shapes, and whether the groups are comparable.

For a difference of medians, inspect each group’s distribution rather than looking only at its median.

Fourth, construct an interval suited to the estimand.

Start with a transparent bootstrap implementation. Compare alternative interval methods or assess empirical coverage through simulation when the problem warrants it.

Fifth, formulate the null hypothesis before running the test.

Decide whether you are testing identical distributions, a sharp randomization null, equality of a specific population statistic, or something else.

Only then choose a valid permutation or other testing procedure.

Sixth, check numerical stability.

Increase resample counts, use independent seeds, and inspect whether reported conclusions depend heavily on Monte Carlo noise or the selected interval method.

Finally, connect inference to the decision.

Report the estimated effect, confidence interval, p-value when appropriate, sampling assumptions, and the magnitude of improvement needed to justify an action.

16. Final takeaways

Bootstrap confidence intervals and permutation tests make statistical inference possible for many metrics where classical analytical calculations are inconvenient, fragile, or unavailable.

The bootstrap approximates the variation of an estimator under repeated sampling. Permutation testing constructs a reference distribution under a particular null hypothesis.

Both can work with custom statistics, including differences of medians, percentiles, trimmed means, and domain-specific scores.

But neither is assumption-free.

A confidence interval can miss its nominal coverage, and a permutation p-value can be invalid if its exchangeability assumptions are violated.

The central best practice is not simply to replace mathematical formulas with computation. It is to combine a clearly defined statistic, a defensible resampling model, and careful interpretation.

For further reference, the current SciPy bootstrap guide, SciPy permutation testing guide, and Efron’s foundational paper provide the API details and statistical background.

Your next step: run the complete SciPy example, then replace the difference of medians with the 90th-percentile latency difference for your own API or service telemetry. Compare its bootstrap interval and permutation p-value with the median-based results. Investigate whether the same deployment decision is supported by both metrics—and, most importantly, whether your resampling assumptions actually match how those measurements were collected.