How are samples drawn from conformal prediction model / quantiles?

[JuanCruzC97]

Hi, I'm using regression models together with a ConformalNaive model to generate probabilistic forecasts. I'm experimenting with two approaches:

1. Generating historical quantile forecasts using predict_likelihood_parameters=True.
2. Generating historical forecasts drawing a sampled distribution with num_samples=1000.

However, I'm confused about how the samples are being generated. For the same prediction date, I get the following quantiles:

q0.05 = 359.11175744  
q0.25 = 380.32542817  
q0.50 = 407.34136384  
q0.75 = 434.35729951  
q0.95 = 455.57097025

But when I inspect the sampled distribution, I get the following stats:

count    1000.000000  
mean      409.045300  
std        30.632939  
min       359.111757  
25%       382.522803  
50%       409.312385  
75%       436.266100  
max       455.570970

The 0.05 and 0.95 quantiles match the minimum and maximum of the sampled values, and their frequencies are unusually high:

455.570970    59  
359.111757    43  
(other values) ~1 each

This leads to a histogram with spikes at the ends of the distribution.

My questions:

• How are the samples generated from the quantiles?
• Is it expected that the extreme quantiles are overrepresented like this in the sample?
• Shouldn’t the samples more closely reflect a smooth distribution? I'm expecting the 0.05 quantile and the 0.95 quantile to cover the 90% of the sampled distribution but this will get me always 100% coverage.
• Am I doing something wrong here?

Reproducible Example

import pandas as pd

from darts import concatenate, metrics, TimeSeries
from darts.datasets import AirPassengersDataset
from darts.models import ConformalNaiveModel, LinearRegressionModel

series = AirPassengersDataset().load()

train_start = pd.Timestamp("1949-01-01")
cal_start = pd.Timestamp("1957-01-01")
test_start = pd.Timestamp("1959-01-01")
test_end = pd.Timestamp("1960-12-01")

train = series[train_start : cal_start - series.freq]
cal = series[cal_start : test_start - series.freq]
test = series[test_start:test_end]
cal_test = concatenate([cal, test])

multi_horizon = 3
quantiles = [0.05, 0.25, 0.50, 0.75, 0.95]
input_length = 10

model = LinearRegressionModel(
    lags=input_length, 
    output_chunk_length=multi_horizon, 
    use_static_covariates=False
)

model.fit(train)

cp_model = ConformalNaiveModel(
    model=model, 
    quantiles=quantiles
)

comformal_samples = cp_model.historical_forecasts(
    series=cal_test,
    start=test_start,
    forecast_horizon=multi_horizon,
    retrain=False,
    num_samples=1000,
    predict_likelihood_parameters=False,
    last_points_only=True,
)

comformal_quantiles = cp_model.historical_forecasts(
    series=cal_test,
    start=test_start,
    forecast_horizon=multi_horizon,
    retrain=False,
    num_samples=1,
    predict_likelihood_parameters=True,
    last_points_only=True,    
)

Thanks!

[dennisbader]

Hi @JuanCruzC97, let me try to answer your questions.

Conformal prediction

First of all, let's look at the requirements / limitations of the conformal methods (see the disadvantages here):

Requires a Calibration Set: Conformal prediction requires another data / hold-out set that is used solely to compute the calibrated prediction intervals. This can be inefficient for small datasets.
Exchangeability of Calibration Data (a): The accuracy of the prediction intervals depends on the representativeness of the calibration data (or rather the forecast errors produced on the calibration set). The coverage is not guaranteed anymore if there is a distribution shift in forecast errors (e.g. series with a trend but forecasting model is not able to predict the trend).

The raw AirPassengers series is a tough example for conformal prediction because the size is rather small (e.g. your training, calibration and test sets), plus it has a trend and increasing seasonality over time. This means that the coverage cannot be guaranteed anymore if your model errors increase over time.

