Bootstrap and Permutation Methods Card Template (Statistics, R and Python)

Resampling Methods: [Your study question or course]

A method card for resampling: bootstrap confidence intervals (percentile, basic and BCa), permutation tests (exact and Monte Carlo) and Monte Carlo simulation for checking whether a method works on data like yours. It is for students meeting the bootstrap for the first time, thesis writers who need an interval for a statistic with no textbook formula, and analysts who want a check that does not lean on normality. It speeds up learning by putting the decisions in order (what to resample, how many times, which interval, what the bootstrap cannot fix), then making you practise on a fully solved, a half solved and an unsolved example, with R and Python run on the same data. The worked examples use fictional data; replace them with your own. Source basis: Efron (1979), Efron and Tibshirani (1993), Davison and Hinkley (1997), the NIST/SEMATECH e-Handbook and the R boot and SciPy stats documentation.

Document control

When to use it

One-sentence rule: resample from your own data to see how much a statistic would vary from sample to sample (the bootstrap, for intervals), or shuffle group labels to see what differences arise when the labels carry no information (the permutation test, for p-values).

The bootstrap treats the sample as a stand-in for the population. Drawing many new samples of the same size, with replacement, from that stand-in gives a cloud of recomputed statistics whose spread approximates the standard error and whose shape suggests an interval. The permutation test answers a narrower question: if group membership were irrelevant, how unusual would the observed difference be among all the ways the same numbers could have been split into groups? Monte Carlo simulation is the third tool here: instead of resampling your data, you generate many artificial data sets from a known truth and check whether a method delivers what it promises (95% coverage, a 5% false-positive rate).

Alternatives: a t interval or Welch t-test when the mean is the target and n is moderate (academic/stats-t-tests); rank-based tests (academic/stats-nonparametric-tests); exact or Wilson intervals for a proportion (academic/stats-proportion-tests); block, cluster or residual bootstraps for dependent data; Bayesian intervals (academic/stats-bayesian-inference).

The question and the data

Question for Example A: how long do students typically wait for a drink at the Lakemoor University campus cafe at lunchtime, and how precise is that estimate? Wait times are right-skewed (most are short, a few are long), so the median is a sensible summary, and the median has no simple standard error formula for small samples.

Question for Example B: does a flipped tutorial format give higher quiz scores than a lecture format? Sixteen volunteers in a Lakemoor University statistics course were randomly assigned, eight to each format.

Example A data (n = 15 cafe visits, one per student, fictional):

Summary: mean 4.593 minutes, median 3.4 minutes. The mean sits well above the median because of the 12.5 and 9.4 minute waits.

Example B data (quiz score, fictional):

Assumptions and how to check each

How it works

The bootstrap in plain language. You have one sample and want to know how much the median would change if you could repeat the study. You cannot repeat the study, so you pretend the sample is the population. You draw 15 values from your 15, with replacement (some values appear twice, some not at all), compute the median, and repeat thousands of times. The spread of those medians estimates the standard error, and the middle 95% of them gives a first interval. Efron (1979) introduced the idea and showed that it reproduces many classical results while extending to statistics with no formula.

Three intervals from the same bootstrap distribution. Let $\hat\theta$ be the statistic in the original sample and $\hat\theta^*_{(\alpha)}$ the $\alpha$ quantile of the $B$ bootstrap replicates.

The percentile interval reads the quantiles directly. The basic interval reflects them around $\hat\theta$, because it assumes the error $\hat\theta^*-\hat\theta$ mimics the error $\hat\theta-\theta$. The BCa (bias corrected and accelerated) interval moves the percentile cut-offs to correct for median bias $z_0$ and for skewness through an acceleration constant $a$ estimated by the jackknife:

and similarly $\alpha_2$ with $z_{1-\alpha/2}$. The BCa interval is $[\hat\theta^*_{(\alpha_1)},\hat\theta^*_{(\alpha_2)}]$. When $z_0=0$ and $a=0$ it collapses to the percentile interval.

The permutation test in plain language. Under H0 the format label does not matter, so any split of the 16 scores into two groups of 8 was as likely as the one observed. List every split (or a large random sample of splits), compute the difference in means for each, and count how often a split gives a difference at least as extreme as the observed one.

