Conformal Prediction & Uncertainty Quantification#

Source Files

twiga/distributions/conformal/core.py
twiga/distributions/conformal/base.py
twiga/distributions/conformal/cqr.py
twiga/distributions/conformal/residual_conformal
twiga/core/config/base.py

Conformal prediction provides distribution-free prediction intervals with finite-sample coverage guarantees. Unlike Bayesian methods, conformal prediction makes no assumptions about the data distribution - it only requires exchangeability of calibration data.

Because calibration needs only a set of predictions and ground truth from a held-out split, conformal prediction wraps any trained Twiga model without retraining - including plain point-forecast ML models like CatBoost or XGBoost. This makes it the recommended path for adding uncertainty estimates to an existing point forecaster.

How It Works#

        graph LR
    A[Train Model] --> B[Calibrate on Held-Out Data]
    B --> C[Compute Non-Conformity Scores]
    C --> D["Calculate Threshold q̂ at (1-α) quantile"]
    D --> E[Generate Intervals on New Data]
    E --> F["[prediction - q̂, prediction + q̂]"]

Train the model on training data
Calibrate using held-out data to compute non-conformity scores
Calculate threshold \(\hat{q}\) as the \((1-\alpha)\)-quantile of the scores
Generate intervals on new data using the calibrated threshold

The resulting intervals have guaranteed marginal coverage of at least \(1-\alpha\).

Class Hierarchy#

        classDiagram
    class BaseConformal {
        <<abstract>>
        +alpha: float
        +calib_method: str
        +lambda_: float
        +q_hat: float | ndarray
        +calibrate(*calib_args, axis)
        +calculate_conformal_quantile(scores, axis, calib_timestamps)
        +get_scores()*
        +generate_intervals()*
    }

    class SplitConformal {
        +score_type: "res" | "sign-res"
        +get_scores(predicts, targets)
        +generate_intervals(mu_pred)
    }

    class ConformalQuantileRegressor {
        +score_type: "scaled" | "unscaled"
        +get_scores(lower_q, upper_q, targets)
        +generate_intervals(lower_q, upper_q)
    }

    class ConformalResidualFitting {
        +score_type: "res" | "sign-res"
        +q_hat_low: float | ndarray
        +get_scores(predicts, sigma, targets)
        +calibrate(predicts, sigma, targets, calib_timestamps)
        +generate_intervals(loc, sigma)
    }

    class Conformal {
        <<factory>>
        +__new__(method, score_type, alpha, calib_method, lambda_)
    }

    BaseConformal <|-- SplitConformal
    BaseConformal <|-- ConformalQuantileRegressor
    BaseConformal <|-- ConformalResidualFitting
    Conformal ..> SplitConformal : creates
    Conformal ..> ConformalQuantileRegressor : creates
    Conformal ..> ConformalResidualFitting : creates

Factory: `Conformal`#

The Conformal class in twiga/distributions/conformal/core.py is a factory that selects the appropriate method:

from twiga.distributions.conformal.core import Conformal

# Creates a SplitConformal instance
cp = Conformal(method="residual", score_type="res", alpha=0.1)

# Creates a ConformalQuantileRegressor instance
cqr = Conformal(method="quantile", score_type="scaled", alpha=0.1)

# Creates a ConformalResidualFitting instance with uniform calibration
crf = Conformal(method="residual-fitting", score_type="sign-res", alpha=0.1)

# Creates a ConformalResidualFitting instance with temporal weighting
crf_t = Conformal(
    method="residual-fitting",
    score_type="sign-res",
    alpha=0.1,
    calib_method="temporal",
    lambda_=1.0,
)

calib_method and lambda_ are forwarded only when method="residual-fitting".

Configuration#

Conformal prediction is configured via ConformalConfig:

from twiga.core.config import ConformalConfig

# Standard uniform calibration:
config = ConformalConfig(
    method="residual-fitting",
    score_type="sign-res",
    alpha=0.1,
)

# Temporal calibration weighting (addresses seasonal distribution shift):
config = ConformalConfig(
    method="residual-fitting",
    score_type="sign-res",
    alpha=0.1,
    calib_method="temporal",
    lambda_=1.0,
)

