Resampling and Selection Stability

Resampling and Selection Stability#

Bootstrap-based stability answers a question worth asking: if I had collected this data again, would symbolic regression have picked the same terms? But the answer is only as honest as the unit you resample, and getting that wrong does not produce an obvious error. It produces a confident, narrow, plausible-looking spread that is far too narrow.

Your rows are…	Resample	Call
Independent measurements	rows	`bootstrap_model_selection(model, X, y)`
Grouped — several rows per condition, curve, subject, batch	groups	`bootstrap_model_selection(model, X, y, groups=...)`
Produced by an upstream fit — smoother, estimated derivative, simulation	the whole pipeline	`bootstrap_model_selection(model, None, None, resample_fn=...)`

The default is rows, which is right only for the first case.

Grouped rows#

If every row of an isotherm shares a temperature, a row bootstrap leaks: a replicate that drops half the isotherm still trains on the other half. The spread it reports reflects within-curve noise, not the between-condition variability you actually care about.

import numpy as np
from jaxsr import bootstrap_model_selection

# One label per row saying which condition / curve / subject it came from
groups = np.repeat([250.0, 275.0, 300.0, 325.0], 40)

stability = bootstrap_model_selection(model, X, y, n_bootstrap=200, groups=groups, seed=0)
stability["resampling"]        # 'groups'
stability["stability_score"]

Whole groups are drawn with replacement, so a replicate never sees part of a group it also trained on. At least two distinct groups are required.

The same concern applies to cross-validation, where a row-level split leaks a group across the train/test boundary and reports an optimistic score:

from jaxsr import cross_validate

scores = cross_validate(model, X, y, groups=groups, cv=4)
scores["mean_test_score"], scores["std_test_score"]

Passing groups promotes strategy to "group-kfold" automatically, so you cannot silently get a row-level split on grouped data. The available strategies are "kfold", "group-kfold" and "leave-one-group-out".

cv cannot exceed the number of distinct groups — with 4 conditions the default cv=5 raises rather than quietly producing an empty fold. Either lower cv, or use strategy="leave-one-group-out" when groups are few:

scores = cross_validate(model, X, y, groups=groups, strategy="leave-one-group-out")

Rows that came out of an upstream fit#

When the regression target is an estimated derivative, each row is a spline evaluation — not a measurement. Resampling those rows perturbs nothing about the smoother that produced them, so the dominant error source is invisible to a row-wise bootstrap.

resample_fn regenerates the data per replicate, including the upstream stage:

import numpy as np
from jaxsr import bootstrap_model_selection

def replicate(rng):
    """Called once per replicate with the bootstrap's RandomState."""
    raw = simulate_experiment(seed=rng.randint(2**31))   # or resample raw measurements
    smoother = fit_smoother(raw)                          # re-run the upstream stage
    return smoother.features(), smoother.derivative()     # -> (X_b, y_b)

stability = bootstrap_model_selection(
    model, None, None, n_bootstrap=90, seed=0, resample_fn=replicate
)
stability["resampling"]              # 'pipeline'
stability["n_distinct_structures"]

X and y go unused — pass None — because resample_fn supplies every replicate. Taking the bootstrap’s RandomState as its only argument is what keeps the ensemble reproducible under seed.

groups and resample_fn are mutually exclusive: resample_fn already decides how a replicate is built.

Reporting replicates you produced yourself#

If you orchestrated the replicates — a loop over simulated datasets, a batch of fits from a cluster job — hand them to the same reporting code:

from jaxsr import summarize_selection_replicates

models = [fit_one_replicate(i) for i in range(90)]
summary = summarize_selection_replicates(models, reference=full_fit, resampling="pipeline")

summary["stability_score"]
summary["n_distinct_structures"]
summary["feature_frequencies"]
summary["parameter_distributions"]   # for parametric bases

It also accepts mappings ({"features": [...], "expression": "..."}) or plain lists of selected term names, so replicates from outside JAXSR report identically. Without a reference, stability_score is the frequency of the most common structure.

Features are keyed by basis identity, not by rendered name, so a parametric basis appears once as the registered "exp(-a*x)" rather than once per refitted value. Before this, a parametric library could score 0.0 stability while every replicate selected the same basis.

Reading the result honestly#

A narrow coefficient interval is misleading when the bootstrap keeps selecting a different symbolic structure; that is why n_distinct_structures and structures are reported alongside feature_frequencies.

And stability is not validity. It measures whether the search is reproducible, not whether the model is right. In at least one well-documented case the two point in opposite directions: see the superposition guide, where a material with no valid transform produced a perfectly stable law while the genuine case disagreed with itself 38% of the time. Whenever a held-out test is available, that is the verdict; stability is context.