Activity 5.3

This page is an ipython notebook, meaning you can download it and run the executable cells. It directly follows Activity 5.3. Reference answers appear after the questions, challenges, predictions, and reflections. They show one reasonable solution or interpretation, but other correct approaches may exist.

Bootstrap & Permutation Procedures from Scratch¶

Tools: Python, NumPy, pandas, matplotlib, SciPy, Jupyter/VS Code

Learning objectives¶

By the end of this tutorial, you should be able to:

• generate bootstrap samples manually;

• construct bootstrap distributions;

• estimate bootstrap standard errors;

• construct percentile bootstrap confidence intervals;

• bootstrap means, medians, and paired statistics;

• preserve observational units during resampling;

• create permutation null distributions;

• calculate one- and two-sided permutation p-values;

• distinguish exact and Monte Carlo permutation logic;

• work with paired resampling structures;

• write reusable resampling functions;

• compare from-scratch results with SciPy.

Part 0 --- Set up¶

Create:

notebooks/week11_resampling.ipynb

python import numpy as np import pandas as pd import matplotlib.pyplot as plt

python rng = np.random.default_rng(556)

python N_BOOT = 10000 N_PERM = 10000

Part 1 --- A small observed sample¶

Create:

python x = np.array([ 4.2, 5.1, 3.8, 6.0, 4.9, 5.4, 4.5, 5.7 ])

python x.mean() np.median(x) x.std(ddof=1)

Part 2 --- One bootstrap sample¶

Run:

python bootstrap_sample = rng.choice( x, size=len(x), replace=True )

python bootstrap_sample

Reference answer¶

The nonparametric bootstrap treats the observed sample as an empirical approximation to the unknown population. We know the observed values but not the true population distribution that generated them.

Part 3 --- Why replacement matters¶

Try:

python rng.choice( x, size=len(x), replace=False )

Question¶

What has this actually done?

Why does it fail to create meaningful bootstrap variation in the sample values?

Reference answer¶

A bootstrap sample has the same size as the original and is drawn with replacement, so observations can repeat and others can be omitted. Replacement creates bootstrap variation.

Part 4 --- One bootstrap statistic¶

Calculate:

python bootstrap_mean = ( bootstrap_sample.mean() )

python x.mean() bootstrap_mean

Reference answer¶

Sampling the whole dataset without replacement only reorders the observations. Statistics such as the mean are unchanged, so it does not approximate sampling variability.

Part 5 --- Bootstrap distribution with a loop¶

Create:

python bootstrap_means = np.empty( N_BOOT )

python for b in range(N_BOOT): sample_b = rng.choice( x, size=len(x), replace=True )

bootstrap_means[b] = (
    sample_b.mean()
)

python bootstrap_means[:10]

Part 6 --- Visualize the bootstrap distribution¶

Plot:

python plt.hist( bootstrap_means, bins=30 )

plt.axvline( x.mean() )

plt.xlabel(“Bootstrap mean”) plt.ylabel(“Frequency”) plt.show()

Part 7 --- Bootstrap standard error¶

Calculate:

python bootstrap_se = ( bootstrap_means.std( ddof=1 ) )

python analytical_se = ( x.std(ddof=1) / np.sqrt(len(x)) )

python bootstrap_se analytical_se

Reference answer¶

The bootstrap distribution is a distribution of a statistic across resamples, not a distribution of raw observations. Its standard deviation estimates the statistic’s sampling standard error.

Part 8 --- Percentile confidence interval¶

Compute:

python bootstrap_ci = np.quantile( bootstrap_means, [0.025, 0.975] )

bootstrap_ci

Part 9 --- Bootstrap bias¶

Calculate:

python bootstrap_bias = ( bootstrap_means.mean() - x.mean() )

bootstrap_bias

Reference answer¶

A percentile interval takes empirical quantiles of the bootstrap statistic distribution, such as the 2.5th and 97.5th percentiles for a nominal 95% interval.

Part 10 --- Vectorized bootstrap¶

Generate:

python bootstrap_samples = rng.choice( x, size=( N_BOOT, len(x) ), replace=True )

python bootstrap_samples.shape

python bootstrap_means_vec = ( bootstrap_samples.mean( axis=1 ) )

python bootstrap_means_vec.mean() bootstrap_means_vec.std(ddof=1)