# Via coverage probability:
config = ConformalConfig.from_coverage(
    coverage=0.9,
    method="residual",
    score_type="res",
)

Parameter	Type	Constraints	Description
`method`	`Literal["residual", "quantile", "residual-fitting"]`	Required	Conformal prediction method
`score_type`	`Literal["scaled", "unscaled", "res", "sign-res"]`	Required	Non-conformity score type
`alpha`	`float`	`0 < alpha < 1`	Significance level (1 minus alpha = target coverage)
`calib_method`	`Literal["uniform", "temporal"]`	Default `"uniform"`	Quantile estimation method for calibration scores
`lambda_`	`float`	`>= 0`, default `1.0`	Exponential decay rate for temporal weighting

ConformalConfig.from_coverage(coverage, **kwargs) is a convenience factory that converts a coverage probability to alpha = 1 - coverage and forwards any remaining keyword arguments to the constructor.

Valid Score Types per Method#

Method	Class	Valid Score Types	Description
`"residual"`	`SplitConformal`	`"res"`, `"sign-res"`	Absolute or signed residuals
`"quantile"`	`ConformalQuantileRegressor`	`"scaled"`, `"unscaled"`	CQR with or without scaling
`"residual-fitting"`	`ConformalResidualFitting`	`"res"`, `"sign-res"`	Scale-adapted residuals

Methods#

Split Conformal (`method="residual"`)#

The simplest method. Computes non-conformity scores as residuals between predictions and targets.

Score types:

"res": \(s_i = |y_i - \hat{y}_i|\) (absolute residuals)
"sign-res": \(s_i = y_i - \hat{y}_i\) (signed residuals)

Intervals: \([\hat{y} - \hat{q}, \; \hat{y} + \hat{q}]\)

from twiga.distributions.conformal.base import SplitConformal

cp = SplitConformal(score_type="res", alpha=0.1)
cp.calibrate(predictions, targets)
lower, upper = cp.generate_intervals(new_predictions)

Tip

Split conformal works with any point prediction model and is the easiest to set up. Use this as a starting point.

Conformal Quantile Regression (`method="quantile"`)#

Requires a quantile regression model that produces lower and upper quantile predictions.

Score types:

"unscaled": \(s_i = \max(q_{lo,i} - y_i, \; y_i - q_{hi,i})\)
"scaled": \(s_i = \max\left(\frac{q_{lo,i} - y_i}{q_{hi,i} - q_{lo,i}}, \; \frac{y_i - q_{hi,i}}{q_{hi,i} - q_{lo,i}}\right)\)

Intervals (unscaled): \([q_{lo} - \hat{q}, \; q_{hi} + \hat{q}]\)

Intervals (scaled): \([q_{lo} - \hat{q} \cdot (q_{hi} - q_{lo}), \; q_{hi} + \hat{q} \cdot (q_{hi} - q_{lo})]\)

from twiga.distributions.conformal.cqr import ConformalQuantileRegressor

cqr = ConformalQuantileRegressor(score_type="scaled", alpha=0.1)
cqr.calibrate(lower_quantile, upper_quantile, targets)
lower, upper = cqr.generate_intervals(lower_quantile_new, upper_quantile_new)

Note

Scaled CQR adapts interval width based on the model’s quantile spread, producing narrower intervals where the model is more confident. Use with quantile regression models.

Conformal Residual Fitting (`method="residual-fitting"`)#

Requires a model that produces both a point prediction (\(\mu\)) and a scale estimate (\(\sigma\)).

Score types:

"res": \(s_i = \frac{|\hat{y}_i - y_i|}{\sigma_i}\) (absolute)
"sign-res": \(s_i = \frac{\hat{y}_i - y_i}{\sigma_i}\) (signed, can be negative)

Intervals (res): \([\mu - \sigma \cdot \hat{q}, \; \mu + \sigma \cdot \hat{q}]\)

