A Large-Scale Study of Probabilistic Calibration in Neural Network
  Regression
Accurate probabilistic predictions are essential for optimal decision making.
While neural network miscalibration has been studied primarily in
classification, we investigate this in the less-explored domain of regression.
We conduct the largest empirical study to date to assess the probabilistic
calibration of neural networks. We also analyze the performance of
recalibration, conformal, and regularization methods to enhance probabilistic
calibration. Additionally, we introduce novel differentiable recalibration and
regularization methods, uncovering new insights into their effectiveness. Our
findings reveal that regularization methods offer a favorable tradeoff between
calibration and sharpness. Post-hoc methods exhibit superior probabilistic
calibration, which we attribute to the finite-sample coverage guarantee of
conformal prediction. Furthermore, we demonstrate that quantile recalibration
can be considered as a specific case of conformal prediction. Our study is
fully reproducible and implemented in a common code base for fair comparisons.
Victor Dheur, Souhaib Ben Taieb
# A Large-Scale Study of Probabilistic Calibration

###### Abstract

Accurate probabilistic predictions are essential for optimal decision making. While neural network miscalibration has been studied primarily in classification, we investigate this in the less-explored domain of regression. We conduct the largest empirical study to date to assess the probabilistic calibration of neural networks. We also analyze the performance of recalibration, conformal, and regularization methods to enhance probabilistic calibration. Additionally, we introduce novel differentiable recalibration and regularization methods, uncovering new insights into their effectiveness. Our findings reveal that regularization methods offer a favorable tradeoff between calibration and sharpness. Post-hoc methods exhibit superior probabilistic calibration, which we attribute to the finite-sample coverage guarantee of conformal prediction. Furthermore, we demonstrate that quantile recalibration can be considered as a specific case of conformal prediction. Our study is fully reproducible and implemented in a common code base for fair comparisons.

Machine Learning, Probabilistic Calibration, Probabilistic Calibration

## 1 Introduction

Neural network predictions affect critical decisions in many applications, including medical diagnostics and autonomous driving (Gulshan et al., 2016; Guizilini et al., 2020). However, effective decision making often requires accurate probabilistic predictions (Gawlikowski et al., 2021; Abdar et al., 2021). For example, consider a probabilistic regression model that produces 90% prediction intervals. An important property would be that 90% of these prediction intervals contain the realizations.

For models that output a predictive distribution, _probabilistic calibration_ is an important property that states that all quantiles must be calibrated, i.e., the frequency of realizations below these quantiles must match the corresponding quantile level. Additionally, predictive distributions should be sufficiently sharp (i.e., concentrated around the realizations) and leverage the information in the inputs.

In the classification setting, Guo et al. (2017) found that common neural architectures trained on image and text data were miscalibrated, sparking increased interest in neural network calibration. In a follow-up study, Minderer et al. (2021) showed that more recent neural architectures demonstrate improved calibration. However, there has been less research on calibration for neural probabilistic regression models compared to classification. Therefore, it remains uncertain whether the same results apply to the regression setting. This paper addresses this gap by conducting a comprehensive study on probabilistic calibration for regression using tabular data. We explore various calibration methods, including quantile recalibration (Kuleshov, Fenner, et al., 2018) and conformalized quantile regression (Romano, Patterson, et al., 2019). We also consider regularization methods, which have been shown to perform well in the classification setting (Karandikar et al., 2021; Popordanoska et al., 2022; Yoon et al., 2023).

We make the following main contributions:

1. We conduct the largest empirical study to date on probabilistic calibration of neural regression models using 57 tabular datasets (Sections 4 and 6). We consider multiple state-of-the-art calibration methods (Section 5), including post-hoc recalibration, conformal prediction, and regularization methods, with various scoring rules and predictive models.
2. Building on quantile recalibration, we propose a new differentiable calibration map using kernel density estimation, which provides improved negative log-likelihood compared to baselines. We also introduce two new regularization objectives based on the probabilistic calibration error (Section 5).
3. We show that quantile recalibration is a special case of conformal prediction, providing an explanation for
their superior performance in terms of probabilistic calibration (Section 6).

## 2 Background

We consider a univariate regression problem where the target variable \(Y\in\mathcal{Y}\) depends on an input variable \(X\in\mathcal{X}\), with \(\mathcal{Y}=\mathbb{R}\) representing the target space and \(\mathcal{X}\) representing the input space. Our objective is to approximate the conditional distribution \(P_{Y|X}\) using training data \(\mathcal{D}=\{\,(X_{i},Y_{i})\,\}_{i=1}^{N}\) where \((X_{i},Y_{i})\stackrel{{\text{i.i.d.}}}{{\sim}}P\equiv P_{X}\times P _{Y|X}\).

A probabilistic predictor \(F_{\theta}:\mathcal{X}\rightarrow\mathcal{F}\) is a function parametrized by \(\theta\in\Theta\) that maps an input \(x\in\mathcal{X}\) to a predictive cumulative distribution function (CDF) \(F_{\theta}(\cdot\mid x)\) in the space \(\mathcal{F}\) of distributions over \(\mathbb{R}\). Additionally, given \(x\in\mathcal{X}\), we denote the predictive quantile function (QF) by \(Q_{\theta}(\cdot\mid x)\), and probability density function (PDF) by \(f_{\theta}(\cdot\mid x)\). Similarly, the marginal CDF, QF, or PDF of a random variable \(R\) is denoted by \(F_{R}\), \(Q_{R}\), or \(f_{R}\), respectively.

