# Fit Quantile Regression

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](../reference/QuantileGAM.md#whittaker.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.


``` python
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](../reference/PredictionResult.md#whittaker.PredictionResult) with a `.values` array.


``` python
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.


``` python
preds[0.1].values[:5]
```


    array([-1.07200688e+09, -1.02097481e+09, -9.88019600e+08, -9.57536026e+08,
           -9.31254712e+08])


``` python
preds[0.9].values[:5]
```


    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.


``` python
model.coverage(data)
```


    0.9624060150375939


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.