Intervals (sign-res): \([\mu - \sigma \cdot \hat{q}^+, \; \mu - \sigma \cdot \hat{q}^-]\) where \(\hat{q}^+ = Q_{1-\alpha/2}\) and \(\hat{q}^- = Q_{\alpha/2}\) of the signed scores, giving asymmetric bounds with \(1-\alpha\) marginal coverage.

from twiga.distributions.conformal.crc import ConformalResidualFitting

# Uniform calibration (default):
crf = ConformalResidualFitting(score_type="sign-res", alpha=0.1)
crf.calibrate(predictions, sigma, targets)
lower, upper = crf.generate_intervals(new_predictions, new_sigma)

# Temporal calibration (pass calibration timestamps):
import pandas as pd

crf_t = ConformalResidualFitting(
    score_type="sign-res",
    alpha=0.1,
    calib_method="temporal",
    lambda_=1.0,
)
calib_timestamps = pd.to_datetime(calib_df["timestamp"])
crf_t.calibrate(predictions, sigma, targets, calib_timestamps=calib_timestamps)
lower, upper = crf_t.generate_intervals(new_predictions, new_sigma)

Note

Use residual-fitting with probabilistic models like GAUSSCATBOOSTModel or neural models with sigma outputs like MLPGAFModel.

Temporal Calibration Weighting#

Standard conformal prediction assumes exchangeability: calibration and test samples come from the same distribution. In rolling-origin forecasting this assumption can be violated when calibration data is drawn from a different season than the test window.

Temporal weighting (Tibshirani et al. 2019) addresses this by assigning higher weight to calibration samples that are more recent and therefore closer in distribution to the test period.

Weight formula:

\[w_i \propto \exp(-\lambda \cdot d_i), \quad d_i = \frac{t_{\max} - t_i}{t_{\max} - t_{\min}}\]

where \(d_i \in [0, 1]\) is the normalised age of sample \(i\) (\(d_i = 0\) for the newest, \(d_i = 1\) for the oldest). Weights are normalised to sum to 1.

The weighted quantile selects the smallest score \(s_{(k)}\) such that:

\[\sum_{j=1}^{k} w_{(j)} \geq 1 - \alpha\]

Decay rate lambda_:

Value	Effect
`0.0`	Uniform weighting (same as `calib_method="uniform"`)
`1.0`	Moderate recency bias (recommended starting point)
`2.0`	Strong recency bias; upweights the most recent ~30 % of calibration data

When to use temporal weighting

Use calib_method="temporal" when the calibration window is drawn from a different season than the test window, for example in expanding-window cross-validation where calibration is always at the training tail while the test window advances into a new season. Temporal weighting does not change the CV protocol and is compatible with any calib_size.

from twiga.core.config import ConformalConfig
from twiga.forecaster.core import TwigaForecaster

conformal_config = ConformalConfig(
    method="residual-fitting",
    score_type="sign-res",
    alpha=0.1,
    calib_method="temporal",
    lambda_=1.0,
)

forecaster = TwigaForecaster(
    data_params=data_config,
    model_params=[model_cfg],
    cv_params=exp_config,
    conformal_params=conformal_config,
)

forecaster.fit(train_df=train_df)
forecaster.calibrate(calibrate_df=calib_df)  # timestamps extracted automatically
lower, point, upper = forecaster.predict_interval(test_df=test_df)

TwigaForecaster.calibrate() extracts per-sample timestamps from the calibration data automatically and passes them to ConformalResidualFitting.calibrate() when calib_method="temporal".

Helper Functions#

The weighting utilities are available as standalone functions for custom workflows:

import pandas as pd
from twiga.distributions.conformal.base import (
    calculate_temporal_weights,
    calculate_weighted_conformal_value,
)

timestamps = pd.date_range("2021-09-01", periods=90, freq="D")
weights = calculate_temporal_weights(timestamps, lambda_=1.0)
# weights: ndarray of shape (90,) summing to 1, newest sample gets highest weight

threshold = calculate_weighted_conformal_value(scores, alpha=0.1, weights=weights)
# threshold: scalar (1-D scores) or ndarray of shape (H,) for horizon-wise scores

Integration with TwigaForecaster#

