2.6 Resampling Methods

Data have natural variability. Any statistic we compute from a sample (a mean, a correlation, a regression slope) inherits that variability.

Resampling methods let us measure it. Broadly, resampling refers to any technique where we repeatedly draw observations from a sample, recompute a statistic on each draw, and study the distribution of results. Applications include hypothesis testing, uncertainty propagation, and confidence intervals.

The lesson has three levels. Level 1 walks through three basic resampling techniques: randomization, bootstrapping, and Monte Carlo. Level 2 applies the bootstrap to model inference for a linear regression, first on synthetic GNSS data with a known answer, then on real GNSS data from the Cascadia subduction zone. Level 3 turns to the other meaning of resampling — changing the sampling of a signal: downsampling without aliasing, interpolating gaps under an explicit policy, aggregating irregular station networks, and the block bootstrap for correlated noise.

First, import the modules we need:

🖥️ Lecture slides — Session 06 (Mon Oct 12)

A note on random numbers. Older NumPy code seeds the global random state with np.random.seed(42) and then calls functions like np.random.normal. Modern NumPy replaces this with an explicit generator object: rng = np.random.default_rng(42). The generator carries its own state, so different parts of a program (or a notebook) do not interfere with each other through a hidden global. We use the generator idiom throughout this book. Fixing the seed makes the notebook reproducible from top to bottom.

1. Examples of resampling techniques (Level 1)¶

1.1 Randomization¶

Given two datasets, $A$ and $B$, and a parameter, $x$, we can randomly reassign observations to either $A$ or $B$, calculate some statistic (e.g., $x_A - x_B$), and repeat to build a distribution of that statistic under the null hypothesis that group labels do not matter.

Now, we resample.

For each run, we shuffle the combined data, split it into two new groups of the original sizes, and recalculate the difference of means.

1.2 Bootstrapping¶

When bootstrapping, you repeatedly draw observations, with replacement, from a sample. Each bootstrap resample has the same size as the original sample. Because you draw with replacement, a resample will contain duplicate observations and omit others. The spread of the statistic across resamples estimates the sampling uncertainty of that statistic.

Bootstrapping does not make strong assumptions about the underlying distribution of the data.

We first create a synthetic “true population” of correlated data using the generator method multivariate_normal (docs here). The mean of both variables is zero and their covariance is -0.75.

Check that the data is indeed anticorrelated by plotting one variable against the other.

We verify that the Pearson correlation coefficient is close to our target using the numpy function corrcoef.

In practice we never observe the full population. We now take a small subset of the data --- think of the subset as the sample we actually collected.

Now we bootstrap. Each resample draws len(subset) pairs from the sample itself, with replacement. The resample size matches the original sample size: that is what makes the spread of the resampled statistic mimic the sampling variability of the original estimate.

The median of the bootstrap estimates sits near the sample correlation coefficient, not the true population value. The bootstrap cannot fix the bias of a small sample; it quantifies the uncertainty around the estimate you have.

What happens if you increase the size of the original sample?

What if the resample size differs from the sample size?¶

The bootstrap is defined with resample size equal to the sample size. As an experiment, we deliberately break that rule and vary the resample size. Watch the spread of the bootstrap distribution.

Smaller resamples produce a wider spread of the statistic: a correlation estimated from 10 pairs is noisier than one estimated from 50. Resamples larger than the sample produce a spread that is too narrow, which understates the true uncertainty. Only a resample size equal to the original sample size reproduces the sampling variability of the estimate you actually made.

1.3 Monte Carlo¶

Named after the casino in Monaco, Monte Carlo methods involve simulating new data based on a known (or assumed!) statistical model. Unlike the previous two examples, we do not take draws from an existing sample.

Monte Carlo techniques have many applications: probabilistic risk assessment, uncertainty propagation in models, and evaluation of systems too complicated for closed-form analysis.

Here, we demonstrate Monte Carlo sampling by estimating π.

We begin with a central conceit: the ratio of the area of a circle to the area of its bounding square is $\pi/4$.

We then imagine a circle of radius 1 inscribed within a square with sides going from -1 to 1.

We now estimate π using the counts of points in the circle and in the square as approximations to their areas.

With a few runs, we do not get a good answer.

A Monte Carlo approach needs more runs to converge.

Explore how changing the number of runs changes your estimate of π (and how your calculation converges).

2. Using resampling for robust model inference (Level 2)¶

