Distributional Metrics¶

← Interval Metrics | Calibration Metrics →

Distributional metrics evaluate probabilistic forecasts — how well does the predicted probability distribution \(F\) match the true (but unknown) data-generating process \(G\)? Unlike point metrics, these assess both calibration (reliability) and sharpness (precision).

Proper Scoring Rules¶

A scoring rule \(S(F, y)\) is proper if its expected value is minimised when the predicted distribution \(F\) is equal to the true distribution \(G\) (Ref. 1).

In torchregress, we prioritise proper scoring rules for evaluating all probabilistic models.

Continuous Ranked Probability Score (CRPS)¶

The CRPS is the most widely used proper scoring rule for univariate regression. It can be viewed as the integral of the pinball loss over all possible quantiles \(\tau \in [0, 1]\).

Properties¶

• Units: Same as the target variable \(y\).
• Point Mass: Reduces to Mean Absolute Error (MAE) if \(F\) is a point mass.
• Duality: Simultaneously rewards calibration (is the truth within the predicted range?) and sharpness (is the predicted range narrow?).

Implementation¶

from torchregress.metrics import crps_gaussian, energy_score

# For Gaussian models — note: argument order is (mean, y_true, std)
loss = crps_gaussian(mu, y_true, sigma)

# For non-parametric models (e.g., Ensembles, BNNs) using samples
loss = energy_score(y_samples, y_true)

→ See Mathematical Foundations for the Gaussian closed-form derivation. API Reference: crps_gaussian.

Estimators and conventions¶

The fair estimator is unbiased for the CRPS of the distribution the members are drawn from. The biased \(1/M^2\) version (properscoring.crps_ensemble, scoringrules estimator="nrg") favours small ensembles. One sample is a Dirac forecast, so the CRPS is the absolute error. The quantile form only approximates the integral: it has no information beyond the outermost levels, so use the sample or closed-form estimator when you have samples or a distribution.