Probabilistic calibration.Given an input \(x\in\mathcal{X}\), the model \(F_{\theta}\) is ideal if it precisely matches the conditional distribution \(P_{Y|X}\). However, learning the ideal model based on finite data is not possible without additional (strong) assumptions (Foygel Barber et al., 2021). To avoid additional assumptions, we can instead enforce certain desirable properties that are attainable in practice and that a good or ideal forecaster should exhibit. One such property is probabilistic calibration.

Let \(Z=F_{\theta}(Y\mid X)\in[0,1]\) denote the probability integral transform (PIT) of \(Y\) conditional on \(X\). The model \(F_{\theta}\) is _probabilistically calibrated_ (also known as PIT-calibrated) if \(\forall\alpha\in[0,1]\),

\[F_{Z}(\alpha)\doteq\text{Pr}(Z\leq\alpha)=\alpha. \tag{1}\]

Let \(U\in[0,1]\) be a uniform random variable independent of \(Z\). The left and right hand sides of (1) can be interpreted as the CDF of \(Z\) and \(U\), respectively, as a function of \(\alpha\). This shows that the uniformity of the PIT is equivalent to probabilistic calibration (Dawid, 1984).

Since the ideal forecaster is probabilistically calibrated, we can require this property from any competent forecaster. However, probabilistic calibration, though necessary, is not sufficient for making accurate probabilistic predictions. Additionally, as discussed by Gneiting and Resin (2021), probabilistic calibration primarily addresses unconditional aspects of predictive performance and is implied by more robust conditional notions of calibration, such as auto-calibration.

Probabilistic calibration error.The most common approach for evaluating probabilistic calibration is to consider distances of the form \(\int_{0}^{1}|F_{Z}(\alpha)-F_{U}(\alpha)|^{p}d\alpha\) where \(p>0\). The particular cases of \(p=1\) and \(p=2\) are known as the 1-Wasserstein distance and Cramer-von Mises distance, respectively. We denote the empirical CDF of the PIT as \(\hat{F}_{Z}(\alpha)=\frac{1}{N}\sum_{i=1}^{N}\mathbb{1}(Z_{i}\leq\alpha)\) where \(Z_{i}=F_{\theta}(Y_{i}\mid X_{i})\) are PIT realizations. A common approach to assess probabilistic calibration using Monte Carlo estimation is to evaluate it at equidistant values \(\alpha_{1}<\dots<\alpha_{M}\) as follows:

\[\text{PCE}_{p}(F_{\theta},\mathcal{D})=\frac{1}{M}\sum_{j=1}^{M}\left|\alpha_ {j}-\hat{F}_{Z}(\alpha_{j})\right|^{p}. \tag{2}\]

This metric has been previously employed in literature such as Zhao et al. (2020) and Zhou et al. (2021) with \(p=1\), and Kuleshov, Fenner, et al. (2018) and Utpala and Rai (2020) with \(p=2\). It is important to note that, unlike the classical definition of the \(p\)-norm, we do not exponentiate \(\frac{1}{p}\) in (2) to maintain consistency with prior literature. In the subsequent sections, we focus our analysis on \(\text{PCE}_{1}\) and use the abbreviation PCE for brevity.

One limitation of scalar metrics like PCE is their inability to provide detailed information regarding calibration errors at individual quantile levels, \(\alpha_{1},\dots,\alpha_{M}\). Instead, PIT reliability diagrams offer a visual assessment of probabilistic calibration across all quantile levels by plotting the empirical CDF of the PIT \(Z\). These diagrams display the right side of (1) against its left side, with a perfectly calibrated model represented by a diagonal line (in the asymptotic case). Figure 2 provides examples of such reliability diagrams, which have been employed in studies by Pinson and Hagedorn (2012) and Kuleshov, Fenner, et al. (2018).

## 3 Related Work

**Post-hoc calibration approaches** involve adjusting the predictions of a trained model using a mapping learned from a separate calibration dataset. In the context of classification, temperature scaling (Guo et al., 2017) is a simple and effective method that adjusts predictive confidence while maintaining accuracy. For regression tasks, quantile recalibration (Kuleshov, Fenner, et al., 2018) aims to achieve probabilistic calibration. Conformal prediction (Vovk et al., 2020) is a general approach that provides prediction sets with a finite-sample coverage guarantee. Notable methods applied with deep learning include Conformal Quantile Regression (Romano, Patterson, et al., 2019) and Distributional Conformal Prediction (Izbicki et al., 2020; Chernozhukov et al., 2021). Furthermore, post-hoc approaches have also been proposed for conditional notions of calibration (Song et al., 2019; Kuleshov and Deshpande, 2022).

**Regularization approaches** aim to improve calibration during training by incorporating regularization techniques. Some methods, proposed by Zhao et al. (2020) and Feldman et al. (2021), utilize regularization to target different
conditional notions of calibration based on the inputs. Zhou et al. (2021) introduced an alternative loss function involving the simultaneous training of two neural networks, while Pearce et al. (2018), Chung et al. (2021), and Thiagarajan et al. (2020) proposed objectives that allow control over the tradeoff between coverage and sharpness of prediction intervals. To our knowledge, the only regularization objective specifically targeting probabilistic calibration is quantile regularization (Utpala and Rai, 2020). Other types of uncertainty quantification methods include ensembling (Lakshminarayanan, Pritzel, et al., 2017) and Bayesian methods (Jospin et al., 2022)

## 4 Are Neural Regression Models

Probabilistically Calibrated?

We conduct an extensive empirical study to evaluate the probabilistic calibration of neural regression models. To this end, we calculate the _probabilistic calibration error_ defined in (2) for various state-of-the-art models across multiple benchmark datasets.