The TwigaForecaster manages conformal prediction end-to-end:

from twiga.core.config import (
    ConformalConfig,
    DataPipelineConfig,
    ExperimentConfig,
)
from twiga.forecaster.core import TwigaForecaster
from twiga.models.ml.xgboost_model import XGBOOSTConfig

# 1. Configure with conformal prediction
conformal_config = ConformalConfig(
    method="residual",
    score_type="res",
    alpha=0.1,
)

forecaster = TwigaForecaster(
    data_params=data_config,
    model_params=[XGBOOSTConfig()],
    cv_params=train_config,
    conformal_params=conformal_config,
)

# 2. Train
forecaster.fit(train_df=train_df, val_df=val_df)

# 3. Calibrate on held-out data
forecaster.calibrate(calibrate_df=calibration_df)

# 4. Generate prediction intervals
interval_dict, times = forecaster.predict_interval(test_df=test_df)

for model_name, (lower, point, upper) in interval_dict.items():
    print(f"{model_name}: coverage target = {1 - conformal_config.alpha:.0%}")

# 5. Evaluate interval quality
predictions_df, metrics_df = forecaster.evaluate_interval_forecast(test_df=test_df)
# metrics_df includes: picp (coverage), winkle-score, ace, nmpi, cwe

Calibration Flow#

        sequenceDiagram
    participant User
    participant TF as TwigaForecaster
    participant Conf as Conformal Factory
    participant Model

    User->>TF: calibrate(calibrate_df)
    TF->>Model: predict(calibrate_df)
    Model-->>TF: predictions
    TF->>TF: get_ground_truth()
    loop For each model
        TF->>Conf: Conformal(method, score_type, alpha)
        Conf-->>TF: conformal_instance
        TF->>Conf: calibrate(predictions, targets)
        Conf->>Conf: get_scores() → calculate_conformal_quantile()
        Conf-->>TF: calibrated (q_hat set)
        TF->>TF: store in self.conformal[model_name]
    end

API Reference#

twiga.distributions.conformal.base.calculate_conformal_value(scores, alpha, axis=0)#

Calculate the 1-alpha quantile of scores for conformal prediction using NumPy.

This function computes the threshold value (quantile) used to construct prediction sets based on the given non-conformity scores and significance level alpha. If the scores are empty or the quantile value exceeds 1, it returns an array of np.inf.

Parameters:

scores (ndarray) – Non-conformity scores with shape (N, …).
alpha (float) – Significance level, must be between 0 and 1.
axis (int) – Axis along which to compute the quantile. Defaults to 0.

Return type:

ndarray

Returns:

np.ndarray – The threshold value used to construct prediction sets, with shape (scores.shape[1:]).

Raises:

ValueError – If alpha is not between 0 and 1.

twiga.distributions.conformal.base.calculate_temporal_weights(timestamps, lambda_)#

Compute exponentially decaying temporal weights for calibration samples.

Weights are proportional to exp(-lambda_ * d_i) where d_i is the normalised age of sample i relative to the most recent calibration point (d_i = 0 for the newest sample, d_i = 1 for the oldest). The resulting weights are normalised to sum to 1.

A lambda_ of 0 recovers uniform weights. Larger values concentrate weight on more recent samples, reducing the influence of calibration samples that may have been drawn from a different seasonal regime.

Parameters:

timestamps (Series | ndarray) – Sequence of calibration timestamps (any datetime-like or numeric type). Must be 1-D with length equal to the number of calibration samples.
lambda_ (float) – Exponential decay rate. Must be >= 0.

Return type:

ndarray

Returns:

1-D array of normalised weights summing to 1, shape (N,).

Raises:

ValueError – If lambda_ is negative or timestamps is empty.

twiga.distributions.conformal.base.calculate_weighted_conformal_value(scores, alpha, weights)#

Compute the weighted (1-alpha) quantile of non-conformity scores.

Implements the weighted conformal prediction quantile of Tibshirani et al. (2019). For each horizon column, scores are sorted ascending and the smallest score whose cumulative normalised weight reaches 1 - alpha is returned as the threshold.