In your example you also only look at the coverage of the last predicted point of each historical forecast (last_points_only=True). In total, I believe the forecasts on your test set ended up with around 12 of these last points. Having a coverage of 100% over 12 points isn't enough to say that the coverage will not converge to 90% (e.g. if the next predicted point was covered then we would are already be much closer to 90%).

So to summarize this part:

• conformal prediction only guarantees coverage if the errors of your model in the calibration set are comparable with those in the test set
• your dataset size is rather small, which also has a deteriorating effect on the coverage guarantees

Sampling

How are the samples generated from the quantiles?

You can find the description and logic here:

darts/darts/utils/utils.py

Line 606 in b71851b

def sample_from_quantiles(

In short, the logic is:

• Generate num_samples uniform random samples between [0, 1]
• Find the lower and upper quantile between which each sample falls into
• Linearly interpolate between the values of these quantiles with the sample's value
• For samples that fall below the lowest or above the highest quantile, repeat these extreme quantile values (since we cannot interpolate as we don't know about the tails).

(...) But when I inspect the sampled distribution, I get the following stats:

If you increase the number of samples, the computed quantile values on the samples will converge towards the quantile values from predict_likelihood_parameters=True.

Is it expected that the extreme quantiles are overrepresented like this in the sample?

Yes, it's expected. Since we don't know about the tails of the distribution (e.g. below q=0.05 and above q=0.95), we repeat the values, so that if you compute the extreme quantiles over the samples, you should get approximately the same quantile values.

Shouldn’t the samples more closely reflect a smooth distribution? I'm expecting the 0.05 quantile and the 0.95 quantile to cover the 90% of the sampled distribution but this will get me always 100% coverage.

Sampling from quantiles is only an approximation of the actual distribution. The more quantiles you add, the more the samples will reflect a smooth distribution. See for example below where I fitted a conformal model on the ElectricityConsumptionZurichDataset with all quantiles [0.01, 0.02, ... 0.99].

I hope this helped :)

[JuanCruzC97]

Hi again, and thanks for the thorough explanation!

I’ll definitely dive deeper into the theory of conformal prediction—your pointers are much appreciated. I initially chose the AirPassengers dataset because it highlights the “heavy spots” in the error distribution, which mirrors the behavior I’m seeing in my own calibration set (large sample, no clear trend, and no increasing variance). I realize it isn’t the perfect example to append here.

Your walkthrough of how the samples are generated—and why those heavy spots arise—really clarified things. That brings me to the one remaining wrinkle I’m encountering:

• Quantile crossing at median (0.50) when symmetric=False.
  In my application, the predicted 0.50 quantile consistently crosses over neighboring quantiles. Is this behavior anticipated under symmetric=False, or does it suggest something wrong in my setup?

I’ve attached two plots illustrating the crossing issues for reference:

Thanks again for all your time and help!

[dennisbader]

Hi @JuanCruzC97 and sorry for the late reply. It is expected in the way that we always return the actual forecasting model forecasts as the 0.5 quantile (not the actual conformalized predictions at quantile 0.5).

For the symmetric case this makes sense, as the upper and lower quantiles have the same absolute distance to the model forecast. In that case the conformalized quantile intervals also give the "uncertainty" of the model. The larger / wider the intervals, the higher the uncertainty.

In the asymmetric case, the upper and lower quantiles are computed independently on the residual distribution. If your model has a bias, it can happen that the quantile prediction can shoot under or above the model forecast (the chance gets larger the closer the conformalized quantiles to 0.5, as you see in your image). We chose to still return the forecasting model prediction at q0.5 to keep a connection to the underlying model.

Maybe we could think of changing this behavior for q0.5 for the asymmetric case though and actually return the conformalized quantile 0.5..

[swertz]

Maybe we could think of changing this behavior for q0.5 for the asymmetric case though and actually return the conformalized quantile 0.5..

This touches on something we were wondering about when using a conformalized quantile regressor: the 50% quantile seemed to not be affected by the calibration step. So you're saying that the 50% one is always taken from the raw quantile regression output without calibration?

If that's case I agree that it's important to change it to handle all quantiles consistently!