**Benchmark datasets**. We analyze a total of 57 datasets, including 27 from the OpenML curated benchmark (Grinsztajn et al., 2022), 18 from the AutoML Repository (Gijsbers et al., 2019), and 12 from the UCI Machine Learning Repository (Dua and Graff, 2017).

These datasets are widely used in the evaluation of deep probabilistic models and uncertainty quantification, as evidenced by previous studies such as Fakoor et al. (2021), Chung et al. (2021), Zhou et al. (2021), Utpala and Rai (2020), and Gal and Ghahramani (2016).

Figure 1 provides an overview of the utilization of these datasets in previous studies. To the best of our knowledge, our study represents the most comprehensive assessment of probabilistic calibration for neural regression models published to date.

**Neural probabilistic regression models**. We consider three state-of-the-art neural probabilistic regression models. The first model predicts a parametric distribution, where the parameters are obtained as outputs of a hypernetwork. Previous studies have often focused on the Gaussian distribution (Lakshminarayanan, Pritzel, et al., 2017; Utpala and Rai, 2020; Zhao et al., 2020). To introduce more flexibility, we consider a mixture of \(K\) Gaussian distributions. Given an input \(x\in\mathcal{X}\), the hypernetwork parametrizes the means \(\mu_{k}(x)\), standard deviations \(\sigma_{k}(x)\), and weights \(w_{k}(x)\) for each component \(k=1,...,K\). To ensure positive standard deviations and that the mixture weights form a discrete probability distribution, we use the Softplus and Softmax activations, respectively. We have two variants of this model depending on the scoring rule used for training: the negative log-likelihood (NLL) or the continuous ranked probability score (CRPS). These models are denoted as MIX-NLL and MIX-CRPS, respectively. It is worth noting that the CRPS of a mixture of Gaussians has a closed-form expression (Grimit et al., 2006).

The second model predicts quantiles of the distribution (Tagasovska and Lopez-Paz, 2019; Chung et al., 2021; Feldman et al., 2021). Specifically, given an input \(x\in\mathcal{X}\) and a quantile level \(\alpha\in[0,1]\), the model outputs a quantile \(Q_{\theta}(\alpha\mid x)\). The full quantile function can be obtained by evaluating the model at multiple quantile levels. The model is trained by minimizing the quantile score at multiple levels, which is asymptotically equivalent to minimizing the CRPS (Bracher et al., 2021). We denote this model as SQR-CRPS, where SQR stands for simultaneous quantile regression (Tagasovska and Lopez-Paz, 2019).

**Experimental setup**. We adopt the large-sized regime introduced by Grinsztajn et al. (2022), which involves truncating the datasets to a maximum of 50,000 examples. Among the 57 datasets, the number of examples ranges from 135 to 50,000, and the number of features ranges from 3 to 3,6111. Each of the 57 datasets is divided into four sets: training (65%), validation (10%), calibration (15%), and test (10%). We normalize the input \(X\) and target \(Y\) using

Figure 1: Multiple regression benchmark datasets with references. Datasets inside parentheses have not been considered in this study. Full dataset names are available in Table 1.
the mean and standard deviation from the training split. The final predictions are then transformed back to the original scale. For our neural network models, we employ the same fully-connected architecture as previous studies conducted by Kuleshov, Fenner, et al. (2018), Chung et al. (2021), and Fakoor et al. (2021). Further details regarding the model hyperparameters can be found in Appendix C.

**Results.** In Figure 2, the first row displays the PCE (averaged over five random train-validation-test splits) for MIX-NLL in blue on each of the 57 datasets. For comparison, the PCE of a perfectly calibrated model, i.e. with uniformly distributed PITs, computed using \(5\times 10^{4}\) simulated values is shown in orange. The second row presents reliability diagrams for five datasets, with 90% consistency bands as in Gneiting, Wolffram, et al. (2023). Similar information is provided for MIX-CRPS and SQR-CRPS in Figures 12 and 13 in Appendix B.4, respectively. Additionally, reliability diagrams for all datasets can be found in Figure 15 in Appendix B.6.

The analysis reveals that the (average) PCE is generally high across many datasets, although there are significant variations between datasets. To test the statistical significance of these results, \(10^{4}\) samples were generated from the sampling distribution of the average PCE under the null hypothesis of probabilistic calibration. The resulting sampling distribution for all datasets is presented in Appendix B.5.

By computing the p-value associated with a one-sided test in the upper tail of the distribution (as illustrated in Appendix B.5), it was observed that most datasets have a p-value of zero. This indicates that the average PCE obtained for the considered model is higher than all the simulated average PCEs of the probabilistically calibrated model. Applying a threshold of 0.01 and a Holm correction for the 57 hypothesis tests, the null hypothesis is rejected for 11 datasets out of the 57.

Overall, the results indicate that the neural models considered in this study are generally probabilistically mis-calibrated on a significant number of benchmark tabular datasets. In Section 6, we will further explore how calibration methods can substantially improve the PCE of neural models.

## 5 Calibration Methods

We begin by discussing the three main approaches to calibration: quantile recalibration, conformal prediction, and regularization-based calibration. Following that, we introduce two novel variants of regularization-based calibration.

Quantile recalibration and conformal prediction are post-hoc methods, meaning they are applied after model training. These approaches utilize a separate calibration dataset \(\mathcal{D}^{\prime}=\{\ (X_{i}^{\prime},Y_{i}^{\prime})\ \}_{i=1}^{N^{\prime}}\), where \((X_{i}^{\prime},Y_{i}^{\prime})\stackrel{{\text{i.i.d.}}}{{\sim}}P _{X,Y}\). On the other hand, regularization-based calibration operates directly during training and relies solely on the training data \(\mathcal{D}\).

