The mcycle dataset records the acceleration of a motorcycle helmet during a simulated crash. Variance in acceleration grows sharply during impact, making the data strongly heteroscedastic. In such cases, quantile regression is more informative than a mean fit: instead of a single conditional mean, you obtain a family of curves that reveal how the full distribution of the response shifts across time.
Fit
Load the dataset and fit a QuantileGAM with the default quantile set (10th, 25th, 50th, 75th, and 90th percentiles). The smooth term s(times, k=15) gives the model enough flexibility to capture the non-linear acceleration trajectory.
import whittaker as wk
# Load data and fit a five-quantile GAM
data = wk.load_dataset("mcycle")
model = wk.QuantileGAM("accel ~ s(times, k=15)").fit(data)
model.summary()
'QuantileGAM summary\n============================================================\nFormula: accel ~ s(times, k=15)\nQuantiles: [0.1, 0.25, 0.5, 0.75, 0.9]\nNon-crossing: True\nSigma: 0.1\n\n tau=0.10: edf=14.3, dev=70278231591.4\n tau=0.25: edf=14.3, dev=222720818495.8\n tau=0.50: edf=14.3, dev=296786551906.6\n tau=0.75: edf=14.3, dev=301880734857.5\n tau=0.90: edf=2.0, dev=63869001446.0'
Predict
Generate predictions on an evenly spaced grid across the observed time range. The result is a dictionary keyed by the exact quantile values; each entry holds a PredictionResult with a .values array.
import numpy as np
# Predict on a fine grid over observed time range
new_data = {"times": np.linspace(data["times"].min(), data["times"].max(), 200)}
preds = model.predict(new_data)
preds[0.5].values[:5]
array([-1.04636168e+09, -1.02097481e+09, -9.88019600e+08, -9.57536026e+08,
-9.31254712e+08])
The 10th and 90th percentile predictions mark the outer limits of the estimated distribution at each time point.
array([-1.07200688e+09, -1.02097481e+09, -9.88019600e+08, -9.57536026e+08,
-9.31254712e+08])
array([-2.76422515e+08, -2.55166163e+08, -2.34069269e+08, -2.09946139e+08,
-1.86002621e+08])
Interpret
Check empirical coverage, which is the fraction of training observations that fall between the outermost predicted quantiles.
The spread between the 10th and 90th quantile curves is widest during the high-impact phase of the crash, directly capturing the heteroscedastic nature of the data. Where the curves are tightly packed, the response is relatively predictable; where they fan out, uncertainty is high. This is information a single conditional mean cannot convey.