Parameters:

scores (ndarray) – Non-conformity scores with shape (N,) or (N, H).
alpha (float) – Significance level in (0, 1).
weights (ndarray) – Normalised sample weights with shape (N,). Must sum to 1.

Return type:

ndarray

Returns:

Threshold(s) with shape () (scalar) for 1-D input or (H,) for 2-D input. Returns np.inf for any horizon where the cumulative weight cannot reach 1 - alpha.

Raises:

ValueError – If alpha is outside (0, 1) or shapes are incompatible.

class twiga.distributions.conformal.base.BaseConformal(alpha=0.1, calib_method='uniform', lambda_=1.0, pooled=False)#

Bases: ABC

Abstract base class for Conformal Prediction in Regression.

Provides core functionality for constructing prediction intervals with finite-sample coverage guarantees.

Variables:

alpha (float) – Significance level for prediction intervals (0 < alpha < 1).
calib_method (str) – Quantile estimation method — "uniform" (default) or "temporal" (exponentially weighted by sample recency).
lambda (float) – Exponential decay rate for temporal weighting. Ignored when calib_method="uniform".
q_hat (float | np.ndarray) – Calibrated threshold(s) for intervals.

Parameters:

alpha (float) – Significance level (0 < alpha < 1). Defaults to 0.1.
calib_method (str) – "uniform" or "temporal". Defaults to "uniform".
lambda_ (float) – Decay rate for temporal weighting. Defaults to 1.0.
pooled (bool) – Pool non-conformity scores across the horizon axis before computing the calibration quantile, producing one global threshold instead of the default one-per-horizon. Requires calib_method="uniform". Defaults to False.

Raises:

ValueError – If alpha is not in (0, 1), calib_method is invalid, or pooled=True is combined with calib_method="temporal".

__init__(alpha=0.1, calib_method='uniform', lambda_=1.0, pooled=False)#

Initializes base conformal predictor with validation.

calculate_conformal_quantile(scores, axis=0, calib_timestamps=None, alpha=None)#

Computes the (1-alpha)-adjusted quantile of non-conformity scores.

Dispatches to the uniform empirical quantile (Lei et al. 2017) or the temporally weighted quantile (Tibshirani et al. 2019) depending on self.calib_method.

Parameters:

scores (ndarray) – Non-conformity scores array with shape (N,) or (N, H).
axis (int) – Quantile computation axis (uniform method only). Defaults to 0.
calib_timestamps (Series | ndarray | None) – Calibration sample timestamps required when calib_method="temporal". Ignored for "uniform".
alpha (float | None) – Significance level override. When None (default) self.alpha is used. Pass self.alpha / 2 for split-quantile calibration (e.g. sign-res score type) to compute each tail at half the nominal level, achieving valid 1-alpha joint coverage.

Return type:

float | ndarray

Returns:

Quantile values for interval construction.

Raises:

ValueError – For empty scores or missing timestamps in temporal mode.

calibrate(*calib_args, axis=0)#

Calibrates conformal thresholds using provided arguments.

Parameters:

*calib_args – Implementation-specific calibration data
axis (int) – Axis for quantile computation. Defaults to 0.

Raises:

ValueError – If calibration data validation fails

Return type:

None

abstractmethod generate_intervals(*pred_args)#

Generates prediction intervals from inputs.

Must be implemented by concrete subclasses.

Return type:: tuple[ndarray, ndarray]
Returns:: Tuple of (lower bounds, upper bounds)

abstractmethod get_scores(*args)#

Computes non-conformity scores from inputs.

Must be implemented by concrete subclasses.

Return type:: ndarray
Returns:: Array of non-conformity scores

class twiga.distributions.conformal.base.SplitConformal(score_type='res', alpha=0.1, calib_method='uniform', lambda_=1.0, pooled=False)#

Bases: BaseConformal

Implements residual-based conformal prediction for regression tasks.

Variables:: score_type (str) – Type of non-conformity score (‘res’ or ‘sign-res’).

generate_intervals(mu_pred)#

Generates prediction intervals.