### Quantile Recalibration

Quantile recalibration aims to transform a potentially mis-calibrated CDF \(F_{\theta}\) into a probabilistically calibrated CDF \(F_{\theta}^{\prime}=F_{Z}\circ F_{\theta}\), using the calibration map \(F_{Z}\) which repre

Figure 2: The top row shows the PCE for different datasets with one standard error (error bar). The bottom row gives examples of PIT reliability diagrams for five datasets.
sents the CDF of the PITs for \(F_{\theta}\). For a given quantile level \(\alpha\in[0,1]\), the recalibrated CDF \(F^{\prime}_{\theta}\) satisfies:

\[\text{Pr}(F^{\prime}_{\theta}(Y\mid X)\leq\alpha) =\text{Pr}(F_{\theta}(Y\mid X)\leq Q_{Z}(\alpha)) \tag{3}\] \[=F_{Z}(Q_{Z}(\alpha))\] (4) \[=\alpha. \tag{5}\]

In practice, \(F_{Z}\) is not directly available and needs to be estimated from data. Kuleshov et al. (2018) proposed estimating it using isotonic regression, while Utpala and Rai (2020) showed that computing the empirical CDF is an equivalent and simpler method. Specifically, given a set of PIT values \(Z^{\prime}_{i}=F_{\theta}(Y^{\prime}_{i}\mid X^{\prime}_{i}),i=1,\ldots,N^{\prime}\), the calibration map \(\phi^{\text{EMP}}\) is computed as:

\[\phi^{\text{EMP}}(\alpha;\{\,Z^{\prime}_{i}\}_{i=1}^{N^{\prime}})=\frac{1}{N^ {\prime}}\sum_{i=1}^{N^{\prime}}\mathbbm{1}(Z^{\prime}_{i}\leq\alpha), \tag{6}\]

where \(\alpha\in[0,1]\).

Similarly to Utpala and Rai (2020), we also consider a linear calibration map \(\phi^{\text{LIN}}\), which is continuous, and corresponds to a linear interpolation between the points \(\Big{\{}\;(0,0),(Z^{\prime}_{(1)},\nicefrac{{1}}{{N^{\prime}+1}}),\ldots,(Z^{ \prime}_{(N^{\prime})},\nicefrac{{N^{\prime}}}{{N^{\prime}+1}}),(1,1)\; \Big{\}}\), where \(Z^{\prime}_{(k)}\) is the \(k\)th order statistic of \(Z^{\prime}_{1},\ldots,Z^{\prime}_{N^{\prime}}\).

In addition, we propose a calibration map based on kernel density estimation (KDE), denoted as \(\phi^{\text{KDE}}\). This calibration map offers the advantage of being differentiable and can lead to improved NLL performance. The key idea is to use a relaxed approximation of the indicator function, which allows us to make the PIT CDF (6) differentiable. Specifically, we compute

\[\mathbbm{1}_{\tau}(a\leq b)=\sigma(\tau(b-a))\approx\mathbbm{1}(a\leq b),\]

where \(\tau>0\) is an hyperparameter and \(\sigma(x)=\frac{1}{1+e^{-x}}\) denotes the sigmoid function. The resulting smoothed empirical CDF is given by

\[\phi^{\text{KDE}}(\alpha;\{\,Z^{\prime}_{i}\,\}_{i=1}^{N^{\prime}})=\frac{1}{ N^{\prime}}\sum_{i=1}^{N^{\prime}}\mathbbm{1}_{\tau}(Z^{\prime}_{i}\leq\alpha). \tag{7}\]

This corresponds to estimating the CDF \(F_{Z}\) using KDE based on \(N^{\prime}\) realizations of \(Z\) (\(\{\,Z^{\prime}_{i}\,\}_{i=1}^{N^{\prime}}\)). Since \(\sigma\) is the CDF of the logistic distribution, we use the PDF of the logistic distribution as the kernel in the KDE. Algorithm 1 summarises this method.

```
Input: Predictive CDF \(F_{\theta}\) and \(\mathcal{D}^{\prime}=\{\;(X^{\prime}_{i},Y^{\prime}_{i})\;\}_{i=1}^{N^{\prime}}\).  Compute \(Z^{\prime}_{i}=F_{\theta}(Y^{\prime}_{i}\mid X^{\prime}_{i})\quad(i=1,\ldots,N ^{\prime})\)  Compute the calibration map \(\phi\), either \(\phi^{\text{EMP}}\), \(\phi^{\text{LIN}}\), or \(\phi^{\text{KDE}}\) Return: Recalibrated CDF \(F^{\prime}_{\theta}=\phi\circ F_{\theta}\).
```

**Algorithm 1** Quantile recalibration

### Conformal Prediction

Let us assume the realizations of our calibration dataset \(\mathcal{D}^{\prime}\) are drawn exchangeably from \(P_{X,Y}\)2. Given a predictive model \(M_{\theta}\) and a coverage level \(\alpha\in[0,1]\), (inductive) conformal prediction allows us to construct a prediction set \(C_{\alpha}(X)\subseteq\mathcal{Y}\) for any input \(X\), satisfying the property:

Footnote 2: This is implied by the common i.i.d. assumption.

\[\text{Pr}(Y\in C_{\alpha}(X)) =\frac{\lceil(N^{\prime}+1)\alpha\rceil}{N^{\prime}+1} \tag{8}\] \[\approx\alpha. \tag{9}\]