The plan: fit a linear trend to GNSS position data and use the bootstrap to put an uncertainty on the slope (the plate velocity). We start with a synthetic series, where we know the true velocity, so we can compare the bootstrap distribution against the right answer. Then we repeat the analysis on real data from station P395 in the Pacific Northwest.

2.1 Synthetic GNSS series with a known velocity¶

The helper mlgeo_synth.gnss_series builds a synthetic daily GNSS displacement series with a known ground truth: a linear trend, seasonal cycles, and realistic noise (white, flicker, and random walk). We set the true velocity to 12 mm/yr.

Fit a straight line with scipy.stats.linregress. The slope is our velocity estimate.

The point estimate is close to the truth, but how confident should we be? Bootstrap: resample the (time, displacement) pairs with replacement --- resample size equal to the data length --- refit the line each time, and collect the slopes.

Two things to notice.

First, the bootstrap distribution is centered on the estimated slope, not on the truth. The bootstrap quantifies the variability of the estimator; it cannot remove the offset between the estimate and the true velocity that this particular noise realization produced.

Second, the distribution is very narrow --- and the true velocity sits several bootstrap standard deviations away from its center. The error bar is too small. The reason: resampling pairs treats the residuals as independent, but GNSS noise is time-correlated (flicker and random walk). The pair bootstrap therefore understates the true uncertainty. More advanced schemes (block bootstrap) resample contiguous chunks of the series to preserve the correlation. Having the ground truth is what exposed this: with real data alone, the tight histogram would have looked reassuring.

2.2 Plate motion from real geodetic data¶

Now the real thing. We use a GNSS time series from station P395 in the Pacific Northwest and estimate the long-term motion due to the Cascadia subduction zone.

We download the time series from the University of Nevada, Reno data center. The tenv3 file is whitespace-delimited with one header line. The columns we need are the decimal year (yyyy.yyyy) and the east, north, and up positions in meters (__east(m), _north(m), ____up(m)).

2.3 Linear regression¶

There is a clean linear trend in the horizontal position data. We can fit the data using:

$E(t) = V_e t + u_e,$
$N(t) = V_n t + u_n,$

where $t$ is time. We regress the data to find the coefficients $V_e$, $u_e$, $V_n$, $u_n$. The displacements are mostly westward, so we focus on the East component $E(t)$ for this exercise. The coefficients $u_e$ and $u_n$ are the intercepts at $t=0$. They are not zero here because $t$ starts in 2006. The coefficients $V_e$ and $V_n$ have the dimension of velocities:

$V_e \sim E(t) / t$, $V_n \sim N(t) / t$,

so this example lets us discuss a simple linear regression and resampling. We use both a SciPy function and a scikit-learn function.

To measure fit performance, we measure how well the variance is reduced by fitting the data (scatter points) against the model. The variance is:

$\text{Var}(x) = 1/n \sum_{i=1}^n (x_i-\hat{x})^2$,

where $\hat{x}$ is the mean of $x$. When fitting the regression, we predict the values $x_{pred}$. The residuals are the differences between the data and the predicted values: $e = x - x_{pred}$. $R^2$ or coefficient of determination is:

$R^2 = 1 - \text{Var}(x-x_{pred}) /\text{Var}(x) = 1 - \text{Var}(e) /\text{Var}(x)$

The smaller the error, the “better” the fit (we will discuss later that a fit can be too good!), and the closer $R^2$ is to one.

We can also use the scikit-learn package:

To evaluate the errors of the model fit using sklearn, we use the following functions:

2.4 Bootstrapping the velocity¶

Now we use bootstrapping to estimate the slope of the regression over many resampled datasets, exactly as we did for the synthetic series.

Scikit-learn provides resample in the utils module. Make sure you use replace=True and a resample size equal to the data length (the default). For reproducible results, you can pass a fixed random_state. Bootstrapping is usually repeated many times (unlike K-fold cross-validation, the model-evaluation scheme of Chapter 3.8, which splits the data into a fixed number of non-overlapping folds).

The bootstrap spread on the real data is small because the trend dominates the noise, just as it did for the synthetic case. The same caveat applies: GNSS noise is time-correlated, so treat this error bar as a lower bound.

3. Signal resampling and irregular data (Level 3)¶

Sections 1 and 2 used “resampling” in the statistical sense: redrawing from a sample to measure uncertainty. The word has a second meaning that every sensor stream forces on you sooner or later: changing the sampling of a time series — downsampling a high-rate record, filling or refusing to fill gaps, and putting irregular observations onto a regular grid. Both meanings share a trap: done naively, they manufacture signal that was never measured.

