Skip to content
Explore Automated Analysis

How to bootstrap a confidence interval for the median in SciPy (with Code)

Editorial · Published

The Short Answer (TL;DR)

Call scipy.stats.bootstrap with np.median, a one-element data tuple and a recorded random generator seed to estimate a median confidence interval. Set method="percentile" explicitly for this transparent example and assess whether independent resampling and its coverage approximation suit the data.

import numpy as np; from scipy import stats
x = np.array([4.2,4.9,5.1,4.7,5.8,6.3,5.5,4.6,5.7,6.0,5.2,4.4,5.9,6.2,5.0])
r = stats.bootstrap((x,),np.median,n_resamples=9999,batch=256,method="percentile",confidence_level=.95,rng=np.random.default_rng(20261027)); print(np.median(x),r.confidence_interval)

Detailed Implementation

This example uses synthetic data and does not need a network connection. Install the recorded package versions in a separate Python environment before running it.

python -m pip install numpy==2.5.4 scipy==1.18.1
from pathlib import Path
import numpy as np
from scipy import stats

# Fifteen independent synthetic observations in common measurement units.
x = np.array([4.2,4.9,5.1,4.7,5.8,6.3,5.5,4.6,5.7,6.0,5.2,4.4,5.9,6.2,5.0])
if x.ndim != 1 or x.size < 3 or not np.isfinite(x).all():
    raise ValueError("Use at least three finite one-dimensional observations.")
if np.ptp(x) == 0:
    raise ValueError("A constant sample cannot estimate sampling uncertainty with this bootstrap.")
seed, n_resamples = 20261027, 9999
result = stats.bootstrap(
    (x,), np.median, vectorized=True, axis=0, n_resamples=n_resamples,
    batch=256, confidence_level=.95, alternative="two-sided",
    method="percentile", rng=np.random.default_rng(seed),
)
low, high = result.confidence_interval
if not np.isfinite([low,high,result.standard_error]).all() or high <= low:
    raise ValueError("The resampled statistic is degenerate; inspect the sample and interval method.")
report = (
    f"n={x.size}; median={np.median(x):.6f} measurement units\n"
    f"95% percentile CI=[{low:.6f}, {high:.6f}] measurement units\n"
    f"bootstrap SE={result.standard_error:.9f}; resamples={n_resamples}; seed={seed}\n"
)
Path("bootstrap-median.txt").write_text(report, encoding="utf-8")
print(report)

The sample median is 5.200000 and the 95% percentile interval is [4.890000, 5.800000] measurement units; the output also records the bootstrap standard error and seed.

Why is the data wrapped as (x,)?

The function expects one array per sample; the comma creates a one-element tuple. Each resample draws 15 values with replacement from these same 15 observations, then computes a median. The seed, resample count, interval method and package versions are part of the reported calculation, and the batch size limits temporary work arrays.

Why choose the percentile method here?

Its endpoints are directly inspectable quantiles of the resampled medians, making it useful for explaining and independently checking the computation. It is an explicit teaching choice, not a claim that it has the best coverage. SciPy defaults to BCa, which applies bias and acceleration adjustments. A median with a small or tied sample can have a coarse bootstrap distribution; choosing another method does not automatically repair that problem.

What is independently verified?

A separate NumPy index-resampling calculation reproduces all 9,999 medians with the recorded seed. Its 2.5th and 97.5th percentiles give [4.89, 5.80], and its sample standard deviation agrees with the returned standard error. Reusing the saved bootstrap result allows a 99% interval without drawing additional resamples; a repeated seeded run also agrees exactly in the recorded environment.

Which sampling assumptions and interpretation matter?

Ordinary observation resampling assumes independent, identically distributed sampling from a population represented by the sample. A time series, clustered observations or a complex survey requires a resampling design that respects that structure. Fifteen synthetic values illustrate mechanics rather than establish coverage for a real study. Increasing the number of resamples reduces simulation variability but does not add observations or remove sampling bias. The interval concerns the population median, not the range containing 95% of future observations, and it is not a posterior probability statement.

Mathematical Context

For each bootstrap sample b, draw n indices independently and uniformly with replacement from 1 through n and calculate the sample median theta*_b. The percentile interval uses empirical quantiles Q_0.025 and Q_0.975 of these medians. The bootstrap standard error is their sample standard deviation with denominator B-1. The point estimate 5.2, interval endpoints 4.89 and 5.80, and standard error all use the original measurement units. Approximate 95% coverage describes repeated sampling behavior of the procedure under suitable assumptions; it is not guaranteed for this small sample or for every underlying distribution.

\[ \mathrm{CI}_{0.95}=\left[Q_{0.025}(\theta^*),Q_{0.975}(\theta^*)\right] \]

Reproducibility

Executed with Python 3.12.14; numpy 2.5.4, scipy 1.18.1. Automated checks confirm the stated output properties.

Primary documentation

Related guides