Conformal prediction achieves this by utilizing a conformity score \(s_{\theta}(Y\mid X)\), which intuitively quantifies the similarity between new samples and previously observed samples. When the conformity score increases with \(Y\), an interval \(C_{\alpha}(X)=(-\infty,s_{\theta}^{-1}(\alpha\mid X)]\) can be constructed, ensuring the conformal guarantee (8) at level \(\alpha\).

Let \(Q^{\prime}_{\theta}(\alpha\mid X)=s_{\theta}^{-1}(\alpha\mid X)\) represent the (revised) model obtained through conformal prediction from \(Q_{\theta}(\alpha\mid X)\). Under the assumption that \(Q^{\prime}_{\theta}\) is continuous and strictly increasing, the conformal guarantee implies that \(\text{Pr}(Y\leq Q^{\prime}_{\theta}(\alpha\mid X))\approx\alpha\), which indicates approximate probabilistic calibration at level \(\alpha\).

Conformalized Quantile Regression (Romano et al., 2019) is an example of a conformal procedure, where the conformity score is defined as \(s_{\theta}(Y\mid X)=Y-Q_{\theta}(\alpha\mid X)\), representing the quantile residual. Another example is Distributional Conformal Prediction (Izbicki et al., 2020; Chernozhukov et al., 2021), which employs the conformity score \(s_{\theta}(Y\mid X)=F_{\theta}(Y\mid X)\), referring to the PIT. Algorithm 2 provides a summary of how to compute calibrated quantiles using inductive conformal prediction.

```
Input: Trained model \(M_{\theta}\), \(\mathcal{D}^{\prime}=\{\;(X^{\prime}_{i},Y^{\prime}_{i})\;\}_{i=1}^{N^{\prime}}\), strictly increasing conformity score \(s\), quantile level \(\alpha\in[0,1]\), input \(X\).  Compute \(S_{i}=s_{\theta}(Y^{\prime}_{i}\mid X^{\prime}_{i})\quad(i=1,\ldots,N^{\prime})\)  Compute \(\hat{q}=S_{(\lceil(N^{+}1)\alpha\rceil)}\) where \(S_{(k)}\) denote the \(k\)th smallest value among \(\{\;S_{1},\ldots,S_{N^{\prime}},+\infty\;\}\) Return: Calibrated quantile \(Q^{\prime}_{\theta}(\alpha\mid X)=s_{\theta}^{-1}(\hat{q}\mid X)\)
```

**Algorithm 2** Calibrated quantiles with conformal prediction

### Regularization-based Calibration

Regularization-based calibration methods aim to enhance model calibration by incorporating a regularization term into the training objective. Although widely used in classification, there are relatively fewer methods specifically designed for regression problems. In this section, we discuss two approaches: quantile regularization (Utpala and Rai, 2020)
and the truncation method (Chung et al., 2021). The main steps of regularization-based calibration are summarized in Algorithm 3.

```
Input: Model \(M_{\theta}\), calibration regularizer \(\mathcal{R}(\theta)\) and tuning parameter \(\lambda\geq 0\).  Compute \(\theta^{*}=\arg\min_{\theta\in\Theta}\mathcal{L}^{{}^{\prime}}(\theta;\mathcal{D})\) where \(\mathcal{L}^{{}^{\prime}}(\theta;\mathcal{D})=\nicefrac{{1}}{{N}}\sum_{i=1}^{N }\mathcal{L}(M_{\theta}(\cdot\mid X_{i}),Y_{i})+\lambda\mathcal{R}(\theta; \mathcal{D})\) Return: Regularized model \(M_{\theta}\).
```

**Algorithm 3** Regularization-based calibration

#### 5.3.1 Quantile Regularization

The regularizer proposed by Utpala and Rai (2020) aims to measure the deviation of the PIT variable \(Z\) from a uniform distribution, which is characteristic of a probabilistically calibrated model. This regularization penalty encourages the selection of calibrated models during training.

The authors observed that the KL divergence between \(Z\) and a uniform random variable is equivalent to the negative differential entropy of \(Z\), denoted as \(H(Z)\). To approximate \(H(Z)\), they employed sample-spacing entropy estimation (Vasicek, 1976), resulting in the following regularizer:

\[\mathcal{R}_{\text{QR}}(\theta;\mathcal{D}) \tag{10}\] \[=\frac{1}{N-k}\sum_{i=1}^{N-k}\log\left[\frac{N+1}{k}(Z_{(i+k)}-Z _{(i)})\right]\] (11) \[\approx H(Z), \tag{12}\]

where \(k\) is a hyperparameter satisfying \(1\leq k\leq N\), and \(Z_{(i)}\) represents the \(i\)th order statistic of \(Z\).

To ensure differentiability during optimization, the authors employed a differentiable relaxation technique called NeuralSort (Grover et al., 2019), as sorting is a non-differentiable operation.

#### 5.3.2 Truncation-based Calibration

The regularization approach introduced by Chung et al. (2021), that we denote Trunc, involves truncating the predictive distribution based on the current level of calibration.

Given a quantile model \(Q_{\theta}\), let \(\hat{F}_{Z}(\alpha)=\frac{1}{N}\sum_{i=1}^{N}\mathbbm{1}(Y_{i}\leq Q_{\theta}( \alpha\mid X_{i}))\) be the estimated PIT CDF evaluated at \(\alpha\) and \(\rho(x,y)=(y-x)\mathbbm{1}(x<y)\). The regularization objective for level \(\alpha\) is defined as follows:

\[\mathcal{R}_{\text{Trunc}}(\theta;\mathcal{D},\alpha) \tag{13}\] \[=\begin{cases}\frac{1}{N}\sum_{i=1}^{N}\rho(Q_{\theta}(\alpha \mid X_{i}),Y_{i})&\text{if }\hat{F}_{Z}(\alpha)<\alpha\\ \frac{1}{N}\sum_{i=1}^{N}\rho(Y_{i},Q_{\theta}(\alpha\mid X_{i}))&\text{ otherwise}\end{cases} \tag{14}\]

This regularization objective adjusts \(\hat{F}_{Z}(\alpha)\) to match \(\alpha\) by increasing it when \(\hat{F}_{Z}(\alpha)<\alpha\), and vice versa. The final regularization objective is computed by averaging \(\mathcal{R}_{\text{Trunc}}(\theta;\mathcal{D},\alpha)\) over multiple quantile levels \(\{\,\alpha_{j}\,\}_{j=1}^{M}\):

\[\mathcal{R}_{\text{Trunc}}(\theta;\mathcal{D})=\frac{1}{M}\sum_{j=1}^{M} \mathcal{R}_{\text{trunc}}(\theta;\mathcal{D},\alpha_{j}). \tag{15}\]

It is worth noting that Chung et al. (2021) combine the previous regularization objective with a sharpness objective that penalizes the path between the quantile predictions, given by \(\frac{1}{M}\sum_{j=1}^{M}\frac{1}{N}\sum_{i=1}^{N}|Q_{\theta}(\alpha_{j}\mid X _{i})-Q_{\theta}(1-\alpha_{j}\mid X_{i})|\). Instead, we combine it with a strictly proper scoring rule.

### New Regularization-based Calibration Methods

Building upon the quantile calibration method discussed in Section 5.3.1, we propose two new regularization objectives which compute a differentiable PCE\({}_{p}\) using alternative statistical distances.

The first approach, named PCE-KDE, leverages the differentiable calibration map \(\phi^{\text{KDE}}\) (7) based on KDE. Given a set of quantile levels \(\{\,\alpha_{j}\,\}_{j=1}^{M}\), the regularization objective is given by

\[\mathcal{R}_{\text{PCE-KDE}}(\theta;\mathcal{D})=\frac{1}{M}\sum_{j=1}^{M} \left|\alpha_{j}-\phi^{\text{KDE}}(\alpha_{j};\{\,Z_{i}\,\}_{i=1}^{N})\right| ^{p}, \tag{16}\]

where \(p>0\). Note that \(\mathcal{R}_{\text{PCE-KDE}}\) reduces to PCE\({}_{p}\) in (2) when \(\tau\) in (7) goes to \(\infty\).

The second approach considers distances of the form \(\int_{0}^{1}|Q_{Z}(\alpha)-Q_{U}(\alpha)|^{p}d\alpha\), where \(Q_{Z}\) and \(Q_{U}\) denote the quantile functions of the true and uniform distributions, respectively. When \(p=1\), this distance reduces to the 1-Wasserstein distance, equivalent to \(\int_{0}^{1}|F_{Z}(\alpha)-F_{U}(\alpha)|d\alpha\), which aligns with PCE (see Proposition 1 in Appendix A.1).

By exploiting the fact that \(\mathbb{E}\left[F_{Z}(Z_{(i)})\right]=\nicefrac{{i}}{{N+1}}\), we approximate \(Q_{Z}(\nicefrac{{i}}{{N+1}})\) using the \(i\)-th order statistic \(Z_{(i)}\). The regularization objective is given by

\[\mathcal{R}_{\text{PCE-Sort}}(\theta;\mathcal{D})=\frac{1}{N}\sum_{i=1}^{N} \left|Z_{(i)}-\frac{i}{N+1}\right|^{p}, \tag{17}\]

where \(p>0\). Differentiable relaxations to sorting, such as those proposed by Blondel et al. (2020) and Cuturi et al. (2019), can be employed to obtain the order statistics.

## 6 A Comparative Study of Probabilistic Calibration Methods

In continuation of the empirical study described in Section 4, we now proceed to evaluate the performance of the probabilistic calibration methods outlined in the previous section.
Specifically, we apply eight distinct calibration methods to the three neural regression models introduced in Section 4. These methods are evenly divided into two categories: post-hoc methods and regularization-based methods.

To assess the effectiveness of these calibration methods, we employ four different evaluation metrics. The evaluation is conducted on a set of 57 datasets, utilizing the same experimental setup detailed in Section 4. To ensure a fair and consistent comparison, all the methods have been implemented within a unified codebase3.