This section works through three data streams that cover most of what you will meet in practice:

• a regular high-rate series (an hourly tide gauge) that we downsample to daily values;
• a regular series with outages (a daily GNSS record with gaps) that we interpolate under an explicit gap policy;
• an irregular sparse point stream (a multi-decade groundwater-well network) where nothing about the sampling is regular and aggregation choices dominate the result.

Each one is synthetic with known ground truth, so every repair gets graded. We close by returning to the too-narrow bootstrap error bar of section 2.1 and fixing it with a block bootstrap.

3.1 Downsampling and aliasing¶

Downsampling looks harmless: keep every $q$-th sample, discard the rest. It is not. A series sampled at interval $\Delta t$ can only represent frequencies up to the Nyquist frequency $f_N = 1/(2\Delta t)$. Any signal above the new Nyquist does not disappear when you subsample — it folds (aliases) into a lower frequency, masquerading as a signal that was never there.

The clean testbed is a tide gauge. mlgeo_synth.tide_gauge_series generates hourly sea level from four astronomical constituents plus a trend, a seasonal cycle, and weather noise — and returns the tidal components as ground truth. Suppose we want a daily sea-level series to study the slow (subtidal) signal.

The dominant constituent is M2, the principal lunar semidiurnal tide, with a period of 12.42 hours — a frequency of 1.93 cycles per day. A daily series has a Nyquist frequency of 0.5 cycles per day, so the entire tide lives above the new Nyquist. If we keep one sample per day (say, the midnight reading), M2 folds to $|1.93 - 2| = 0.068$ cycles per day: a spurious oscillation with a 14.8-day period and the full ~0.8 m tidal amplitude.

The fix is the rule every downsampling must follow: low-pass filter below the new Nyquist frequency first, then subsample. A daily mean is a crude low-pass filter (a 24-hour boxcar) and already suppresses most of the tide; scipy.signal.decimate applies a proper anti-alias filter before subsampling. We grade all three against the true tide-free daily sea level, which we can compute exactly because the generator returned the tide as a separate column.

The midnight samples carry a ~15-day oscillation of over half a meter that does not exist in the subtidal ocean — that is the aliased M2 tide, and its RMS error is thirty times larger than either anti-aliased version. Nothing about the naive series looks wrong; the fortnightly wiggle even resembles a plausible ocean signal. That is what makes aliasing dangerous: the artifact is physically dressed. Satellite altimetry lives with exactly this problem — the TOPEX/Jason orbit samples each point every ~10 days, aliasing M2 to a 62-day signal that must be modeled away.

pandas.DataFrame.resample('D').mean() — the one-liner you will reach for most often — is already a decent anti-alias filter for this purpose. The rule to internalize: before reducing the sampling rate, ask what lives above the new Nyquist frequency and remove it. If the answer is “nothing”, say so explicitly in your data card.

3.2 Gaps: interpolation is a decision, not a default¶

Real daily streams arrive with outages. mlgeo_synth.degrade_series injects gaps into the synthetic GNSS series from section 2.1 — plus one 150-day outage that we place, deliberately, across an earthquake (a station knocked out by the shaking it was supposed to record is not a hypothetical). The function returns the uncensored truth, so we can grade any repair.

The tempting one-liner is s.interpolate(): connect the dots across every gap. Before trusting it, measure what it costs — interpolate everything and compare against the truth, gap by gap.

Short gaps interpolate at or below the noise level: over a few days the trend and seasonal cycle barely move, so a straight line is as good as the data. The 150-day gap is different — its error is more than three times the noise, because the interpolation drew a smooth ramp through a 25 mm coseismic step it had no way of knowing about. The filled values are not noisy; they are confidently wrong, and any downstream code (a trend fit, an ML feature window) will treat them as measurements.

So we state a gap policy as an explicit decision rather than a library default:

• Interpolate gaps of 10 days or shorter. Rationale, from the table above: at 12 mm/yr and a ~3 mm seasonal amplitude, the deterministic signal moves well under the 2.4 mm noise level in 10 days, so interpolation error is bounded by noise. The threshold is set by this station’s signal rates — a station with faster motion or larger seasonal swings earns a shorter threshold, and the threshold belongs in the data card.
• Mask longer gaps as missing. A NaN is an honest statement: we do not know what happened in there — and in this record, something did.