Reference answer¶

Bootstrap bias is estimated as the mean bootstrap statistic minus the original-sample statistic. Its sign gives the estimated direction of bias.

Part 11 --- Memory discussion¶

Estimate how many values would be created by:

text 1,000,000 bootstrap samples × 10,000 observations

Part 12 --- Bootstrap the median¶

Construct:

python bootstrap_medians = np.empty( N_BOOT )

python median_se = ( bootstrap_medians.std( ddof=1 ) )

median_ci = np.quantile( bootstrap_medians, [0.025, 0.975] )

Reference answer¶

Full vectorization can require a huge resample matrix. A loop or batched implementation may use much less memory and can therefore be preferable for large problems.

Part 13 --- Write a generic bootstrap function¶

Write:

python def bootstrap_statistic( data, statistic, n_boot, rng ): results = np.empty( n_boot )

for b in range(n_boot):
    sample = rng.choice(
        data,
        size=len(data),
        replace=True
    )

    results[b] = statistic(
        sample
    )

return results

python means = bootstrap_statistic( x, np.mean, N_BOOT, rng )

python medians = bootstrap_statistic( x, np.median, N_BOOT, rng )

Reference answer¶

Bootstrap is attractive when an analytic sampling distribution or SE formula is unavailable or inconvenient. A generic bootstrap function should accept data, a statistic, number of resamples, and an RNG.

Part 14 --- Bootstrap rows, not columns¶

Create:

python df = pd.DataFrame({ “age”: [22, 25, 31, 35, 41, 48], “score”: [70, 75, 80, 86, 91, 94] })

python observed_r = ( df[“age”] .corr( df[“score”] ) )

Part 15 --- Bootstrap paired rows¶

Generate bootstrap indices:

python indices = rng.choice( len(df), size=len(df), replace=True )

python sample_b = df.iloc[ indices ]

python sample_b[ “age” ].corr( sample_b[“score”] )

Part 16 --- Bootstrap a correlation¶

Repeat row-resampling N_BOOT times.

Store:

python bootstrap_correlations

Reference answer¶

For correlation, resample row indices so paired (x,y) observations stay together. Resampling columns independently would destroy the joint structure.

Part 17 --- Two-group data¶

Create:

python group_a = np.array([ 82, 88, 91, 85, 90, 87, 93, 89 ])

group_b = np.array([ 75, 80, 78, 84, 79, 77, 81, 76 ])

python observed_difference = ( group_a.mean() - group_b.mean() )

observed_difference

Part 18 --- State the null hypothesis¶

Write in words:

text H0:

Reference answer¶

The null hypothesis is that group membership is not associated with the outcome/statistic in the relevant sense; under the null, group labels are exchangeable.

Part 19 --- Pool the observations¶

Create:

python combined = np.concatenate([ group_a, group_b ])

python combined len(combined)

python n_a = len(group_a)

Reference answer¶

If H0 is true, the observed labels carry no systematic outcome information. Pooling and randomly reallocating labels simulates that null world.

Part 20 --- One permutation¶

Run:

python permuted = rng.permutation( combined )

python a_perm = permuted[ :n_a ]

b_perm = permuted[ n_a: ]

python perm_difference = ( a_perm.mean() - b_perm.mean() )

Part 21 --- Build the permutation distribution¶

Create:

python null_distribution = np.empty( N_PERM )

python for b in range(N_PERM): permuted = rng.permutation( combined )

a_perm = permuted[
    :n_a
]

b_perm = permuted[
    n_a:
]

null_distribution[b] = (
    a_perm.mean()
    - b_perm.mean()
)

Reference answer¶

Compare the observed statistic with the permutation distribution. An observed value deep in a tail is unusual under the null.

Part 22 --- Visualize the null distribution¶

Plot:

python plt.hist( null_distribution, bins=30 )

plt.axvline( observed_difference )

plt.xlabel( “Difference in means under H0” )

plt.ylabel(“Frequency”) plt.show()

Part 23 --- One-sided p-value¶

Suppose:

text H1: mean_A > mean_B

python p_greater = ( 1 + np.sum( null_distribution >= observed_difference ) ) / ( N_PERM + 1 )

p_greater

Reference answer¶

A one-sided p-value is the fraction of null statistics at least as extreme as the observed statistic in the prespecified direction.