Footnote 3: The code can be accessed at the following GitHub repository: [https://github.com/Vekteur/probabilistic-calibration-study](https://github.com/Vekteur/probabilistic-calibration-study)

### Experimental Setup

Base probabilistic models and calibration methods.We consider the three probabilistic models presented in Section 4, namely MIX-NLL, MIX-CRPS, and SQR-CRPS. For the MIX models, when applying quantile recalibration, we transform the CDF using the empirical CDF estimator (Rec-EMP), the linear estimator (Rec-LIN), or the KDE estimator (Rec-KDE). For SQR-CRPS, we transform multiple quantiles using conformalized quantile regression (CQR). For the three models, we consider the four regularization objectives presented in Sections 5.3 and 5.4 (with \(p=1\)), namely \(\mathcal{R}_{\text{PCE-KDE}}\) (PCE-KDE), \(\mathcal{R}_{\text{PCE-Sort}}\) (PCE-Sort), \(\mathcal{R}_{\text{QR}}\) (QR), and \(\mathcal{R}_{\text{Trune}}\) (Trunc). PCE-Sort is only shown in the Appendix because it performs similarly to PCE-KDE.

Metrics.We measure the accuracy of the probabilistic predictions using NLL and CRPS. For the SQR model, we estimate CRPS by averaging the quantile score at 64 equidistant quantile levels. Probabilistic calibration is measured using PCE, defined in (1). Finally, we measure sharpness using the mean standard deviation of the predictive distributions, denoted by STD.

Hyperparameters.In our experiments, MIX-NLL and MIX-CRPS output a mixture of 3 Gaussians, and SQR-CRPS outputs 64 quantiles. We justify the choice of these hyperparameters in Appendix C. The hyperparameter \(\tau\) of Rec-KDE and PCE-KDE is fixed at 100, which was found to perform well empirically. For regularization methods, an important hyperparameter is the regularization factor \(\lambda\). As previously observed in classification (Karandikar et al., 2021), we found that higher values of \(\lambda\) tend to improve calibration but worsen NLL, CRPS, and STD. Karandikar et al. (2021) proposed to limit the loss in accuracy by a maximum of 1%. We adopt a similar strategy by selecting \(\lambda\) which minimizes PCE with a maximum increase in CRPS of 10% in the validation set. For each dataset, we select \(\lambda\) in the set \(\{0,0.01,0.05,0.2,1,5\}\), which corresponds to various degrees of calibration regularization.

Comparison of multiple models over many datasets.As in Karandikar et al. (2021), and since NLL, CRPS and STD have different scales across datasets, we report Cohen's d, which is a standardized effect size comparing the mean of one method (over 5 runs, in our case) against a baseline. Values of \(-0.8\) and \(-2\) are considered large and huge, respectively. Due to the heterogeneity of the datasets that we consider, the performance of our models can vary widely across datasets. To visualize the results, we show the distribution of Cohen's d using letter-value plots (Hofmann et al., 2011), which indicate the quantiles at levels \(\nicefrac{{1}}{{8}}\), \(\nicefrac{{1}}{{4}}\), \(\nicefrac{{1}}{{2}}\), \(\nicefrac{{3}}{{4}}\) and \(\nicefrac{{7}}{{8}}\), as well as outliers. A median value below zero indicates that the model improved the metric on more than half the datasets.

In order to assess whether significant differences exist between different methods, we follow the recommendations of Ismail Fawaz et al. (2019), which are based on (Demsar, 2006). First, we test for a significant difference among model performances using the Friedman test (Friedman, 1940). Then, we use the pairwise post-hoc analysis recommended by Benavoli et al. (2016) using a Wilcoxon signed-rank test (Wilcoxon, 1945) with Holm's alpha correction (Holm, 1979). The results of this procedure are shown using a critical difference diagram (Demsar, 2006). The lower the rank (further to the right), the better performance of a model. A thick horizontal line shows a group of models whose performance is not significantly different, with a significance level of 0.05.

### Results

Figure 3 shows the letter-values plots for the Cohen's d of PCE (top panel) as well as the associated critical diagram (bottom panel), for all methods and datasets. The reference model is MIX-NLL. The results with other models as reference are available in Appendix B.1. Blue, green, and red colors are used for the post-hoc methods, the regularization-based methods, and the base models, respectively. The same information is given in Figures 4 and 5 for the CRPS and the NLL, respectively.

Comparison of PCE.As expected, Figure 3 shows that the PCE of calibration methods is improved compared to the base models. Furthermore, independently of the base model, we can see that post-hoc methods achieve significantly better PCE than regularization methods. When comparing PCE-KDE with QR, we can see that there is a significantly larger decrease in PCE with the MIX-CRPS base model compared to MIX-NLL. Finally, both PCE-KDE and Trunc decrease PCE for SQR-CRPS, without a significant difference between them.
Comparison of CRPS.While post-hoc methods outperform regularization methods in terms of PCE, Figure 4 shows they have a higher CRPS (except for the SQR base model). This can be explained by the fact that regularization methods prevent the CRPS from increasing exceedingly due to the selection criterion for \(\lambda\).

Comparison of NLL.Figure 5 shows the importance of the calibration map. In fact, quantile recalibration with a linear map significantly increases the NLL, while smooth interpolation decreases PCE without a large increase in NLL. Note that we only consider MIX models since we cannot compute the NLL for SQR.

On the choice of a calibration method.If probabilistic calibration is critical to the application, our experiments suggest that post-hoc methods such as quantile recalibration and conformal prediction should be preferred. However, when we also want to control the CRPS or the NLL, regularization methods can offer a better trade-off in terms of calibration and sharpness. In fact, as shown in Figure 6 in Appendix B.1, when the base model is MIX-NLL, all regularization methods provide a significant improvement in probabilistic calibration without deteriorating the CRPS, NLL or STD. For the MIX-CRPS model, Figure 7 shows that QR has limited impact on CRPS and NLL, while providing better calibration. For the SQR-CRPS base model, Figure 8 shows that the SQR-CRPS + CQR conformal method significantly outperforms the Trunc and PCE-KDE regularization methods both in terms of PCE and CRPS. Overall, Appendix B.1 suggests that MIX-NLL + PCE-KDE, MIX-CRPS + QR and SQR-CRPS + CQR are good choices for practitioners aiming to improve PCE without significantly impacting other aspects of the conditional distribution. Finally, since both regularization and post-hoc methods are able to improve calibration, we investigate whether a combination of these two methods can lead to better performance. Figure 9 in Appendix B.2 shows that such a combination does not significantly improve probabilistic calibration, with an increase in CRPS and NLL. This indicates that practitioners should exercise caution when applying regularization to a model that is already well-calibrated.

### Link between Quantile Recalibration and Conformal Prediction

Conformal prediction methods are well-known for their finite-sample coverage guarantee. Interestingly, a specific implementation of quantile recalibration can be considered a special case of conformal prediction. This implies that quantile recalibration can also provide a finite-sample coverage guarantee. This observation could potentially explain why both methods, conformal prediction and quantile recalibration, are effective in improving probabilistic calibration.

**Theorem 1**.: _Quantile recalibration is equivalent to Distributional Conformal Prediction (DCP) of left intervals at each coverage level \(\alpha\in[0,1]\). The equivalence is obtained

Figure 4: Comparison of CRPS with multiple base losses and calibration methods.

Figure 3: Comparison of PCE with multiple base losses and calibration methods.
_when the estimator of the calibration map is defined by a slightly different estimator than the conventional one in (6), namely \(\phi_{\text{DCP}}(\alpha)=\frac{1}{N^{\prime}+1}\sum_{i=1}^{N^{\prime}}\mathbbm{1} \left(Z_{i}^{\prime}\leq\alpha\right)\)._

Proof.: Given a predictive distribution \(F_{\theta}\) learned from a training dataset \(\mathcal{D}=\{\left(X_{i},Y_{i}\right)\}_{i=1}^{N}\) where \(\left(X_{i},Y_{i}\right)\stackrel{{\text{i.i.d.}}}{{\sim}}P_{X,Y}\), let \(Z_{i}^{\prime}=F_{\theta}(Y_{i}^{\prime}\mid X_{i}^{\prime})\) represent the PIT values computed on a separate calibration dataset \(\mathcal{D}^{\prime}=\{\left(X_{i}^{\prime},Y_{i}^{\prime}\right)\}_{i=1}^{N^ {\prime}}\), where \(\left(X_{i}^{\prime},Y_{i}^{\prime}\right)\stackrel{{\text{i.i.d.}}}{{ \sim}}P_{X,Y}\).

In the DCP approach, as outlined in Algorithm 2, the conformal scores are given by the PIT values \(Z_{i}^{\prime}\). DCP first computes the \(\alpha\) empirical quantile of the scores as \(\hat{q}=Z_{(\lceil(N^{\prime}+1)\alpha\rceil)}^{\prime}\), where \(Z_{(k)}^{\prime}\) represents the \(k\)th smallest value among \(\{\left.Z_{1}^{\prime},\ldots,Z_{N^{\prime}}^{\prime},+\infty\right.\}\). Then, the conformalized quantile is computed as \(Q_{\theta}^{\prime}(\alpha\mid X)=Q_{\theta}(\hat{q}\mid X)\), which corresponds to conformal prediction with coverage \(\alpha\) for the left interval \((-\infty,Q_{\theta}^{\prime}(\alpha\mid X)]\).

Let us consider quantile recalibration with the calibration map \(\phi\) in Algorithm 1 given by \(\phi_{\text{DCP}}(\alpha)=\frac{1}{N^{\prime}+1}\sum_{i=1}^{N^{\prime}} \mathbbm{1}(Z_{i}^{\prime}\leq\alpha)\). It computes a recalibrated CDF \(F_{\theta}^{\prime}\) by composing the original CDF \(F_{\theta}\) with \(\phi_{\text{DCP}}\), yielding \(F_{\theta}^{\prime}(y\mid X)=\phi_{\text{DCP}}(F_{\theta}(y\mid X))\).

We observe that \(\phi_{\text{DCP}}\) is the CDF of a discrete random variable, with \(\phi_{\text{DCP}}^{-1}(\alpha)=Z_{(\lceil(N^{\prime}+1)\alpha\rceil)}^{\prime}\) representing its empirical quantile function. Furthermore, the composition \(\phi_{\text{DCP}}\circ F_{\theta}(\cdot\mid X)\) acts as the inverse function of \(Q_{\theta}(\cdot\mid X)\circ\phi_{\text{DCP}}^{-1}\). As a result, both the DCP approach and quantile recalibration yield QFs and CDFs that correspond to the same underlying distribution.

Quantile recalibration with other recalibration maps (e.g., \(\phi^{\text{EMP}}\), \(\phi^{\text{LIN}}\), or \(\phi^{\text{KDE}}\)) would correspond to DCP where the empirical quantile \(\hat{q}\) is selected using other strategies which does not provide the exact conformal guarantee (8).

## 7 Conclusion

The observation that neural network classifiers tend to be miscalibrated (Guo et al., 2017) has prompted the development of various approaches for calibrating these models. In this paper, we present the largest empirical study conducted to date on the probabilistic calibration of neural regression models. Our study provides valuable insights into their performance and the selection of calibration methods. Notably, we introduce a novel differentiable calibration map based on kernel density estimation for quantile recalibration, as well as two novel regularization objectives derived from the PCE.

Our study reveals that regularization methods can provide a favorable tradeoff between calibration and sharpness. However, post-hoc methods demonstrate superior performance in terms of PCE. We attribute this finding to the finite-sample coverage guarantee offered by conformal prediction and demonstrate that quantile recalibration can be viewed as a specific case of conformal prediction.

Future investigations may extend the study of probabilistic calibration to other models, such as tree-based models, and explore alternative notions of calibration (Gneiting and Resin, 2021). Notably, distribution calibration represents a promising direction, as it has inspired the development of calibration methods (Song et al., 2019; Kuleshov and Deshpande, 2022).

## Acknowledgement

This work was supported by the Fonds de la Recherche Scientifique - FNRS under Grants T.0011.21 and J.0011.20.
