Computing prediction intervals using quantile regression

Environment setup

We need to install some extra dependencies for this notebook if needed (when running jupyterlite).

%pip install -q https://pypi.anaconda.org/ogrisel/simple/polars/1.24.0/polars-1.24.0-cp39-abi3-emscripten_3_1_58_wasm32.whl
%pip install -q skrub altair holidays plotly nbformat

ERROR: polars-1.24.0-cp39-abi3-emscripten_3_1_58_wasm32.whl is not a supported wheel on this platform.

Note: you may need to restart the kernel to use updated packages.

Note: you may need to restart the kernel to use updated packages.

import functools
import importlib
import warnings

import altair
import numpy as np
import polars as pl
import skrub

import tutorial_helpers

importlib.reload(tutorial_helpers)

from feature_engineering_lib import feature_engineering_outputs, time_range
from multiple_horizons_prediction_lib import neg_mape
from single_horizon_prediction_lib import TimeSeriesSplitter, get_cv_results, get_hgb_params
from tutorial_helpers import (
    binned_coverage,
    collect_cv_predictions,
    plot_lorenz_curve,
    plot_predictions,
    plot_reliability_diagram,
    plot_residuals_vs_predicted,
)

# Ignore warnings from pkg_resources triggered by Python 3.13's multiprocessing.
warnings.filterwarnings("ignore", category=UserWarning, module="pkg_resources")

Define the quantile regressors

In this section, we show how one can use a gradient boosting but modify the loss function to predict different quantiles and thus obtain an uncertainty quantification of the predictions.

In terms of evaluation, we reuse the MAPE score. However, they it is not helpful to assess the reliability of quantile models. For this purpose, we use a derivate of the metric minimized by the quantile regressors: the pinball loss. We use the D2 score that is easier to interpret since the best possible score is bounded by 1 and a score of 0 corresponds to constant predictions at the target quantile.

from sklearn.metrics import d2_pinball_score, mean_absolute_percentage_error


def split_by_quantile(pred):
    quantile_cols = {}
    for c in pred.columns:
        quantile_cols.setdefault(c.split("__")[1], []).append(c)
    return {
        q: pred.select(cols).rename(lambda c: c.split("__")[0]) for q, cols in quantile_cols.items()
    }


def neg_mape_median(y_true, y_pred):
    "return the negative MAPE for the prediction for the quantile 0.5"
    return neg_mape(y_true, split_by_quantile(y_pred)["q_0.5"])


def neg_mape_median_scorer(estimator, X, y):
    return neg_mape_median(y, estimator.predict(X))


def pinball(y_true, y_pred):
    quantile_predictions = split_by_quantile(y_pred)
    scores = {}
    for q, q_pred in quantile_predictions.items():
        scores[f"d2_pinball_score__average__{q}"] = d2_pinball_score(
            y_true, q_pred, alpha=float(q.removeprefix("q_"))
        )
        detail = d2_pinball_score(y_true, q_pred, multioutput="raw_values")
        scores.update(
            {f"d2_pinball_score__{c}__{q}": float(s) for c, s in zip(y_true.columns, detail)}
        )
    return scores


def pinball_scorer(estimator, X, y):
    return pinball(y, estimator.predict(X))

TIME_HORIZONS = (1, 12, 24)
features, y = feature_engineering_outputs(TIME_HORIZONS, TimeSeriesSplitter())

We follow a multiple regressor approach and define separate models per quantile as follows:

• a model predicting the 5th percentile of the load

• a model predicting the median of the load

• a model predicting the 95th percentile of the load

We evaluate the performance of the quantile regressors via cross-validation.

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.ensemble import HistGradientBoostingRegressor


class HGBQuantileRegressor(RegressorMixin, BaseEstimator):
    def __init__(self, quantiles=(0.05, 0.5, 0.95), hgb_params=None):
        self.quantiles = quantiles
        self.hgb_params = hgb_params

    def fit(self, X, y):
        self.quantiles_ = sorted(self.quantiles)
        params = (self.hgb_params or {}) | {"loss": "quantile"}
        self.estimators_ = {
            q: HistGradientBoostingRegressor(quantile=q, **params).fit(X, y)
            for q in self.quantiles_
        }
        return self

    def predict(self, X):
        result = np.asarray([e.predict(X) for e in self.estimators_.values()])
        result.sort(axis=0)
        return pl.DataFrame(result, schema=[f"q_{q}" for q in self.quantiles_])


