# 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 | Field | Value | |---|---| | Course or project | [STAT 2xx Computational statistics / thesis chapter 4] | | Learner or analyst | [Your name] | | Data set | [Name, source, n, date extracted] *Example:* Lakemoor University cafe wait times, n = 15 (fictional) | | Software and versions | R 4.5.1 with boot 1.3.31; Python 3.14 with numpy 2.5.4, scipy 1.18.1 | | Seeds and resamples | [Seed, B] *Example:* seed 2026, B = 9,999 | | Date | YYYY-MM-DD | | Version | [0.1] | ## 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). | Use when | Do not use when | |---|---| | The statistic has no convenient SE formula (median, trimmed mean, ratio, correlation, a difference in medians) | The sample is very small (n below about 10 to 15): the bootstrap only reshuffles the few values you have | | You want an interval that can be asymmetric for a skewed statistic | Observations are dependent (repeated measures, clusters, time series) and you would resample single rows | | Two groups were randomly assigned and you want a p-value that does not assume normality (permutation) | The question is about extremes (the maximum, a 99th percentile in a small sample) where the bootstrap is known to fail | | You want to check coverage or Type I error of a method under realistic data (Monte Carlo) | You need to fix bias from a bad design or a non-random sample: resampling cannot create information that was never collected | **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. | Variable | Role | Measurement level | Unit or coding | |---|---|---|---| | wait | Outcome (Example A) | Ratio, continuous | Minutes | | score | Outcome (Example B) | Interval, treated as continuous | Quiz points out of 100 | | format | Group (Example B) | Nominal, two levels | flipped, lecture | Example A data (n = 15 cafe visits, one per student, fictional): | Visit | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | |---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---| | wait (min) | 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 | 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): | Format | Scores | Mean | |---|---|---| | flipped | 78, 85, 81, 90, 74, 88, 83, 86 | 83.125 | | lecture | 72, 80, 69, 77, 84, 70, 75, 73 | 75.000 | ## Assumptions and how to check each | Assumption | Check | R | Python | What to do if violated | |---|---|---|---|---| | Observations are independent and identically distributed (bootstrap) | Design: one row per independent unit? Any clusters, repeated measures or time order? | `table(id)` to find repeated units; `acf(x)` if time ordered | `df["id"].value_counts()`; `statsmodels.graphics.tsaplots.plot_acf` | Resample clusters (cluster bootstrap) or blocks of consecutive values (block bootstrap) | | Sample is large enough to stand in for the population | n, number of distinct values, plot of the bootstrap distribution | `length(unique(x))`; `hist(b$t)` | `np.unique(x).size`; histogram of `res.bootstrap_distribution` | Report the limitation; prefer a method with known small-sample behaviour; do not trust BCa when n is tiny | | Statistic is smooth (bootstrap) | Is the statistic a median, quantile or maximum? Is the bootstrap distribution lumpy? | `table(b$t)` | `np.unique(res.bootstrap_distribution, return_counts=True)` | Use more data, a smoothed bootstrap, or report that the interval is coarse | | Exchangeability under H0 (permutation) | Were labels randomly assigned, or are groups exchangeable when H0 holds (same distribution, not only same mean)? | Design check; compare spreads with `var()` | `np.var(x, ddof=1)` | If spreads differ a lot, a permutation test of means can mislead: use Welch t or a studentised permutation statistic | | Enough resamples | Monte Carlo error of the p-value or interval endpoints | Rerun with another seed | Rerun with another `rng` | Increase B (9,999 or more for p-values near α; 10,000 or more for BCa) | ## 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. $$ \text{Percentile: } \left[\hat\theta^*_{(\alpha/2)},\ \hat\theta^*_{(1-\alpha/2)}\right] \qquad \text{Basic: } \left[2\hat\theta-\hat\theta^*_{(1-\alpha/2)},\ 2\hat\theta-\hat\theta^*_{(\alpha/2)}\right] $$ 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: $$ z_0=\Phi^{-1}\!\left(\frac{\#\{\hat\theta^*_b<\hat\theta\}}{B}\right),\qquad \alpha_1=\Phi\!\left(z_0+\frac{z_0+z_{\alpha/2}}{1-a\,(z_0+z_{\alpha/2})}\right) $$ 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. $$ p_{\text{exact}}=\frac{\#\{\text{splits with } |D_\pi|\ge |D_{\text{obs}}|\}}{\binom{n_1+n_2}{n_1}}, \qquad p_{\text{MC}}=\frac{1+\#\{b: |D_b|\ge |D_{\text{obs}}|\}}{1+B} $$ 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. | Symbol | Meaning | |---|---| | $n$ | Sample size (number of independent resampling units) | | $B$ | Number of bootstrap resamples or random permutations | | $\hat\theta$ | Statistic computed on the original sample | | $\hat\theta^*_b$ | Statistic computed on resample $b$ | | $\hat\theta^*_{(q)}$ | $q$ quantile of the $B$ replicates | | $z_0$, $a$ | BCa bias correction and acceleration | | $\Phi$, $z_q$ | Standard normal CDF and its $q$ quantile | | $D_{\text{obs}}$, $D_\pi$ | Observed difference, difference under permutation $\pi$ | ```mermaid flowchart TD Q["What do you need?"] --> I["An interval or SE for a statistic"] Q --> P["A p-value comparing groups"] Q --> C["To check a method on data like mine"] I --> D1["Units independent?"] D1 -->|No| BL["Cluster or block bootstrap"] D1 -->|Yes| N1["n at least about 15 and statistic smooth?"] N1 -->|No| W["Report limitation or use exact or Bayesian method"] N1 -->|Yes| BCA["BCa interval, B at least 9999, check percentile too"] P --> X["Labels randomly assigned or exchangeable under H0?"] X -->|No| WT["Welch t or studentised permutation"] X -->|Yes| SZ["Number of splits small enough to list?"] SZ -->|Yes| EX["Exact permutation test"] SZ -->|No| MC["Monte Carlo permutation, B at least 9999"] C --> SIM["Simulate from a known truth, record coverage or Type I error"] ``` **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")`): ```r 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): ```text 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: ```python 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): ```text 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: ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Bootstrap distribution of the median wait, 9999 resamples, R seed 2026", "width": 400, "height": 220, "data": {"values": [ {"median": 2.1, "count": 1}, {"median": 2.2, "count": 38}, {"median": 2.6, "count": 207}, {"median": 2.9, "count": 641}, {"median": 3.0, "count": 1287}, {"median": 3.3, "count": 1820}, {"median": 3.4, "count": 2090}, {"median": 3.7, "count": 1820}, {"median": 4.2, "count": 1247}, {"median": 4.8, "count": 613}, {"median": 5.1, "count": 192}, {"median": 7.9, "count": 41}, {"median": 9.4, "count": 2} ]}, "mark": "bar", "encoding": { "x": {"field": "median", "type": "ordinal", "title": "Bootstrap median (minutes)"}, "y": {"field": "count", "type": "quantitative", "title": "Number of resamples"} } } ``` **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: ```r 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): ```text 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: ```python 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): ```text 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 | Output line | Meaning | What to report | |---|---|---| | `original` (R `boot`) | Statistic on the original sample, $\hat\theta$ | Yes, as the point estimate | | `bias` (R `boot`) | Mean of replicates minus $\hat\theta$ | Usually no; a large bias relative to SE is a warning sign | | `std. error` (R) / `standard_error` (SciPy) | SD of the bootstrap replicates | Yes, as the bootstrap SE if you report an SE | | `Based on 9999 bootstrap replicates` | B | Yes, always with the seed and software | | `Percentile`, `Basic`, `BCa` columns | Three intervals from the same replicates | Report one, chosen in advance (usually BCa), not the one you like best | | `confidence_interval(low, high)` (SciPy) | The interval for the chosen `method` | Yes | | Warning "extreme order statistics used as endpoints" (R) | B is too small for BCa at that level | Increase B before reporting | | `null_distribution.size` (SciPy) | Number of permutations used; equals all splits when exact | Yes: say "exact" or the number of random permutations | | `pvalue` (permutation) | Proportion of permutations as or more extreme | Yes, three decimals | ## 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])." | Statistic | Estimate | Interval type | 95% CI | B | Seed | Software | |---|---|---|---|---|---|---| | *Example:* Median wait (min) | 3.4 | BCa | [2.6, 4.8] | 9,999 | 2026 | R 4.5.1, boot 1.3.31 | | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | **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. ```r 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): ```text coverage t interval: 0.9205 coverage percentile bootstrap: 0.894 Welch Type I error rate at alpha = .05: 0.0343 ``` ```python 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): ```text 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 | Item | Your entry | |---|---| | Confidence with this method (1 to 5) | [ ] | | What I can now do without notes | [ ] | | Open questions | [ ] | | Mistakes to copy to my error log (academic/stats-error-log) | [ ] | | Next review date (spacing: 2 days, 1 week, 3 weeks) | YYYY-MM-DD | | Next template to study | [academic/stats-time-series-arima or academic/stats-power-sample-size] | ## 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 | Version | Date | Author | Change | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Your name] | First version from the template | | [ ] | YYYY-MM-DD | [ ] | [ ] |