Parameters:: mu_pred (ndarray) – Model predictions.
Return type:: tuple[ndarray, ndarray]
Returns:: tuple[np.ndarray, np.ndarray] – Lower and upper bounds.

get_scores(predicts, targets)#

Computes non-conformity scores based on residuals.

Parameters:

predicts (ndarray) – Model predictions.
targets (ndarray) – Ground truth targets.

Return type:

ndarray

Returns:

np.ndarray – Non-conformity scores.

class twiga.distributions.conformal.cqr.ConformalQuantileRegressor(score_type='scaled', alpha=0.1, target_coverage='marginal', *, pool_targets=None)#

Bases: BaseConformal

Conformal quantile regression for uncertainty estimation.

Implements conformal prediction intervals for quantile regression models using either scaled or unscaled non-conformity scores.

For multi-target quantile regression (lower/upper/targets shape (N, H, D)), three coverage modes are available via target_coverage:

"marginal" (default): calibrate each target independently. Each target’s interval covers that target with probability at least 1 - alpha. The joint probability that every target is covered is not controlled and is lower than 1 - alpha – roughly (1 - alpha)^D under independence, somewhere between (1 - alpha)^D and 1 - alpha under dependence.
"pooled": pool the per-target score by max across the trailing target axis before calibrating, giving one shared threshold per horizon (shape (H,) instead of (H, D)). Because max_d s_{i,h,d} <= q_hat[h] iff every term does, the resulting intervals cover all targets simultaneously with probability at least 1 - alpha. The cost is that the shared threshold is driven by the target whose normalised error is largest, so every other target’s interval is wider than it strictly needs to be. This is the same max-pooling construction as MultivariateConformal’s "max-norm" score, applied to CQR’s quantile-band score instead of a location/scale residual, and inherits the same valid-but-conservative trade-off.
"bonferroni": keep per-target calibration but reduce the effective significance level to alpha / D, where D is the number of targets. Each target is covered with probability at least 1 - alpha / D; by the union bound the joint probability that every target is covered is at least 1 - alpha. This is more conservative than "pooled" under strong negative dependence and less conservative under strong positive dependence. Neither "bonferroni" nor "pooled" is universally tighter; the choice depends on the dependence structure of the targets and the coverage statement the consumer requires.

In all three modes, score_type="scaled" keeps the per-target interval-width scaling factor upper - lower; only the conformal correction q_hat and its shape differ.

Variables:

score_type (str) – Non-conformity score type, "scaled" or "unscaled".
alpha (float) – Significance level for prediction intervals.
target_coverage (str) – One of "marginal", "pooled", "bonferroni".
pool_targets (bool) – Read-only derived property, target_coverage == "pooled". Retained for backward compatibility with code that read the old flag; setting it is not supported.
q_hat (float | np.ndarray) – Calibrated threshold(s). Shape depends on target_coverage: (H, D) for "marginal" and "bonferroni", (H,) for "pooled".
_EPS (float) – Numerical stability constant.

Parameters:

score_type (str) – Score calculation method. Defaults to "scaled".
alpha (float) – Significance level (0 < alpha < 1). Defaults to 0.1.
target_coverage (str) – Coverage mode. "marginal" (default) gives per-target coverage; "pooled" gives simultaneous coverage via max-pooling; "bonferroni" gives simultaneous coverage via a union-bound correction. See the class docstring for the trade-offs.

Raises:

ValueError – For invalid score_type, alpha, or target_coverage values.

Conformal Prediction & Uncertainty Quantification#

How It Works#

Class Hierarchy#

Factory: Conformal#

Configuration#

Valid Score Types per Method#

Methods#

Split Conformal (method="residual")#

Conformal Quantile Regression (method="quantile")#

Conformal Residual Fitting (method="residual-fitting")#

Temporal Calibration Weighting#

Helper Functions#

Integration with TwigaForecaster#

Calibration Flow#

API Reference#

Factory: `Conformal`#

Split Conformal (`method="residual"`)#

Conformal Quantile Regression (`method="quantile"`)#

Conformal Residual Fitting (`method="residual-fitting"`)#