quantiles = (0.05, 0.5, 0.95)

hgb_params = get_hgb_params()

hgb_q_regressor = HGBQuantileRegressor(quantiles=quantiles, hgb_params=hgb_params)

def concat_horizons(all_pred, mode=skrub.eval_mode()):
    """
    Consolidate predictions of models for different horizons in one dataframe.
    """
    if mode == "fit":
        return all_pred
    return pl.concat(
        [v.rename(f"{h}h__{{}}".format) for h, v in all_pred.items()], how="horizontal"
    )


def make_multi_horizon_pred(features, y, regressor):
    """
    Create a full DataOp for predicting the specified horizons and quantiles.
    """
    predictions = {
        h: feat.skb.drop(["prediction_time", "target_time"])
        .skb.apply(regressor, y=y[f"{h}h"])
        .skb.set_name(f"pred_{h}h")
        for h, feat in features.items()
    }
    return skrub.deferred(concat_horizons)(predictions)


pred = (
    make_multi_horizon_pred(features, y, regressor=hgb_q_regressor)
    .skb.with_scoring(neg_mape_median_scorer)
    .skb.with_scoring(pinball_scorer)
)

pred

<Scoring <Call 'concat_horizons'> (2 scorers)>

    This DataOp will be scored with:
      - <function neg_mape_median_scorer at 0x701231769080>
      - <function pinball_scorer at 0x701231769300>
    Use .skb.cross_validate(…) or .skb.make_learner(…).score(…) to compute scores.

Show/Hide graph

