Computing prediction intervals using quantile regression

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 re
import warnings
from datetime import datetime

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 single_horizon_prediction_lib import TimeSeriesSplitter
from tutorial_helpers import (
    binned_coverage,
    collect_cv_predictions,
    plot_lorenz_curve,
    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(y_true, y_pred, quantile_regression=False):
    if quantile_regression:
        quantile_predictions = split_by_quantile(y_pred)
        scores = {}
        for q, q_pred in quantile_predictions.items():
            q_neg_mape = neg_mape(y_true, q_pred, quantile_regression=False)
            scores.update({f"{k}__{q}": v for k, v in q_neg_mape.items()})
            if q == "q_0.5":
                # Pick the median if available for comparison with non-quantile
                # models
                scores.update(q_neg_mape)
        return scores
    average = mean_absolute_percentage_error(y_true, y_pred)
    detail = mean_absolute_percentage_error(y_true, y_pred, multioutput="raw_values")
    return {"neg_mape__average": -average} | {
        f"neg_mape__{c}": -float(s) for c, s in zip(y_true.columns, detail)
    }


def neg_mape_scorer(estimator, X, y, quantile_regression=False):
    return neg_mape(y, estimator.predict(X), quantile_regression=quantile_regression)


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)

learning_rate = skrub.choose_float(0.01, 0.7, default=0.1, log=True, name="learning_rate")
max_leaf_nodes = skrub.choose_int(3, 300, default=30, log=True, name="max_leaf_nodes")
hgb_params = dict(
    random_state=0,
    learning_rate=learning_rate,
    max_leaf_nodes=max_leaf_nodes,
)

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.
    """
    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(functools.partial(neg_mape_scorer, quantile_regression=True))
    .skb.with_scoring(pinball_scorer)
)

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

def concat_X_y_predictions(X_test, y_test, prediction):
    return pl.concat(
        [
            X_test,
            y_test,
            prediction.rename("pred_{}".format),
        ],
        how="horizontal",
    )


def cross_val_predict(data_op, environment=None):
    """
    Get cross-validated predictions for different horizons.
    """
    all_predictions, all_scores = [], []
    for i, split in enumerate(data_op.skb.iter_cv_splits(environment=environment)):
        learner = data_op.skb.make_learner().fit(split["train"])
        score, predictions = learner.score(split["test"], return_predictions=True)
        all_predictions.append(
            concat_X_y_predictions(
                split["X_test"], split["y_test"], predictions["predict"]
            ).with_columns(split=pl.lit(i)),
        )
        print(split["X_test"]["prediction_time"].min().isoformat())
        print(score)
        all_scores.append(score | {"split": i})
    all_predictions = pl.concat(all_predictions, how="vertical")
    all_scores = pl.DataFrame(all_scores)
    return all_predictions, all_scores


cv_predictions_hgbr = cross_val_predict(
    pred, environment={"start": "2023-01-01", "end": "2025-05-31"}
)

2025-02-01T00:00:00+00:00
{'neg_mape__average__q_0.05': -0.0547254425058962, 'neg_mape__1h__q_0.05': -0.0496482227650003, 'neg_mape__12h__q_0.05': -0.05830689772690905, 'neg_mape__24h__q_0.05': -0.05622120702577925, 'neg_mape__average__q_0.5': -0.01937438657698719, 'neg_mape__1h__q_0.5': -0.01182505784162459, 'neg_mape__12h__q_0.5': -0.022765864066416557, 'neg_mape__24h__q_0.5': -0.023532237822920424, 'neg_mape__average': -0.01937438657698719, 'neg_mape__1h': -0.01182505784162459, 'neg_mape__12h': -0.022765864066416557, 'neg_mape__24h': -0.023532237822920424, 'neg_mape__average__q_0.95': -0.055956368677507064, 'neg_mape__1h__q_0.95': -0.039386390651328526, 'neg_mape__12h__q_0.95': -0.06436989198542618, 'neg_mape__24h__q_0.95': -0.0641128233957665, 'd2_pinball_score__average__q_0.05': 0.7614808833563985, 'd2_pinball_score__1h__q_0.05': 0.6082885133111983, 'd2_pinball_score__12h__q_0.05': 0.5660010980441968, 'd2_pinball_score__24h__q_0.05': 0.5884211951020477, 'd2_pinball_score__average__q_0.5': 0.8713531164805067, 'd2_pinball_score__1h__q_0.5': 0.9187021455168717, 'd2_pinball_score__12h__q_0.5': 0.8494418664369818, 'd2_pinball_score__24h__q_0.5': 0.8459153374876669, 'd2_pinball_score__average__q_0.95': 0.8127335915718662, 'd2_pinball_score__1h__q_0.95': 0.742578004547708, 'd2_pinball_score__12h__q_0.95': 0.5868203808992172, 'd2_pinball_score__24h__q_0.95': 0.5931879817743396}

2025-05-01T00:00:00+00:00
{'neg_mape__average__q_0.05': -0.03635870000916958, 'neg_mape__1h__q_0.05': -0.024633488977422787, 'neg_mape__12h__q_0.05': -0.042201122511490365, 'neg_mape__24h__q_0.05': -0.04224148853859561, 'neg_mape__average__q_0.5': -0.02043916263735127, 'neg_mape__1h__q_0.5': -0.01383459120008417, 'neg_mape__12h__q_0.5': -0.023988061465548327, 'neg_mape__24h__q_0.5': -0.02349483524642131, 'neg_mape__average': -0.02043916263735127, 'neg_mape__1h': -0.01383459120008417, 'neg_mape__12h': -0.023988061465548327, 'neg_mape__24h': -0.02349483524642131, 'neg_mape__average__q_0.95': -0.0777483429416056, 'neg_mape__1h__q_0.95': -0.07499432756700079, 'neg_mape__12h__q_0.95': -0.08150426735762355, 'neg_mape__24h__q_0.95': -0.07674643390019245, 'd2_pinball_score__average__q_0.05': 0.7053171352283405, 'd2_pinball_score__1h__q_0.05': 0.7421214488915335, 'd2_pinball_score__12h__q_0.05': 0.5546607349584718, 'd2_pinball_score__24h__q_0.05': 0.5490494326767557, 'd2_pinball_score__average__q_0.5': 0.7963716818065506, 'd2_pinball_score__1h__q_0.5': 0.8650488535736427, 'd2_pinball_score__12h__q_0.5': 0.7610019276242612, 'd2_pinball_score__24h__q_0.5': 0.7630642642217479, 'd2_pinball_score__average__q_0.95': 0.6117939290717587, 'd2_pinball_score__1h__q_0.95': 0.30237361184783385, 'd2_pinball_score__12h__q_0.95': 0.21438511630211277, 'd2_pinball_score__24h__q_0.95': 0.2516684383103205}

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

from datetime import UTC, timedelta

import plotly.graph_objects as go


def plot_predictions(results, horizons=None, start="2025-03-01"):
    if start is not None:
        results = results.filter(
            pl.col("prediction_time") > datetime.fromisoformat(start).replace(tzinfo=UTC)
        )
    if horizons is None:
        horizons = sorted(
            {
                int(m.group(1))
                for c in results.columns
                if (m := re.match(r"^pred_(\d+)h.*$", c)) is not None
            }
        )
    fig = go.Figure()
    for i, h in enumerate(horizons):
        target_time = (results["prediction_time"] + timedelta(hours=h)).to_list()
        if not i:
            fig.add_trace(
                go.Scatter(
                    x=target_time,
                    y=results[f"{h}h"].to_list(),
                    mode="lines+markers",
                    line={"dash": "dash"},
                    name="true_load_mw",
                    hovertemplate="%{x|%Y-%m-%d} (%{x|%A}): %{y}<extra></extra>",
                )
            )
        for col in filter(lambda c: f"pred_{h}h" in c, results.columns):
            fig.add_trace(
                go.Scatter(
                    x=target_time,
                    y=results[col].to_list(),
                    mode="lines+markers",
                    name=col,
                    hovertemplate="%{x|%Y-%m-%d} (%{x|%A}): %{y}<extra></extra>",
                )
            )
    fig.update_layout(height=600, title="CV predicted load mw")
    return fig


plot_predictions(cv_predictions_hgbr[0], horizons=(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_time1h12h24hpred_1h__q_0.05pred_1h__q_0.5pred_1h__q_0.95pred_12h__q_0.05pred_12h__q_0.5pred_12h__q_0.95pred_24h__q_0.05pred_24h__q_0.5pred_24h__q_0.95split
datetime[μs, UTC]f64f64f64f64f64f64f64f64f64f64f64f64i32
2025-02-01 00:00:00 UTC58241.064739.7561082.2556857.58974959216.58679460217.74010960563.6943665181.62269667076.15511856837.59077661667.31837462975.9057090
2025-02-01 01:00:00 UTC56251.561283.2559910.2555428.49169355925.88347358192.69023759830.03110162855.50549565522.30669355967.57631861087.42873861206.908440
2025-02-01 02:00:00 UTC55244.7558958.7557888.554131.32274654452.66401457875.86626658929.75873560880.08618462754.41479254910.60709658820.60106860000.5598670
2025-02-01 03:00:00 UTC56038.7558091.556874.053986.16211556533.19792958706.23017356997.42718260108.99735162301.99783754061.04441458108.63789759166.4797620
2025-02-01 04:00:00 UTC57954.059819.7557455.054782.41701558510.96015561420.31721859600.32585362535.19106665447.03385654559.5363258542.95984960163.2502460

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: 90.2% coverage (90% interval)
Split 1: 87.1% 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"
)