vario_score / VarioScore (Zamo & Naveau) is positively oriented, \(\nu_\rho = \tfrac12\mathbb E\lvert X - X'\rvert^\rho - \mathbb E\lvert X - y\rvert^\rho\) (equal to \(-\text{CRPS}\) at \(\rho = 1\)), so VarioScore.higher_is_better is True. variogram_score sums over pairs \(i < j\), which is half of scoringrules.variogram_score. All reference conventions are pinned in tests/metrics/test_reference_parity.py.

Multivariate: Energy Score¶

The Energy Score (ES) (Ref. 2) is the multivariate generalisation of CRPS to \(\mathbb{R}^d\). It evaluates the joint distribution of multiple targets, capturing correlations that univariate CRPS misses.

where \(\beta \in (0, 2)\) (default is \(\beta=1\)).

Implementation¶

from torchregress.metrics import energy_score

# y_samples: [num_samples, batch_size, num_targets]
score = energy_score(y_samples, y_true)

API Reference: energy_score.

Calibration: Probability Integral Transform (PIT)¶

A model is perfectly calibrated if its predictive CDF \(F(y \mid x)\), when evaluated at the true value \(y\), is uniformly distributed on \([0, 1]\) (Ref. 3).

Diagnosing Miscalibration¶

• U-Shaped: The model is overconfident (true values fall in the tails too often).
• Hump-Shaped: The model is underconfident (true values fall in the center too often).
• Skewed: The model has a consistent bias (predicting too high or too low).

Implementation¶

To visualise calibration, use the PIT Histogram diagnostic from the visualization module:

from torchregress.viz import plot_pit_histogram

# Generate a PIT histogram to visualize calibration for Gaussian predictions
plot_pit_histogram(y_pred, y_pred_std, y_true, n_bins=20)

API Reference: plot_pit_histogram.

Highest Posterior Density (HPD) Coverage¶

For a 1D predictive density \(p(y \mid x)\) on a grid, the HPD level of an observation \(y\) is the mass of the smallest region that contains \(y\) and consists of the highest-density points:

For a calibrated predictive density, \(\ell(y)\) is uniform on \([0, 1]\). highest_posterior_density_level returns \(\ell\) for each target, and highest_posterior_density_coverage returns the fraction of targets with \(\ell(y) \le \alpha\).

import torch
from torchregress.metrics import highest_posterior_density_coverage

support = torch.linspace(-6.0, 6.0, 1201)
mu = torch.randn(2000)  # predicted means
density = torch.exp(-0.5 * (support[None, :] - mu[:, None]) ** 2)  # unit-variance Gaussians (unnormalised)
y_true = mu + torch.randn(2000)  # calibrated: targets drawn from the predictions

cov_90 = highest_posterior_density_coverage(support, density, y_true, alpha=0.9)  # approx. 0.90
cov_10 = highest_posterior_density_coverage(support, density, y_true)  # alpha=0.1 default: approx. 0.10

Sample-based (multivariate) HPD coverage¶

A grid does not scale past one or two targets. For any model that can evaluate \(\log p(y \mid x)\) and draw \(y_s \sim p(\cdot \mid x)\) (conditional flows, DensityEnsemble), the HPD level is estimated by ranking the observed density among the draws:

If \(y\) is exchangeable with the draws, \(\hat\ell\) is uniform on \([0, 1]\) up to \(1/S\) discretisation, in any target dimension. This tests the joint region, so it detects wrong inter-target correlations that per-dimension PIT misses. sample_hpd_coverage takes credibility as the nominal region mass (0.9 = 90% region).

import torch
from torchregress.metrics import sample_hpd_coverage

dist = torch.distributions.MultivariateNormal(torch.zeros(3), torch.eye(3))
y_true = dist.sample((2000,))
draws = dist.sample((2000, 256))  # [N, S, D] predictive draws
coverage = sample_hpd_coverage(dist.log_prob(y_true), dist.log_prob(draws))
# approx. tensor([0.50, 0.68, 0.90, 0.95])

Unified Metrics Report¶

For comprehensive evaluation, use the distribution_metrics_report helper. It consolidates NLL, CRPS, Energy Score, PIT uniformity, and coverage into a single dictionary.

from torchregress.metrics import distribution_metrics_report

# dist: torch.distributions.Distribution
# y_true: Ground truth tensor
results = distribution_metrics_report(dist=dist, y_true=y_true)

print(f"CRPS: {results['crps']:.4f}")
print(f"PIT KS: {results['pit_ks']:.4f}")
print(f"90% Coverage: {results['coverage_90']:.2%}")

This is the recommended way to evaluate complex probabilistic models, as it provides a multi-faceted view of model performance.

results["crps"] uses the best representation available. A Normal dist gets the closed form. samples (or samples drawn from a non-Normal dist) get the fair ensemble estimator. Only y_pred_quantiles use the quantile approximation. Before 0.3 the report always reduced samples to 7 quantiles, which underestimated the CRPS by 3-4%.

Binned-CDF and PIT helpers¶

For discrete probability masses on fixed bin edges, binned_cdf converts the masses to CDF values at the edges and crps_from_binned_cdf evaluates the exact piecewise-linear CRPS:

pit_cvm is the Cramér–von Mises distance of PIT values from Uniform(0, 1), and paired_bootstrap_ci returns a paired percentile interval for the mean difference of two metric vectors. These are generic calibration/diagnostic helpers; binned mass is not a continuous density and its cross-entropy is not a continuous NLL. CRPS includes the distance from a reference target outside the finite support; PIT values and bin masses must be finite and valid. These reporting helpers return CPU scalars/NumPy arrays, including for tensor inputs.

import numpy as np
from torchregress.metrics import binned_cdf, crps_from_binned_cdf, pit_cvm

probabilities = np.array([[0.2, 0.5, 0.3]])
edges = np.array([0., 1., 2., 3.])
y = np.array([1.2])
pit_values = np.array([0.2, 0.5, 0.8])
cdf = binned_cdf(probabilities)                 # [N, K+1]
score = crps_from_binned_cdf(edges, cdf, y)     # exact binned CRPS
cvm = pit_cvm(pit_values, normalized=True)

Summary Matrix¶

References¶

Limitations¶

1. CRPS sample requirements: crps_from_samples requires \(S \ge 50\) MC samples per test point for reliable estimates. Using fewer samples introduces Monte Carlo noise that can obscure model comparisons. For Gaussian models, use the analytic crps_gaussian instead — it is exact and faster.
2. Energy Score is \(\mathcal{O}(S^2)\): Computing the pairwise-expectation term requires evaluating \(\|Y - Y'\|^\beta\) for all \(S(S-1)/2\) sample pairs. For \(S > 500\), this becomes a runtime bottleneck.
3. PIT is necessary but not sufficient: A uniform PIT histogram indicates marginal calibration but does not guarantee conditional calibration \(F(y \mid X = x)\). A model can pass the PIT test while being poorly calibrated for specific subpopulations.
4. NLL sensitivity: NLL is dominated by low-probability events. A single point with \(p(y \mid x) \approx 0\) produces an arbitrarily large NLL contribution. CRPS is more robust to these tail events.

Recommendations¶

• Default suite: Report CRPS (interpretable, robust) + NLL (sensitive to tail calibration) + PIT (calibration diagnostic) for a complete picture.
• For Gaussian models: Always use crps_gaussian (analytic, exact) over crps_from_samples.
• For multivariate targets: Use energy_score (\(\beta=1\)) as the multivariate CRPS analogue.
• PIT bin count: Use 20–50 bins for the PIT histogram. Too few bins hide miscalibration patterns; too many create noise.
• distribution_metrics_report consolidates CRPS, NLL, PIT, and coverage into a single call.

Next Steps¶

• Learn about Calibration Metrics
• View the Distributional Conformal Tutorial
• Explore Normalizing Flow Examples