Result (preview of the learner's output, not the scores):

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Let’s first collect all the cross-validated predictions to make further inspection.

cv_predictions_hgbr = get_cv_results(pred)

split 0: 2023-06-01T00:00:00+00:00
{'neg_mape_average': -0.022408965298349386, 'neg_mape_1h': -0.01583170542186339, 'neg_mape_12h': -0.02526267673039749, 'neg_mape_24h': -0.02613251374278728, 'd2_pinball_score__average__q_0.05': 0.6543313224198686, 'd2_pinball_score__1h__q_0.05': 0.779995946869003, 'd2_pinball_score__12h__q_0.05': 0.6413444656356533, 'd2_pinball_score__24h__q_0.05': 0.6276760027378196, 'd2_pinball_score__average__q_0.5': 0.797151195587183, 'd2_pinball_score__1h__q_0.5': 0.8593740568457844, 'd2_pinball_score__12h__q_0.5': 0.7697151514846812, 'd2_pinball_score__24h__q_0.5': 0.7623643784310833, 'd2_pinball_score__average__q_0.95': 0.6062372120615593, 'd2_pinball_score__1h__q_0.95': 0.3291586735848807, 'd2_pinball_score__12h__q_0.95': 0.2140171007660604, 'd2_pinball_score__24h__q_0.95': 0.23640719515650466}

split 1: 2023-09-01T00:00:00+00:00
{'neg_mape_average': -0.01843180078475345, 'neg_mape_1h': -0.014660201729158334, 'neg_mape_12h': -0.02000987590253614, 'neg_mape_24h': -0.020625324722565876, 'd2_pinball_score__average__q_0.05': 0.7907005196666222, 'd2_pinball_score__1h__q_0.05': 0.7368284550531783, 'd2_pinball_score__12h__q_0.05': 0.5893220923649396, 'd2_pinball_score__24h__q_0.05': 0.5881306256374493, 'd2_pinball_score__average__q_0.5': 0.8592470866963954, 'd2_pinball_score__1h__q_0.5': 0.8888361239908931, 'd2_pinball_score__12h__q_0.5': 0.8455533213276223, 'd2_pinball_score__24h__q_0.5': 0.8433518147706711, 'd2_pinball_score__average__q_0.95': 0.8185086691860827, 'd2_pinball_score__1h__q_0.95': 0.657819637266799, 'd2_pinball_score__12h__q_0.95': 0.5926998554803007, 'd2_pinball_score__24h__q_0.95': 0.6045767978526273}

split 2: 2023-12-01T00:00:00+00:00
{'neg_mape_average': -0.02147396070661371, 'neg_mape_1h': -0.01432850902117237, 'neg_mape_12h': -0.02506058437807546, 'neg_mape_24h': -0.025032788720593294, 'd2_pinball_score__average__q_0.05': 0.6244317779512799, 'd2_pinball_score__1h__q_0.05': 0.308744109714861, 'd2_pinball_score__12h__q_0.05': 0.23516667616139286, 'd2_pinball_score__24h__q_0.05': 0.28557771371536134, 'd2_pinball_score__average__q_0.5': 0.8015114582069321, 'd2_pinball_score__1h__q_0.5': 0.8694191673178604, 'd2_pinball_score__12h__q_0.5': 0.7683206214376956, 'd2_pinball_score__24h__q_0.5': 0.7667945858652403, 'd2_pinball_score__average__q_0.95': 0.7501683342903854, 'd2_pinball_score__1h__q_0.95': 0.721595647279887, 'd2_pinball_score__12h__q_0.95': 0.44936907338691456, 'd2_pinball_score__24h__q_0.95': 0.4625535008274011}

split 3: 2024-03-01T00:00:00+00:00
{'neg_mape_average': -0.023645029838124627, 'neg_mape_1h': -0.015857442444054386, 'neg_mape_12h': -0.027090510256754297, 'neg_mape_24h': -0.027987136813565204, 'd2_pinball_score__average__q_0.05': 0.7355070368383654, 'd2_pinball_score__1h__q_0.05': 0.7416140035373288, 'd2_pinball_score__12h__q_0.05': 0.5937248309067196, 'd2_pinball_score__24h__q_0.05': 0.5813806753969646, 'd2_pinball_score__average__q_0.5': 0.8102380014711067, 'd2_pinball_score__1h__q_0.5': 0.8703058324681592, 'd2_pinball_score__12h__q_0.5': 0.7845788335378543, 'd2_pinball_score__24h__q_0.5': 0.7758293384073065, 'd2_pinball_score__average__q_0.95': 0.7407556724771331, 'd2_pinball_score__1h__q_0.95': 0.6072058958523814, 'd2_pinball_score__12h__q_0.95': 0.3971662443515662, 'd2_pinball_score__24h__q_0.95': 0.39973672219579504}

split 4: 2024-06-01T00:00:00+00:00
{'neg_mape_average': -0.014531345512832734, 'neg_mape_1h': -0.009387012101896184, 'neg_mape_12h': -0.017197332812909222, 'neg_mape_24h': -0.017009691623692792, 'd2_pinball_score__average__q_0.05': 0.8241017390136228, 'd2_pinball_score__1h__q_0.05': 0.7820711572838047, 'd2_pinball_score__12h__q_0.05': 0.6517621015558843, 'd2_pinball_score__24h__q_0.05': 0.6337144560535894, 'd2_pinball_score__average__q_0.5': 0.8623396111518488, 'd2_pinball_score__1h__q_0.5': 0.9115092090327068, 'd2_pinball_score__12h__q_0.5': 0.8369587822423226, 'd2_pinball_score__24h__q_0.5': 0.838550842180517, 'd2_pinball_score__average__q_0.95': 0.6693224918571827, 'd2_pinball_score__1h__q_0.95': 0.4087330333070992, 'd2_pinball_score__12h__q_0.95': 0.3579286903540869, 'd2_pinball_score__24h__q_0.95': 0.3789800771107785}

split 5: 2024-09-01T00:00:00+00:00
{'neg_mape_average': -0.016574871291635106, 'neg_mape_1h': -0.01060111924957616, 'neg_mape_12h': -0.0196485731477636, 'neg_mape_24h': -0.019474921477565556, 'd2_pinball_score__average__q_0.05': 0.7922947774256607, 'd2_pinball_score__1h__q_0.05': 0.73764235207102, 'd2_pinball_score__12h__q_0.05': 0.5864008316458218, 'd2_pinball_score__24h__q_0.05': 0.6010896364082216, 'd2_pinball_score__average__q_0.5': 0.8666102553350941, 'd2_pinball_score__1h__q_0.5': 0.9153144677617081, 'd2_pinball_score__12h__q_0.5': 0.8415941093070951, 'd2_pinball_score__24h__q_0.5': 0.842922188936479, 'd2_pinball_score__average__q_0.95': 0.8230022444012137, 'd2_pinball_score__1h__q_0.95': 0.6867030989240508, 'd2_pinball_score__12h__q_0.95': 0.6215526542236005, 'd2_pinball_score__24h__q_0.95': 0.627078679176025}

split 6: 2024-12-01T00:00:00+00:00
{'neg_mape_average': -0.02002954456179126, 'neg_mape_1h': -0.011794786446836607, 'neg_mape_12h': -0.02462592559785109, 'neg_mape_24h': -0.02366792164068608, 'd2_pinball_score__average__q_0.05': 0.5422515688939357, 'd2_pinball_score__1h__q_0.05': -0.07019308402453239, 'd2_pinball_score__12h__q_0.05': 0.04139411870380938, 'd2_pinball_score__24h__q_0.05': 0.08689864395107971, 'd2_pinball_score__average__q_0.5': 0.7834765896845122, 'd2_pinball_score__1h__q_0.5': 0.8721777103959318, 'd2_pinball_score__12h__q_0.5': 0.7346604742897471, 'd2_pinball_score__24h__q_0.5': 0.7435915843678574, 'd2_pinball_score__average__q_0.95': 0.7421152370795165, 'd2_pinball_score__1h__q_0.95': 0.702277126611818, 'd2_pinball_score__12h__q_0.95': 0.44304532678156394, 'd2_pinball_score__24h__q_0.95': 0.47471138977020944}

split 7: 2025-03-01T00:00:00+00:00
{'neg_mape_average': -0.020647583404155474, 'neg_mape_1h': -0.012746538051496864, 'neg_mape_12h': -0.024130027202063517, 'neg_mape_24h': -0.025066184958906043, 'd2_pinball_score__average__q_0.05': 0.7870409791714095, 'd2_pinball_score__1h__q_0.05': 0.7675272792873958, 'd2_pinball_score__12h__q_0.05': 0.6459032937670941, 'd2_pinball_score__24h__q_0.05': 0.6421235435312523, 'd2_pinball_score__average__q_0.5': 0.8498500354300744, 'd2_pinball_score__1h__q_0.5': 0.9070055963308842, 'd2_pinball_score__12h__q_0.5': 0.8251259668325958, 'd2_pinball_score__24h__q_0.5': 0.8174185431267433, 'd2_pinball_score__average__q_0.95': 0.8087457151876744, 'd2_pinball_score__1h__q_0.95': 0.6311001403225274, 'd2_pinball_score__12h__q_0.95': 0.47483796748808904, 'd2_pinball_score__24h__q_0.95': 0.4770705022666911}

Now, we can inspect the cross-validated predictions and plot them for the different quantiles.

plot_predictions(cv_predictions_hgbr[0], horizon=12, start="2023-01-01").show()

Now, let’s collect the cross-validated predictions and plot the residual vs predicted values for the different models into a report.

cv_predictions_hgbr[0].head(5)

shape: (5, 14)

prediction_time	1h	12h	24h	pred_1h__q_0.05	pred_1h__q_0.5	pred_1h__q_0.95	pred_12h__q_0.05	pred_12h__q_0.5	pred_12h__q_0.95	pred_24h__q_0.05	pred_24h__q_0.5	pred_24h__q_0.95	split
datetime[μs, UTC]	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	i32
2023-06-01 00:00:00 UTC	36635.0	48413.0	37880.0	35530.837285	35809.494355	40971.274655	44549.282685	45328.905999	47568.327205	36605.478511	37636.537264	40650.877915	0
2023-06-01 01:00:00 UTC	34759.0	47436.0	36351.0	34418.467282	35171.849994	40971.274655	43228.708684	44225.202966	47174.362691	35450.140842	35735.475731	40465.624785	0
2023-06-01 02:00:00 UTC	34081.0	46226.0	34385.0	34352.651495	35556.415574	40971.274655	42480.998412	43170.894509	46053.099837	33035.506552	34065.711341	40465.624785	0
2023-06-01 03:00:00 UTC	35248.0	43971.0	33693.0	34816.00478	36086.052601	41463.293895	42738.755651	44139.337298	46396.404838	33398.491357	34567.096174	40465.624785	0
2023-06-01 04:00:00 UTC	38346.0	44276.0	34955.0	37063.841556	38434.2166	42030.748727	43312.276785	44262.820291	46704.71722	34589.629576	36387.745458	41235.856864	0

Focusing on the different D2 scores, we observe that each model minimize the D2 score associated to the target quantile that we set. For instance, the model predicting the 5th percentile obtained the highest D2 pinball score with alpha=0.05. It is expected but a confirmation of what loss each model minimizes.

Now, let’s collect the cross-validated predictions and plot the residual vs predicted values for the different models.

plot_residuals_vs_predicted(cv_predictions_hgbr[0], 1, quantile=0.05).interactive().properties(
    title=("Residuals vs Predicted Values from cross-validation predictions" " for quantile 0.05")
)

plot_residuals_vs_predicted(cv_predictions_hgbr[0], 1, quantile=0.5).interactive().properties(
    title=("Residuals vs Predicted Values from cross-validation predictions for median")
)

plot_residuals_vs_predicted(cv_predictions_hgbr[0], 1, quantile=0.95).interactive().properties(
    title=("Residuals vs Predicted Values from cross-validation predictions" " for quantile 0.95")
)

We observe an expected behaviour: the residuals are centered and symmetric around 0 for the median model while not centered and biased for the 5th and 95th percentiles models.

Now, we assess if the actual coverage of the models is close to the target coverage of 90%. In addition, we compute the average width of the bands.

import altair

from tutorial_helpers import binned_coverage, coverage

preds = cv_predictions_hgbr[0]
horizon = 1

# --- Overall coverage per fold ---
for (split_idx,), fold_df in preds.group_by("split", maintain_order=True):
    cov = coverage(
        fold_df[f"{horizon}h"].to_numpy(),
        fold_df[f"pred_{horizon}h__q_0.05"].to_numpy(),
        fold_df[f"pred_{horizon}h__q_0.95"].to_numpy(),
    )
    print(f"Split {split_idx}: {cov:.1%} coverage (90% interval)")

# --- Binned coverage plot ---
folds = [fold_df for (_,), fold_df in preds.group_by("split", maintain_order=True)]
binned = binned_coverage(
    y_true_folds=[f[f"{horizon}h"].to_numpy() for f in folds],
    y_quantile_low=[f[f"pred_{horizon}h__q_0.05"].to_numpy() for f in folds],
    y_quantile_high=[f[f"pred_{horizon}h__q_0.95"].to_numpy() for f in folds],
)

altair.Chart(binned).mark_line(point=True).encode(
    x=altair.X("bin_center:Q", title="True load (MW)"),
    y=altair.Y("coverage:Q", title="Coverage", scale=altair.Scale(domain=[0, 1])),
    color=altair.Color("fold_idx:N"),
).properties(title=f"Binned coverage — {horizon}h horizon, 90% interval")

Split 0: 76.5% coverage (90% interval)
Split 1: 81.0% coverage (90% interval)
Split 2: 87.6% coverage (90% interval)
Split 3: 85.3% coverage (90% interval)
Split 4: 90.3% coverage (90% interval)
Split 5: 86.6% coverage (90% interval)
Split 6: 90.5% coverage (90% interval)
Split 7: 90.0% coverage (90% interval)

We observe that the lower and higher bins, so low and high load, have the worse coverage with a high variability.

Reliability diagrams and Lorenz curves for quantile regression

plot_reliability_diagram(
    cv_predictions_hgbr[0], 1, forecast_quantile=0.50
).interactive().properties(
    title="Reliability diagram for quantile 0.50 from cross-validation predictions"
)

plot_reliability_diagram(
    cv_predictions_hgbr[0], 1, forecast_quantile=0.05
).interactive().properties(
    title="Reliability diagram for quantile 0.05 from cross-validation predictions"
)

plot_reliability_diagram(
    cv_predictions_hgbr[0], 1, forecast_quantile=0.95
).interactive().properties(
    title="Reliability diagram for quantile 0.95 from cross-validation predictions"
)

plot_lorenz_curve(cv_predictions_hgbr[0], 1, quantile=0.50).interactive().properties(
    title="Lorenz curve for quantile 0.50 from cross-validation predictions"
)

plot_lorenz_curve(cv_predictions_hgbr[0], 1, quantile=0.05).interactive().properties(
    title="Lorenz curve for quantile 0.05 from cross-validation predictions"
)

plot_lorenz_curve(cv_predictions_hgbr[0], 1, quantile=0.95).interactive().properties(
    title="Lorenz curve for quantile 0.95 from cross-validation predictions"
)