The graded result: short-gap interpolation costs less than the noise, and the mask refuses to invent the 150 days we never saw. When this series later feeds a model, the NaNs force a documented choice (drop the window, flag it, impute with an uncertainty) instead of silently feeding fiction forward. That is the pattern for every repair in this section: fix what the data constrain, mask what they do not, and write the threshold down.

3.3 Irregular sparse points: the multi-well network¶

The third stream has no grid to start from. mlgeo_synth.well_table mimics forty years of water-level measurements across 25 monitoring wells: every well shares one regional signal (seasonal recharge on top of a slow decline), but each has its own datum offset (meters!), its own noise level, its own active period, a multi-year gap, and uneven visit dates. A few wells were only ever visited a handful of times. This is what operational hydrology, geotechnical monitoring, and legacy archives actually look like.

The goal: recover the regional head signal — the shared trend and seasonal cycle — as a quarterly series. The naive move is resample('QS').mean() over all observations. Watch what it does: whenever a well with a high datum offset enters or leaves the record (and they all enter and leave at different times), the mean jumps by a chunk of that offset. The “regional signal” it produces is mostly a history of which wells were being visited.

The gap-aware version makes two moves:

1. Work in anomalies. Subtract each well’s own mean first, so the meter-scale datum offsets cancel before any averaging. (This is exactly how global temperature series are built from weather stations.)
2. Weight by measurement quality. Each observation carries a reported sigma_m; weighting by $1/\sigma^2$ — the inverse-variance weight, which minimizes the variance of the combined estimate — keeps a steel-tape reading from diluting a transducer record.

Anomalies still leave a small bias — each well samples a different piece of the 40-year decline, so its own mean absorbs a slightly different trend segment — but that residual is decimeters, not the meters the offsets would inject. We grade both against truth["regional"], up to a constant (the datum is arbitrary, so we compare all series with their means removed).

The naive mean is wrong by a factor of a few — meter-scale jumps produced entirely by wells entering and leaving the record — while the anomaly-weighted series tracks the true regional decline and its seasonal cycle to about half a meter RMS — mostly the residual trend-segment bias noted above. Nothing sophisticated happened: the entire improvement came from refusing to average incompatible things. Before aggregating any irregular multi-site stream, ask what enters and leaves the average as the composition changes, and remove per-site levels first.

3.4 Closing the loop: the moving-block bootstrap¶

Section 2.1 ended with a warning: the pair bootstrap gave a velocity error bar so narrow that the true velocity sat far outside it, because resampling individual days destroys the time correlation of GNSS noise, and correlated noise is exactly what makes a trend uncertain. The fix is the moving-block bootstrap: instead of resampling days, resample contiguous blocks of residuals, so that each resample preserves the noise correlation up to the block length. We fit the line once, resample blocks of its residuals, add them back to the fitted line, and refit.

The block length is — again — a stated decision: it must exceed the correlation time of the noise you care about. We use 100 days, comfortably longer than the flicker-noise correlation at the periods that matter for a decade-long trend; you can check the sensitivity by rerunning with 50 or 200.

The block bootstrap widens the error bar by roughly a factor of seven — and now the true velocity sits about one standard deviation from the estimate, which is what an honest error bar looks like. Nothing about the data changed; only the resampling respected the correlation the noise actually has. The narrow pair-bootstrap histogram was not conservative, it was wrong — and on real data, with no ground truth to flag it, it would have been published.

The same caveat now transfers to the real P395 record of section 2.4: its pair-bootstrap error bar is a lower bound, and a block bootstrap on its residuals is the follow-up exercise.

4. Exercise¶

1. Alias hunting. Regenerate the tide gauge with n_days=365 and downsample by keeping the noon sample instead of midnight. Does the aliased period change? Explain why or why not from the folding formula.
2. Gap policy sensitivity. Rerun section 3.2 with max_gap_days of 3 and of 30. Report the RMS error and the number of filled samples for each. Where would you set the threshold for a station moving at 50 mm/yr, and why?
3. Well weighting. In section 3.3, drop the $1/\sigma^2$ weights (plain mean of anomalies). How much of the improvement over the naive mean survives? What does that tell you about which of the two moves (anomalies, weights) carries the load for this network?
4. Block length. Rerun the block bootstrap with block lengths of 10, 50, 200, and 500 days and plot the bootstrap standard deviation against block length. Explain the trend at both extremes.