The "+1" in the Monte Carlo version counts the observed split itself, which keeps the p-value from being exactly zero and keeps the test valid.

Monte Carlo simulation for checking a method. Choose a population you can write down (for example, exponential waits with true mean 4 minutes), draw thousands of samples of your size, apply the method to each, and record whether the interval caught the true value or whether the test rejected a true H0. A 95% interval should catch the truth in about 95% of runs; a test at α = .05 should reject a true H0 about 5% of the time. The simulation has its own uncertainty: with R runs, the Monte Carlo standard error of a proportion near 0.95 is about $\sqrt{0.95 \times 0.05 / R}$, which is 0.0049 for 2,000 runs.

Seeds and reproducibility. Every resampling result depends on the random stream. Setting a seed (set.seed(2026) in R, np.random.default_rng(2026) in Python) makes a run repeatable on the same software. It does not make R and Python agree: they use different generators (R defaults to Mersenne Twister with its own sampling algorithm, NumPy's default_rng uses PCG64), so the same seed produces different resamples. Different quantile rules also shift the endpoints. Expect R and Python bootstrap intervals to differ in the second decimal place, and more for lumpy statistics such as a small-sample median. Report the software, seed and B so a reader can reproduce your exact numbers, and choose B large enough that a different seed would not change your conclusion.

Worked example A (full)

Question: a 95% interval for the median (and, for comparison, the mean) wait time at the campus cafe.

Hypotheses: none; this is estimation. The target is the population median wait.

Checks: each wait comes from a different student on a different day, so independence is plausible. The data are right-skewed with two long waits. There are 15 distinct values, so the bootstrap median can take at most 15 values: expect a lumpy bootstrap distribution and coarse interval endpoints.

R code (install the package once with install.packages("boot")):

wait <- c(2.1, 3.4, 1.8, 7.9, 2.6, 4.2, 3.0, 12.5, 2.2, 5.1, 3.7, 2.9, 9.4, 3.3, 4.8)
library(boot)
med_fun <- function(x, i) median(x[i])   # i = row indices of one resample
set.seed(2026)
b <- boot(wait, statistic = med_fun, R = 9999)
b
boot.ci(b, conf = 0.95, type = c("perc", "basic", "bca"))
table(b$t)                                # how lumpy is the bootstrap median?

mean_fun <- function(x, i) mean(x[i])
set.seed(2026)
bm <- boot(wait, statistic = mean_fun, R = 9999)
boot.ci(bm, conf = 0.95, type = c("perc", "bca"))
t.test(wait)$conf.int

Output (abridged and relabelled):

Bootstrap Statistics :
    original   bias    std. error
t1*      3.4 0.169507   0.6250787

Intervals (median):
Level      Basic              Percentile            BCa
95%   ( 2.0,  3.9 )   ( 2.9,  4.8 )   ( 2.6,  4.8 )

 2.1  2.2  2.6  2.9    3  3.3  3.4  3.7  4.2  4.8  5.1  7.9  9.4
   1   38  207  641 1287 1820 2090 1820 1247  613  192   41    2

Intervals (mean):
Level     Percentile            BCa
95%   ( 3.267,  6.200 )   ( 3.440,  6.635 )

t interval for the mean: 2.905634 6.281033

Python code:

import numpy as np
from scipy import stats

wait = np.array([2.1, 3.4, 1.8, 7.9, 2.6, 4.2, 3.0, 12.5, 2.2, 5.1, 3.7, 2.9, 9.4, 3.3, 4.8])
for m in ["percentile", "basic", "BCa"]:
    r = stats.bootstrap((wait,), np.median, n_resamples=9999, method=m,
                        rng=np.random.default_rng(2026))
    print(m, r.confidence_interval, round(r.standard_error, 4))
for m in ["percentile", "BCa"]:
    r = stats.bootstrap((wait,), np.mean, n_resamples=9999, method=m,
                        rng=np.random.default_rng(2026))
    print("mean", m, r.confidence_interval)

Output (abridged and reformatted):

median percentile CI = (2.600, 5.100)  SE = 0.6288
median basic      CI = (1.700, 4.200)  SE = 0.6288
median BCa        CI = (2.900, 5.100)  SE = 0.6288
mean   percentile CI = (3.260, 6.220)
mean   BCa        CI = (3.467, 6.647)

The bootstrap distribution of the median from the R run shows why the endpoints jump between software:

Calculation by hand for the key statistic. The basic interval comes straight from the percentile interval: with $\hat\theta = 3.4$ and R's percentile limits 2.9 and 4.8, the basic limits are $2(3.4) - 4.8 = 2.0$ and $2(3.4) - 2.9 = 3.9$, exactly the R output. Note how the basic interval (2.0, 3.9) is shifted left of the percentile interval: reflecting a right-skewed bootstrap distribution pushes the interval the opposite way. That is correct only if the bootstrap error $\hat\theta^*-\hat\theta$ really mimics the sampling error $\hat\theta-\theta$, which is doubtful for a lumpy small-sample median, so the basic interval is a risky choice here. The bootstrap SE (0.625 in R, 0.629 in Python) is simply the standard deviation of the 9,999 bootstrap medians.

Why R and Python disagree. Both runs use seed 2026 and B = 9,999, yet the percentile interval is (2.9, 4.8) in R and (2.6, 5.1) in Python. Two reasons: the random streams differ, and the bootstrap median can only take the 13 values in the table, so a quantile that falls near the border between 2.6 and 2.9 (about 2.5% of resamples lie at 2.6 or below) flips to the neighbouring value with a small change in counts or in the quantile rule (boot uses an order-statistic rule at $(B+1)q$, SciPy uses NumPy's percentile). With n = 15 the honest conclusion is that the median interval is coarse: roughly 2.6 to 5 minutes.

Effect size with CI. For a one-sample estimation problem the "effect" is the estimate itself: median 3.4 minutes, 95% BCa CI [2.6, 4.8] (R). For the mean, 4.59 minutes, BCa CI [3.44, 6.64], compared with the t interval [2.91, 6.28]. The BCa interval for the mean is shifted right and is asymmetric, matching the right skew; the t interval is symmetric by construction.

Interpretation. A typical lunchtime wait is about three and a half minutes, and the data are consistent with a population median between roughly 2.6 and 4.8 minutes. The mean is higher (about 4.6 minutes) because a few long waits pull it up; for planning staff cover, the long tail matters more than the median.

APA-style report sentence. The median wait was 3.4 minutes (n = 15), 95% BCa bootstrap CI [2.6, 4.8], based on 9,999 resamples (R boot 1.3.31, seed 2026); the mean wait was 4.59 minutes, 95% BCa CI [3.44, 6.64].

Worked example B (half solved)

Question: do flipped-format students score higher than lecture-format students? Two-sided permutation test of the difference in means, α = .05.

Hypotheses: H0: the format label is unrelated to scores (the 16 scores would be the same whichever format each student received). H1: [ ].

Observed statistic: $D_{\text{obs}} = 83.125 - 75.000 =$ [ ].

Number of possible splits: $\binom{16}{8} =$ [ ], small enough to list them all, so an exact test is possible.

R code:

flipped <- c(78, 85, 81, 90, 74, 88, 83, 86)
lecture <- c(72, 80, 69, 77, 84, 70, 75, 73)
obs <- mean(flipped) - mean(lecture)
pooled <- c(flipped, lecture)
idx <- combn(16, 8)                                    # every split, one column each
diffs <- apply(idx, 2, function(g) mean(pooled[g]) - mean(pooled[-g]))
mean(abs(diffs) >= abs(obs) - 1e-9)                    # exact two-sided p
set.seed(2026)                                         # Monte Carlo version
mc <- replicate(9999, { g <- sample(16, 8); mean(pooled[g]) - mean(pooled[-g]) })
(sum(abs(mc) >= abs(obs) - 1e-9) + 1) / (9999 + 1)
t.test(flipped, lecture)                               # Welch, for comparison

Output (abridged and relabelled; obs is 8.125):

observed difference = 8.125
exact splits = 12870  exact two-sided p = 0.01134421
Monte Carlo two-sided p = 0.0125
Welch Two Sample t-test: t = 3.1155, df = 13.984, p-value = 0.007604
95 percent confidence interval: 2.53093 13.71907

Python code:

import numpy as np
from scipy import stats
flipped = np.array([78, 85, 81, 90, 74, 88, 83, 86])
lecture = np.array([72, 80, 69, 77, 84, 70, 75, 73])
def diff(x, y, axis=-1):
    return np.mean(x, axis=axis) - np.mean(y, axis=axis)
exact = stats.permutation_test((flipped, lecture), diff, permutation_type="independent",
                               n_resamples=np.inf, alternative="two-sided")
mc = stats.permutation_test((flipped, lecture), diff, permutation_type="independent",
                            n_resamples=9999, alternative="two-sided",
                            rng=np.random.default_rng(2026))
print(exact.pvalue, exact.null_distribution.size, mc.pvalue)

Output (abridged and relabelled):

exact: statistic = 8.125, p = 0.01134, n null = 12870
Monte Carlo: p = 0.0116

Your turn.

1. How many of the 12,870 splits gave $|D_\pi| \ge 8.125$? (Multiply the exact p by 12,870.) [ ]
2. Why do the exact p-values agree between R and Python but the Monte Carlo p-values (0.0125 and 0.0116) do not? [ ]
3. Decision at α = .05: [ ]
4. Report sentence: "Flipped-format students scored [ ] points higher on average (M = 83.1 vs 75.0), exact permutation test p = [ ], 12,870 splits." [ ]

Exercise C (unsolved)

Using the Example B data, build a 95% percentile bootstrap interval for the difference in mean quiz scores (flipped minus lecture). Resample within each group separately (8 from the flipped scores, 8 from the lecture scores, each with replacement), because the design fixed the group sizes. Use B = 9,999 and seed 2026.

1. Write the R code (hint: replicate with two sample(..., replace = TRUE) calls) and the Python code (hint: stats.bootstrap((flipped, lecture), statistic, method="percentile") with a two-argument statistic).
2. Report the interval and the bootstrap SE from each language.
3. Compare with the Welch 95% CI [2.53, 13.72]. Which is narrower, and why might the bootstrap be too narrow with only 8 per group?
4. Does the interval exclude zero, and does that agree with the permutation test?

Reading software output

Reporting template

Fill-in sentences.

• Bootstrap: "The [statistic] of [variable] was [value] (n = [n]), 95% [BCa / percentile] bootstrap CI [[lower], [upper]], based on [B] resamples of [unit] ([software and version], seed [seed])."
• Permutation: "[Group 1] scored [difference] [units] higher than [group 2] (M = [m1] vs [m2]); [exact / Monte Carlo, B = [B]] permutation test of the difference in means, two-sided p = [p]."
• Simulation check: "In [R] simulated data sets from [population], the [method] interval covered the true value in [coverage]% of runs (Monte Carlo SE [se])."

Results checklist.

• The resampling unit matches the independent unit of the design.
• The interval type was chosen before seeing the results.
• B, seed, software and version are reported.
• The bootstrap distribution was plotted and is not badly lumpy or multimodal.
• For permutation tests: exact or Monte Carlo is stated, and exchangeability is justified by the design.
• A rerun with a different seed does not change the conclusion.
• Limitations for small n or dependence are stated.

Common mistakes and misconceptions

• "The bootstrap creates new data." It does not. It reuses your n values. With n = 8, the bootstrap world has at most 8 distinct values, and intervals tend to be too narrow (the simulation below found about 88% to 89% coverage for a nominal 95% percentile interval at n = 15 with skewed data).
• Resampling rows when units are grouped. If each student contributes five quiz scores, resampling scores treats them as independent and gives an SE that is too small. Resample students.
• Bootstrapping a time series row by row. It destroys the autocorrelation. Use a block bootstrap (academic/stats-time-series-arima).
• Picking the interval that looks best. Percentile, basic and BCa can give different answers; choosing after the fact is a form of p-hacking. Decide in the analysis plan.
• "A permutation test makes no assumptions." It assumes exchangeability under H0. If the two groups have very different spreads, a permutation test of the mean difference can reject because of the spread, not the means.
• Forgetting the +1. A Monte Carlo p-value of 0/9,999 should be reported as p < .001 using (0 + 1)/(9,999 + 1), never as p = 0.
• Expecting the same numbers in R and Python with the same seed. The generators differ. Agreement to about the Monte Carlo error is the right expectation.
• Too few resamples. B = 200 is fine for a rough SE but not for BCa endpoints or for a p-value near .05.
• Using the bootstrap to rescue a biased sample. Resampling a convenience sample reproduces its bias faithfully.

Monte Carlo check you can run. This simulation asks how well two 95% intervals for a mean work with n = 15 from a skewed (exponential) population with true mean 4, and whether the Welch t-test keeps its 5% Type I error with n = 8 per group from the same skewed population.

set.seed(2026)
reps <- 2000; n <- 15; mu <- 4; B <- 999
hit_t <- hit_b <- logical(reps)
for (r in 1:reps) {
  x <- rexp(n, rate = 1 / mu)
  ci_t <- t.test(x)$conf.int
  bm <- replicate(B, mean(sample(x, replace = TRUE)))
  ci_b <- quantile(bm, c(0.025, 0.975))
  hit_t[r] <- ci_t[1] <= mu && mu <= ci_t[2]
  hit_b[r] <- ci_b[1] <= mu && mu <= ci_b[2]
}
mean(hit_t); mean(hit_b)
set.seed(2026)
p <- replicate(10000, t.test(rexp(8), rexp(8))$p.value)
mean(p < 0.05)

Output (abridged and relabelled):

coverage t interval: 0.9205
coverage percentile bootstrap: 0.894
Welch Type I error rate at alpha = .05: 0.0343

import numpy as np
from scipy import stats
rng = np.random.default_rng(2026)
reps, n, mu, B = 2000, 15, 4.0, 999
hit_t = np.empty(reps, bool); hit_b = np.empty(reps, bool)
for r in range(reps):
    x = rng.exponential(scale=mu, size=n)
    se = x.std(ddof=1) / np.sqrt(n); tq = stats.t.ppf(0.975, n - 1)
    hit_t[r] = x.mean() - tq * se <= mu <= x.mean() + tq * se
    bm = rng.choice(x, size=(B, n), replace=True).mean(axis=1)
    lo, hi = np.percentile(bm, [2.5, 97.5])
    hit_b[r] = lo <= mu <= hi
print(hit_t.mean(), hit_b.mean())
rng = np.random.default_rng(2026)
x = rng.exponential(size=(10000, 8)); y = rng.exponential(size=(10000, 8))
print((stats.ttest_ind(x, y, axis=1, equal_var=False).pvalue < 0.05).mean())

Output (abridged and relabelled):

coverage t interval: 0.9070
coverage percentile bootstrap: 0.8835
Welch Type I error rate at alpha = .05: 0.0389

Both languages tell the same story even though the numbers differ: at n = 15 with skewed data, neither interval reaches 95%, and the percentile bootstrap is the worse of the two (about 88% to 89%). The R and Python coverage of the t interval differ by 0.0135, about 1.5 Monte Carlo standard errors of a difference (each coverage estimate near 0.91 has SE about $\sqrt{0.91 \times 0.09/2000} = 0.0064$, so a difference has SE about 0.009), which is ordinary simulation noise from different random streams. The Welch test was conservative here (rejecting a true H0 about 3.4% to 3.9% of the time, not 5%). Lesson: the bootstrap is not a small-sample cure, and a quick simulation shaped like your own data tells you more than a rule of thumb.

Self-test questions

1. You bootstrap the mean of 40 weekly sales figures from one shop, in time order. What is wrong, and what would you do instead?
2. Explain in two sentences why the basic interval can be a poor choice when the bootstrap distribution is right-skewed.
3. A Monte Carlo permutation test with B = 999 finds 0 permutations as extreme as the observed one. What p-value do you report?
4. Why can a seed reproduce your result in R but not make Python give the same numbers?
5. Your simulation gives 93.1% coverage for a nominal 95% interval with 2,000 runs. Is the interval clearly undercovering? Use the Monte Carlo SE.
6. True or false: a permutation test of a difference in means is valid whenever the observations are independent.

My notes and next steps

Sources and further reading

Structure aligned with the method-card shape of this library; content based on the original method papers and the public-domain NIST/SEMATECH e-Handbook, written in our own words.

• Efron, B. (1979). Bootstrap methods: another look at the jackknife. The Annals of Statistics, 7(1), 1–26.
• Efron, B. (1987). Better bootstrap confidence intervals. Journal of the American Statistical Association, 82(397), 171–185.
• Efron, B., and Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman and Hall.
• Davison, A. C., and Hinkley, D. V. (1997). Bootstrap Methods and Their Application. Cambridge University Press (the basis of the R boot package).
• Good, P. I. (2005). Permutation, Parametric and Bootstrap Tests of Hypotheses (3rd ed.). Springer.
• Phipson, B., and Smyth, G. K. (2010). Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, 9(1).
• NIST/SEMATECH e-Handbook of Statistical Methods, section 1.3.3.4 Bootstrap Plot, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• R package boot documentation (boot, boot.ci); SciPy documentation (scipy.stats.bootstrap, scipy.stats.permutation_test).
• Further reading: de Smith, M. J., Statistical Analysis Handbook (2024 edition, statsref.com); OpenIntro Statistics (CC BY-SA).

Answer key

Example B.

• H1: the format label is related to scores (the two population distributions differ in location; two-sided).
• $D_{\text{obs}} = 8.125$ points. $\binom{16}{8} = 12{,}870$ splits.
• 1. $0.01134421 \times 12{,}870 = 146$ splits were as or more extreme.
• 2. The exact test lists every split, so there is no randomness: both languages count the same 146 splits. The Monte Carlo versions draw different random splits (different generators), so their p-values differ by Monte Carlo error.
• 3. p = .011 < .05: reject H0. The data are unlikely under "format does not matter".
• 4. "Flipped-format students scored 8.1 points higher on average (M = 83.1 vs 75.0), exact permutation test p = .011, 12,870 splits." The Welch t-test agrees: t(14.0) = 3.12, p = .008, 95% CI [2.53, 13.72].

Exercise C.

R (seed 2026): d <- replicate(9999, mean(sample(flipped, replace = TRUE)) - mean(sample(lecture, replace = TRUE))), then quantile(d, c(0.025, 0.975)) gives [3.244, 12.625] with bootstrap SE 2.437. Python (stats.bootstrap with method "percentile", rng=np.random.default_rng(2026)) gives [3.250, 12.875] with SE 2.438. Both are narrower than the Welch interval [2.53, 13.72], partly because the percentile bootstrap ignores the extra uncertainty of estimating the SD from 8 values per group (the bootstrap uses divide-by-n variability, and there is no t-like widening). Both intervals exclude zero, in line with the permutation test.

Self-test.

1. Weekly sales in time order are autocorrelated; resampling single weeks breaks that dependence and understates the SE. Use a block (moving-block or stationary) bootstrap, or model the series (academic/stats-time-series-arima).
2. The basic interval reflects the bootstrap quantiles around the estimate. With a long right tail in the replicates, reflection puts the long side of the interval on the left. That is right only if the bootstrap error mimics the real sampling error (a pivot); when it does not, as with a small-sample median, the interval points the wrong way, and the percentile or BCa interval, which respect the shape of the replicates, are safer.
3. (0 + 1)/(999 + 1) = 0.001, so report p = .001 (or p ≤ .001), with B = 999 stated; increase B if the exact value matters.
4. A seed fixes the starting state of one generator. R and NumPy use different generators and different sampling algorithms, so the same seed produces different streams.
5. Monte Carlo SE ≈ $\sqrt{0.95 \times 0.05/2000} = 0.0049$. 0.931 is about 3.9 SE below 0.95, so yes, the undercoverage is clear (not just noise).
6. False. It also needs exchangeability under H0: the groups must have the same distribution when H0 is true, which random assignment guarantees but observational groups with different spreads may not.

Revision history