Part 24 --- Two-sided p-value¶

Calculate:

python p_two_sided = ( 1 + np.sum( np.abs( null_distribution ) >= np.abs( observed_difference ) ) ) / ( N_PERM + 1 )

p_two_sided

Reference answer¶

A two-sided p-value counts null statistics whose magnitude is at least as large as the observed statistic’s magnitude (or an equivalent valid two-tail definition).

Part 25 --- Why add one?¶

Compare:

python np.mean( null_distribution >= observed_difference )

Reference answer¶

‘At least as extreme’ means at least as incompatible with H0 according to the chosen statistic and alternative. Adding one avoids reporting an estimated Monte Carlo p-value of exactly zero.

Part 26 --- Change the number of permutations¶

Repeat the p-value estimation using:

text 100 1,000 10,000

Reference answer¶

A finite permutation count introduces Monte Carlo error into the p-value estimate. More permutations make the estimate more stable.

Part 27 --- Write a statistic function¶

Create:

python def difference_in_means( x, y ): return ( np.mean(x) - np.mean(y) )

python observed_difference

text statistic

text resampling mechanism

Reference answer¶

A statistic function isolates the quantity of interest so the same permutation engine can be reused with means, medians, correlations, and other statistics.

Part 28 --- Generic independent permutation function¶

Write:

python def permutation_distribution( x, y, statistic, n_perm, rng ): combined = np.concatenate([ x, y ])

n_x = len(x)

results = np.empty(
    n_perm
)

for b in range(n_perm):
    permuted = rng.permutation(
        combined
    )

    x_perm = permuted[
        :n_x
    ]

    y_perm = permuted[
        n_x:
    ]

    results[b] = statistic(
        x_perm,
        y_perm
    )

return results

Reference answer¶

A generic independent permutation test pools outcomes, reallocates them into groups of the original sizes, recomputes the statistic, repeats, and compares the observed statistic with the resulting null distribution.

Part 29 --- Compare statistics¶

Repeat the permutation procedure using:

python def difference_in_medians( x, y ): return ( np.median(x) - np.median(y) )

Reference answer¶

Changing the statistic changes what aspect of the groups is being tested, so the null distribution and p-value can change even with the same data.

Part 30 --- Small exact permutation example¶

Use:

python x_small = np.array([ 1, 2 ])

y_small = np.array([ 4, 5 ])

text 1, 2, 4, 5

text C(4,2)

Reference answer¶

For tiny samples, enumerate every distinct label allocation to obtain an exact permutation distribution rather than relying on random permutations.

Part 31 --- Exact vs. Monte Carlo¶

Compare:

text all possible assignments

text randomly sampled assignments

Reference answer¶

Exact permutation evaluates all possible allocations and has no Monte Carlo error; Monte Carlo permutation samples a feasible subset of allocations.

Part 32 --- Paired data¶

Create:

python before = np.array([ 80, 72, 88, 91, 77, 85 ])

after = np.array([ 84, 75, 90, 94, 76, 89 ])

python differences = ( after - before )

python differences differences.mean()

Reference answer¶

Paired data must preserve the pairing. Use a paired null operation such as random sign flips of within-pair differences rather than pooling all observations.

Part 33 --- Preserve the pairing¶

Incorrect¶

Pool all before/after observations and randomly split them.

Better for a paired sign-randomization logic¶

For each difference, randomly choose:

text +d_i or -d_i

python signs = rng.choice( [-1, 1], size=len(differences) )

permuted_difference = ( differences * signs )

Part 34 --- Paired null distribution¶

Repeat the random sign flip N_PERM times.

Store mean differences.

Compare the observed paired mean difference with the resulting null distribution.

Question¶

Why does this preserve more of the original study design than pooling all observations?

Reference answer¶

The paired procedure preserves subject-level dependence and therefore respects the study design more closely than an independent-label permutation.

Part 35 --- Validate with SciPy: permutation¶

Import:

python from scipy.stats import ( permutation_test )

python def statistic( x, y ): return ( np.mean(x) - np.mean(y) )

python result = permutation_test( ( group_a, group_b ), statistic, n_resamples=N_PERM, alternative=“two-sided”, rng=np.random.default_rng( 556 ) )

python result.statistic result.pvalue

Part 36 --- Validate with SciPy: bootstrap¶

Import:

python from scipy.stats import bootstrap

Reference answer¶

SciPy should give a similar bootstrap result when statistic, confidence level, method, and randomization choices are aligned. Small numerical differences can arise from implementation details.

Part 37 --- Bootstrap vs. permutation challenge¶

For each goal, choose bootstrap or permutation.

A¶

Estimate the standard error of the median.

B¶

Test whether two randomized groups differ in mean outcome under a no-group-effect null.

C¶

Estimate uncertainty in a regression coefficient.

D¶

Construct a null distribution for a difference in groups.

E¶

Estimate a confidence interval for a correlation.

Explain each choice.

Reference answer¶

Bootstrap estimates sampling uncertainty around an observed statistic; permutation builds a null distribution under an exchangeability assumption for hypothesis testing.

Part 38 --- Mini resampling study¶

Create:

python rng = np.random.default_rng(556)

control = rng.normal( 70, 10, size=30 )

treatment = rng.normal( 76, 10, size=30 )

Part 39 --- Reproducibility and Monte Carlo stability¶

Repeat the mini study using:

text 1,000 resamples 10,000 resamples

Part 40 --- Resampling audit¶

Before trusting a procedure, answer:

text What is the observational unit? What statistic is being computed? What is resampled? Is replacement used? What null assumption justifies permutation? Is pairing/clustering preserved? How many resamples are used? How is randomness controlled? What empirical distribution is being constructed? What inferential quantity is extracted?

Reference answer¶

A resampling audit should state the observational unit, statistic, what is resampled, replacement rule, null hypothesis if relevant, number of resamples, and RNG/seed strategy.

Part 41 --- Git checkpoint¶

Run:

bash git status

bash git add .

bash git commit -m “Complete Week 11 bootstrap and permutation exercises”

bash git push

Part 42 --- Final reflection¶

Answer in Markdown.

1. Bootstrap¶

Why does a bootstrap sample use replacement?

2. Empirical population¶

What population does the nonparametric bootstrap sample from?

3. Bootstrap distribution¶

What does it approximate?

4. Standard error¶

How do you estimate a bootstrap SE?

5. Percentile interval¶

How is it constructed?

6. Permutation¶

What does a permutation null distribution represent?

7. Exchangeability¶

Why is it central to a permutation test?

8. p-value¶

How is a Monte Carlo permutation p-value calculated?

9. Pairing¶

Why can we not arbitrarily shuffle observations in a paired study?

10. Comparison¶

What is the most important conceptual difference between bootstrap and permutation resampling?

Completion checklist¶

• Created Week 11 notebook

• Generated one bootstrap sample

• Explained why replacement is required

• Constructed a bootstrap distribution with a loop

• Visualized the bootstrap distribution

• Calculated bootstrap SE

• Compared bootstrap and analytical SE

• Constructed a percentile CI

• Estimated bootstrap bias

• Implemented vectorized bootstrap sampling

• Discussed memory tradeoffs

• Bootstrapped a median

• Wrote a generic bootstrap function

• Preserved row pairings when resampling

• Bootstrapped a correlation

• Calculated an observed group difference

• Constructed one random permutation

• Built a permutation null distribution

• Calculated one-sided p-value

• Calculated two-sided p-value

• Used a finite-resampling correction

• Examined permutation-count stability

• Wrote reusable statistic functions

• Wrote a generic permutation function

• Compared mean and median test statistics

• Constructed a small exact permutation distribution

• Compared exact and Monte Carlo permutation tests

• Preserved paired design using sign randomization

• Validated permutation output with SciPy

• Validated bootstrap output with SciPy

• Completed the bootstrap-vs-permutation challenge

• Completed the mini resampling study

• Compared results across resampling counts

• Audited the resampling design

• Committed Week 11 work to Git

• Pushed work to GitHub

What you should now understand¶

text Observed data ↓ Choose statistic ↓ Bootstrap sample with replacement ↓ sampling uncertainty ↓ SE / confidence interval

or

Permutation rearrange under H0 ↓ null distribution ↓ p-value

Next week we move into **numerical optimization**, including root
finding, Newton--Raphson, gradient descent, and statistical objective
functions.

Reference answer¶

Bootstrap samples with replacement from the empirical distribution to approximate sampling variability. Permutation rearranges labels/signs under a null. Both require resampling the correct observational unit.