Distribution-Free Model-Agnostic Regression Calibration via Nonparametric Methods
Authors: Shang Liu
Zhongze Cai
Xiaocheng Li
Venue: NeurIPS 2023 poster (2023)
Keywords: Regression calibration
model recalibration
conditional quantile
nonparametric method

Abstract: In this paper, we consider the uncertainty quantification problem for regression models. Specifically, we consider an individual calibration objective for characterizing the quantiles of the prediction model. While such an objective is well-motivated from downstream tasks such as newsvendor cost, the existing methods have been largely heuristic and lack of statistical guarantee in terms of individual calibration. We show via simple examples that the existing methods focusing on population-level calibration guarantees such as average calibration or sharpness can lead to harmful and unexpected results. We propose simple nonparametric calibration methods that are agnostic of the underlying prediction model and enjoy both computational efficiency and statistical consistency. Our approach enables a better understanding of the possibility of individual calibration, and we establish matching upper and lower bounds for the calibration error of our proposed methods. Technically, our analysis combines the nonparametric analysis with a covering number argument for parametric analysis, which advances the existing theoretical analyses in the literature of nonparametric density estimation and quantile bandit problems. Importantly, the nonparametric perspective sheds new theoretical insights into regression calibration in terms of the curse of dimensionality and reconciles the existing results on the impossibility of individual calibration. To our knowledge, we make the first effort to reach both individual calibration and finite-sample guarantee with minimal assumptions in terms of conformal prediction. Numerical experiments show the advantage of such a simple approach under various metrics, and also under covariates shift. We hope our work provides a simple benchmark and a starting point of theoretical ground for future research on regression calibration.

Distribution-Free Model-Agnostic Regression Calibration via Nonparametric Methods
Shang Liu"
Zhongze Cai"
Imperial College Business School Imperial College London 1.lli211el w.llu21@@mperiin mete1emmmeera....mmmm ac.
Imperial College Business School Imperial College London cai 22@imperiall...m.mmmmmaem..miiii ac uk
Xiaocheng Li
Imperial College Business School
Imperial College London
xiaocheng.. xiaocheng.11@imppee wwwiimemriaall...mmmmmmmm. ac uk
Abstract
In this paper. we consider the uncertainty quantification problem for regression models. Specifically, we consider an individual calibration objective for characteriz- ing the quantiles of the prediction model. While such an objective is well-motivated from downstream tasks such as newsvendor cost, the existing methods have been largely heuristic and lack of statistical guarantee in terms of individual calibration. We show via simple examples that the existing methods focusing on population- level calibration guarantees such as average calibration or sharpness can lead to harmful and unexpected results. We propose simple nonparametrie calibration methods that are agnostic of the underlying prediction model and enjoy both com- putational efficiency and statistical consistency. Our approach enables better understanding of the possibility of individual calibration, and establish match- ing upper and lower bounds for the calibration error our proposed methods. Technically, our analysis combines the nonparametric analysis with covering number argument for parametric analysis, which advances the existing theoretical analyses in the literature of nonparametric density estimation and quantile ban- dit problems. Importantly, the nonparametric perspective sheds new theoretical insights into regression calibration in terms of the curse of dimensionality and reconciles the existing results on the impossibility of individual calibration. To knowledge, make the first effort reach both individual calibration and finite-sample guarantee with minimal assumptions in terms of conformal predic- tion. Numerical experiments show the advantage of such simple approach under various metrics, and also under covariates shift. We hope our work provides a simple benchmark and starting point of theoretical ground for future research on regression calibration.
Introduction
Modern machine learning methods have witnessed great success on a wide range of tasks in the past decade for their high accuracy in dealing with various kinds of complicated data. Uncertainty quantification has played an important role in the interpretation and risk assessment of machine learning models for downstream tasks such as optimization and decision-male makiin.... People observe
"Equal contribution.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

that the output of deep learning models tends to be over over-confident [Guo et al. 2017 Amodei aall||2011 |016] which inspires many recent works on uncertainty calibration for classification problems [Kumar et al.) 2019 Wenger et al. 2020 Luo et al. 2022]. Comparatively, there has been less systematic theoretical understanding of the uncertainty calibration for regression problems. For a regression problem, the output of a prediction model can be the prediction of the conditional mean (by minimizing mean squared error) or of the conditional median (by minimizing mean absolute error). But such a single-point prediction cannot characterize the uncertainty. which motivates the development of uncertainty calibration methods for regression problems. Some efforts have been made in order to estimate the entire distribution, which includes Bayesian ways (MacKay] 1992 Damianou and Lawrence 2013]. the frequentists' ways (Kendall and Gal 2017 Lakshminaaayaan
et 2017| Cui al., 2020] et al. 2020]] and the nonparametric ways Lei and Wasserman 2014 Lei et al. 2018 Bilodeau et al., 2021 Song et al.] 2019]. Another easier task is to predict
the quantiles rather than the whole distribution. Existing ways mainly focus completing the task
of quantile prediction in one go [Pearce et al., 2018 Thiagarajan et al. 2020 Chung et al.] [202]]
Takeuchi al. 2006 Stainwarr and Christmann 2011] However, all those methods suffer from either statistical inconsistency under model misspecification or computational intractability during
the training phase, or sometimes even both (see the detailed review in Appendix A}.
In this paper, we suggest that previous ways of training a new quantile prediction model from scratch while discarding the pre- pre-trained (mean) regression model may not be the best way for a quantile calibration objective, because the pre-trained model (though designed for a mean estimation objective) can be helpful for the quantile calibration. We propose simple, natural, but theoretically non-trivial method that divides the whole quantile calibration task into two steps: (i) train a good regression model and (ii) estimate the conditional quantiles of its prediction error. Although similar idea is
applied in split conformal prediction [Papadopoulos et al.) [2002] Vovk et all] [2005] [ei et al. 2018). the theoretical justification is still in lack since the conformal prediction only requires an average calibration guarantee (see Definition[1 T that can also be achieved without this first training step at all (where the detailed review given in Appendix A By careful analysis of the individual calibration objective (see Definition 3), we capture the intuition behind the two-step procedure and formalize it in mathematical guarantee. After a comprehensive numerical experiment on real-world datasets against existing benchmark algorithms, we suggest that one neither needs to estimate the whole distribution nor to train new quantile model from scratch, while a pre-trained regression model and a split-and- gellt-and-sutene aub-subttact suffice, both theoretically and empirically. Our contribution can be summarized below:
First, we propose simple algorithm that can estimate all percent conditional quantiles simultaneously. We provide the individual consistency of our algorithm and prove the minimax optimal convergence rate with respect to the mean squared error (Theorem T and 2 Our analysis is new and it largely relaxes the assumptions in the existing literature on kernel density estimation and order statistics. By showing the necessity of the Lipschitz assumption (Theorem 6), our result uses minimal assumptions to reach both finite sample guarantee and individual calibration, and our paper is the first to keep the latter two goals simultancously up to our knowledge.
Second, we propose a two-step procedure of estimating "mean quantile of error" rather than directly estimating the conditional quantiles, which enables a faster convergence rate both theoretically and empirically. Specifically, our convergence rate is of order Ö(L 442 dia ), where L is the Lipschitz constant of the conditional quantile with respect to the features, n is the number of samples, and d is
the dimension of feature. Since the conditional mean function and the conditional quantile function are highly correlated, one can greatly reduce the Lipschitz constant by subtracting the mean from the quantile.
Moreover, we construct several simple examples to show the unexpected behavior of the existing calibration methods, suggesting that a population-level criterion such as sharpness or MMD can be misleading. As our analysis works as a positive result on individual calibration, we also provide a detailed discussion about the existing results on the impossibility of individual calibration, illustrating that their definitions of individual calibration are either impractical or too conservative.
2

---

2 Problem Setup
Consider the regression calibration problem for given pre-trained model f. We given a dataset Y=) X y independent of the original data that trains f()). Here X' E [0, 1] TR denotes the covariate/fcature space and y C RR denotes the response/label space. The i.i.d samples {(X,, {. follow an unknown distribution P on X' x y. We aim characterize the uncertainty of the regression error U, = Y, f(X,) a model agnostic manner, i.e., without any further restriction on the choice of the underlying prediction model f.
Now we introduce several common notions of regression calibration considered in the literature. We first state all the definitions with respect to the error Ug and then establish equivalence with respect to the original Y, Define the 7th quantile of random variable Z Q- (Z) = inf F(t) 2 r for any T E [0, 1] where F is the c.d.f. of Z. Accordingly, the quantile the conditional distribution U z is denoted by Q-
z). A quantile prediction/caliiration model is denoted by (x) X + which gives the 7-quantile prediction of U given X
Definition (Marginal/Average Calibration). The model Q (x) is marginally calibrated if
P(U S (x()) =
for any r E [0, 1]. Here the probability distribution is with respect to the joint distribution of (U, X).
As noted by existing works, the marginal calibration requirement is too weak to achieve the goal of uncertainty quantification. A predictor that always predicts the marginal quantile of U for any X = z, i... (z) = (U)). is always marginally calibrated, but such a model does not capture the heteroscedasticity of the output variable Y the error U with respect X.
Definition 2 (Group Calibration). For some pre-specified partition X' Xi U...U Xk. the model (x) is group calibrated if
P(U s Qx())X E =T
for any T E [0, 1] and & L,..., K.
Group calibration is stronger notion of calibration than marginal calibrations. often considered
in a related but different problem called conformal prediction [Vovk] 2012 Lei and Wasserman 2014 Alaa et al. 2023) where the goal is to give sharp covering set C(X) (or covering band if unimodality is assumed) such that IPYY Č(X) X E Xx) 2 1 Vk. Foygel Barber et al. 2021 consider the case where the guarantee is made for all IP(X xx) 2 for some d > O.
Definition 3 (Individual Calibration). The model (x(2) is individually calibrated if
P(U < (x)XX = =) =T
for any T E [0, 1] and x E X. Here the probability distribution is with respect to the conditional distribution U X.
In this paper, we consider the criteria of individual calibration, an even stronger notion of calibration than the previous two. It requires the calibration condition hold with respect to any x E X' in pointwise manner. The difference between the conformal prediction problem with group calibration and the regression calibration problem with individual calibration is that conformal prediction does not make any assumptions on the data distribution, while our lower bound result of Theorem 6 suggests that without any assumption, the convergence rate of the conditional quantile estimator can be arbitrarily slow. The conformal prediction literature considers weaker notion of calibration than the individual calibration setting considered in this paper where some mild assumptions are made.
Finally, the following proposition says that any quantile prediction model on the error U is equivalent corresponding quantile prediction on Y. The result is stated for the individual calibration, while a similar statement can be made for marginal calibration and group calibration as well.
Proposition For any predictor (x) on the th quantile of U|X I, we have
PP PY s f(x)+ = =)) = =P(U (x(x) = =x).
3

---

Oracle
Deep Ensemble
Our Algorithm NRC
106
150
79
25-
-25-
-50
Observations
Observations tions
Observations
Predittions
Predictions 951% Interval
Predictions
95,% Interval
Interval
20.0
12.0
120
Figure 1: Synthetic data where the underlying distribution is obtained by a combination of sine functions. The solid lines denote predicted means, the shaded area denotes predicted intervals between 97.5% and 2.5% quantiles, and the yellow dots denote a subset of real observations. The leftmost plot gives the real mean as well as oracle quantile values, while the rest two plots are predictions from different calibration models. The middle plot is produced by a Deep Ensemble [Lakshminarayan al. 2017 of HNNs trained with 40,000 samples, which is both common benchmark and building block for several regression calibration methods. The rightmost plot is produced by our proposed nonparametric calibration method Algorithm 2 of which the base regression model an ordinary feed-forward regression network. The detailed setup is given in Appendix B.2
2.1 Motivating examples
In addition to the above calibration objective, existing literature also considers sharpness (the average length of confidence intervals (Zhao al.) 2020, Chung et al. (02T]) and maximum mean discrepancy (MMD) principle [Cui et al.] 2020]. To the best of our knowledge, all the existing works consider population-leve objectives such as marginal calibration, sharpness, and MMD when training the calibration model, while we are the first work that directly aims for individual calibration. Here we use a simple example to illustrate that such population-level objectives may lead to undesirable calibration results, and in Appendix B we elaborate with more examples the inappropriateness of existing methods for regression calibration.
Example Consider X ~ (n 1] and the conditional distribution Y|X = x ~ Unif æ=] Then if one aims for 90%, the outcome from sharpness maximization subject to marginal calibration is Q (x) for [0, 0.9] and (x) 0 for E (0.9, 1]. Consequently, the quantile prediction has a 100% coverage over x E [0, 0.9] but 0% coverage over x E (0.9, 1].
Generally, the population-level criteria such as sharpness or MMD may serve as performance metrics to monitor the quality of calibration outcome. However, the inclusion of such criteria in the objective function will encourage the calibration result to be over-conservative for the low variance region (x: E [0, 0.9) while giving up the high variance (x: E (0.9, 1)), which highly undesirable for risk-sensitive applications and/or fairness considerations.
Besides, individual calibration is sometimes also necessitated by the downstream application. For example, in some decision-making problems, regression calibration serves as a downstream task where loss is measured by the pinball loss; say, the newsvendor problem Kullehhvv et al. 2018 Ban and Rudin 2019) in economics and operations research uses pinball loss to trade-off the backorder cost and the inventory cost.
Proposition 2. The pinball loss with respect to is defined by
(ui,, u2) = =(1 (1-r)(u - u2) l(u; > u2) + T(u2 wi) Iful < u2).
We have
Q2
arg arg mig (y,Y z|.
yER
Proposition serves the original justification of pinball loss for quantile regression. More importantly, says that calibration method that is inconsistent with respect to individual calibration can be substantially suboptimal for downstream tasks such as the newsvendor problem.

---

Main Results
3.1 Generic nonparametric quantile estimator
In this subsection, we present and analyze simple algorithm of a nonparametric quantile estimator, which can be of independent interest itself and will be used as a building block for the regression calibration. The algorithm takes the dataset {Xx,, ())=. input and outputs r-quantile estimation of the conditional distribution U U|X = x for any 7 and æ.
Algorithm Simple Nonparametric Quantile Estimator
Input: 7 (0, 1), dataset {(X,, U.) kernel choice (x(:.)) I: Estimate the distribution of U|X by (where 6, is point mass distribution at u)
Di=l (ar X,)ou,
PUjx=r
D Ky (2,
2: Output the conditional 7 quantile by minimizing the pinball loss on PU|X=r
SNQ = arg min E UN x-n Pr (u, U)].
(1)
The minimum of the optimization problem Algorithm i indeed achieved by the 7-quantile of the empirical distribution Pujx=z- In the following, we analyze the theoretical properties of the algorithm under the naive kernel
Ky (=1, x2) = 1[ll=1 - x2ll S h}.
where h is the hyper parameter for window size, and || - æ2l denotes Euclidean distance.
Assumption 1. We assume the following on the joint distribution (U, and T E (0, 1) of interest:
(a) (Lipschitzness). The conditional quantile function is L-Lipschitz with respect z,
IQ- =x1) - -Q-(UX = =x22] S L. L| - -xall, V01, E X.
(b) (Boundedness). The quantile Q+UU| x) is bounded within MM all x E X.
(c) (Density). There exists a density function (u) for the conditional probability distribution U|X I. and the density function is uniformly bounded away from in neighhorhood of the interested quantile. That is, there exist constants p and r.
(u) P.
Vx E X and |u Q- =x)] Sr.
In Assumption 1 part (a) and part (b) impose Lipschitzness and boundedness for the conditional quantile, respectively. Part (c) requires the existence density function and lower bound for density function around the quantile of interest. Its aim is to ensure locally strong convexity for the expected risk (4 U)). Under Assumption rate for Algorithm
we establish consistency and convergence
Theorem 1. Under Assumption 7 Algorithm D statistically consistent, i.e., for any >> 0,
lim P @SNe (X) Q-(UIX) 12+000
26)=0. =0.
26)
Furthermore, by choosing
O(L uंz atz), ), we have for sufficiently large 2 C, when L 0,
E (x)) Q-(U1x) < a' Joं.
where and depends polynomially on / y. M, and log(n). In addition, whens = 0. the hight-hand-side becomes Ô(L).
5

---

While Algorithm I seems to be a natural algorithm for nonparametric quantile regression, Theorem provides the first convergence rate for such a nonparametric quantile estimator, to the best of our knowledge, in the literature of nonparametric regression; also, it is the first theoretical guarantee towards individual calibration for regression calibration. The standard analysis of nonparametric mean estimation [Györfi et al. 2002) cannot be applied here directly because the algorithm involves a local optimization QI). Even more challenging, the optimization problem (T) aims to optimize the quantile of the conditional distribution U X = but the samples used in (1 are from distributions
U|X = Zi. which causes a non-ii.d. issue. To resolve these issues, we combine the idea of bias and variance decomposition in nonparametric analysis with covering concentration argument for the local optimization problem. The detailed proof is deferred to Appendix H
Theorem I also complements the existing results on finite- sample analysis for quantile estimators. One line of literature [Szorenyi al. 2015 Altschuler et al. 2019 establishes the quantile convergence from the convergence of empirical processes, and this requires additional assumptions on the density function and does not permit the non-i.i.d. structure here. The quantile bandits problem also entails the convergence analysis of quantile estimators; for example, Zhang and Ong (2021) utilize the analysis of order statistics [Boucheron and Thomas 2012], and the analysis inevitably requires a non-decreasing hazard rate for the underlying distribution. Other works [Bilodeau et al. [02T] that follow the kernel density estimation approach require even stronger conditions such as realizability. We refer to Appendix for more detailed discussions. In comparison to these existing analyses, analysis is new, and it only requires a minimal set of assumptions. The necessity of the assumptions can be justified from the following lower bound.
Theorem 2 rephrases the lower bound result in the nonparametric (mean) regression literature [Györti et al. 2002 for the quantile estimation problem. It states that the convergence rate in Theorem is minimax optimal (up to poly-logarithmic terms). It is easy to verify that the class p satisfies establish the criticalness of Assumption Assumption [ In this light, Theorem I and Theorem 2
| Furthermore, in Theorem ij Appendix C. we establish that the convergence rate of any estimator can be arbitrarily slow without any Lipschitzness assumption the conditional statistics.
Theorem 2. Let PL be the class of distributions of (X, U) such that X ~ Uniform [0, 1]s. U- = μ(X)= + N. and p(x) is L-Lipschitz where N is an independent standard Gaussian random variable. For any algorithm, there exists a distribution in the class pt. such that the convergence rate of the algorithm is Q2L 1:2 = 2i2).
Our result here provides a guarantee for the estimation of Q, (UI and thus a positive result for individual calibration with respect to a specific т. We note that the result does not contradict the negative results on individual calibration (Zhao et al. 2020 Lei and Wasserman 2014]. Zhao et al.
(0200 measure the quality of individual calibration via the closeness of the distribution Fx (Y) to the uniform distribution. In fact, such a measurement is only meaningful for a continuous distribution of YJX, but they prove the impossibility result based on discrete distribution of Y So, their negative result on individual calibration does not exclude the possibility of individual calibration for a continuous Y|X. Lei and Wasserman (2014) require an almost surely guarantee based on finite observations and prove an impossibility result. This is a very strong criterion; our results in Theorem [ are established for two weaker but more practical settings: either the strong consistency in the asymptotic case as n > 3o or a mean squared error guarantee under fimite-sample. We defer more discussions to Appendix C
3.2 Regression calibration
Now we return to regression calibration and build a calibrated model from scratch. Specifically, Algorithm 2 spllits the data into two parts; it uses the first part to train a prediction model and the second part to calibrate the model with Algorithm [ The first part can be skipped if there is already a well-trained model f where it becomes a recalibration problem.
The theoretical property of Algorithm 2 can be established by combining the results from Section 3.1 with Proposition[] In Algorithm 2 the quantile calibration allows full flexibility in choosing the regression model f and does not require an associated uncertainty model to proceed. The motivation for calibrating the error U; instead of the original Y; is two-fold: First, if one treats Y, itself as U; and applies Algorithm|] []] then essentially restricts the prediction model to be a nonparametric one. Second, the reduction from Y, to Ug gives better smoothness of the conditional distribution (a

---

Algorithm 2 Nonparametric Regression Calibration
Input: Dataset {(X,, r. kernel choice Kn (;.
1: Split the dataset into half: D, = =(X;, mi D2 = ((X,, Yi))
2: Use Di to train regression model f.
3: Calculate the estimation error of f on D2: U4 Y2 f(Xi), i= 1, n.
4: Run Algorithm 1 on the data {(X,, m-))= and obtain Q7
5: Return f(.) + (.).
smaller Lipschitz constant in Assumption (a). And thus it will give a faster convergence rate. We remark that Algorithm [ induces independence between the training of the prediction model and that the calibration model. The design is intuitive and important because otherwise, it may result in predicting smaller confidence intervals both theoretically and empirically.
3.3 Implications from the nonparametric perspective
The nonparametric approach gives positive result in terms of the possibility of individual calibration, but it pays a price with respect to the dimensionality d. The following theorem states that such a dimensionality problem can be addressed under conditional independence. And more generally, it provides a guideline on which variables one should use to perform calibration tasks.
Y -(()). function, and Z C X is do-dimensional subspace of X. If X and U are mutually independent conditioned on Z (i.e, X 7 - U a Markov chain), then is lossless to perform calibration using only (Ut, Z.Y's. s. In particular, applying Algorithm
Theorem 3. Suppose U = =
Suppose = =m(X) E Z, where m some measurable
Ianiiimoonn oo (U,, Z,) will yield the same consistency but with a faster convergence rate of ( 1613).
Theoretically, one can identify such by independence test for each component of and Y f(X) on the validation set. Some other practical methods such as random projection and correlation screening selection are shown in Section When the conditional independence Theorem 3 does not hold, it is still not the end of the world if one calibrates (U;, Z,). The following theorem tells that conditional calibration with respect to Z, as an intermediate between marginal calibration and individual calibration, can be achieved.
Theorem 4. Suppose U Y f(x) and a sub-feature of X. Without loss of generality, we assume that Z = Ndd (X), where Na, projects from X onto do-dimensional subspace 2. Then
P(U S (Q-(UZZ)
Z) =,
which implies that IPY S-YY$/F +
2) =
A notable implication from the above theorem is that if we want to calibrate the model against certain sub-features such as age, gender, geo-location, etc., we can summarize these features in and use Algorithm | with respect to (U:, Z.).
However, if return to the original pinball loss, the following theorem inherits from Proposition 2 and states that the inclusion of more variables will always be helpful in terms of improving the pinball loss. highlights a trade-off between the representation ability of the quantile calibration model (inclusion of more variables)and its generalization ability (slowing down the convergence rate with higher dimensionality).
Theorem 5. Suppose U =Y f(X). Suppose ((d1) and ((d2) are the first di and da components of X. Suppose Q+ (=(Zzia)) is the true th quantile of U conditioned Qr (Fmargiu) is the true 7tth marginal quantile of U. If0 d1 < d2 < d, then
%(d.), i 1.2. Suppose
Y)) 2 (1 Y) (LX(())) x))
7

---

4 Numerical Experiments
In this section, we evaluate our methods against a series of benchmarks. We first introduce the evaluation metries that are widely considered in the literature of uncertainty calibration and conformal prediction, including Mean Absolute Calibration Error (MACE) which calculates the average absolute error for quantile predictions from 0.01 to 0.99; Adversarial Group Calibration Error (AGCE) which finds the sub-group of the test data with the largest MACE; Check Score which shows the empirical pinball loss; Length which measures the average length for the constructed 0.05-0.95 intervals; and Coverage which reflects the empirical coverage rate for those 0.05-0.95 intervals. Formal definitions of these calibration measurements are provided in Appendix F.T
The benchmark methods are listed as follows (the details can be found in Appendix E.4 For the Gaussian- based methods, we implement the vanilla Heteroscedastic Neural Network (HNN)
[endall and Gal 2017]. Deep Ensembles (DeepEnsemble) [Lakshminarayanan et al.) 2017]. and Maximum Mean Discrepancy method (MMD) [Cui et al., 2020J. We also implement MC MC-Dropout (MCDrop) Gal and Ghahramani 2016) and Deep Gaussian Process model (DGP) [Damianou and Lawrence 2013 For models that do not rely on Gaussian assumption, we implement the combined calibration loss (CCL) introduced in Sec. 3.2 of [Chung et al.] (021]. We also implement post-hoc calibration methods, such as isotonic regression (ISR) suggested by [Kuleshov et al. 2018] with a pre-trained HNN model. Another post-hoe method named Model Agnostic Quantile Regression (MAQR) method in Sec. 3.1 Chung al. (202T) is also implemented. For conformal prediction algorithms, we implement Conformalized Quantile Regression (CQR) [Romano et al. 2019 and Orthogonal Quantile Regression (OQR) [Feldman et al. 2021).
As our proposed nonparametric regression calibration (NRC) algorithm is agnostic with respect to the underlying prediction model, we implement feed-forward neural network and random forest model as the prediction model for NRC. Also, we apply the dimension reduction idea (called NRC-DR) including random projection and covariate selection where all technical details given in Appendix E.3 and E.5 All codes are available at https: //tttugb com /ww/////nonn///iiiii /NRC
In the following subsections, we summarize our results on UCI datasets, time series, and high- dimensional datasets. Additional experiments on covariate shift are deferred E.7
4.1 UCI dataset experiments
We evaluate our NRC algorithms the standard 8 UCI datasets Dua and Grafi (2017). Part of the experiment results on representative benchmark models is given in Table ! and the full results are presented in Appendix E.5 due to space limit. Additional details on dataset description, hyperparameter selection, and early stopping criteria are given in Appendix E..5
The result supports the competitiveness of our algorithm. In terms of MACE, our methods achieve the minimum loss on 6 out of 8 datasets, with the non-optimal results close to the best performance seen on other benchmarks. The rest metrics also witness our algorithms outperforming the benchmark models, with a superior result on datasets in terms of AGCE, and 4 in terms of CheckScore. Note that even for the rest 4 datasets where our NRC algorithm does not achieve the best CheckScore (the empirical pinball loss), our performance is still competitive compared to the best score, while for some specific datasets (such as Energy and Yacht) we achieve significant improvement over the CheckScore, which verifies the potential of our methods for the related downstream tasks. As for the Length metric (which evaluates the average interval length of a symmetric 90% confidence interval), it is not easy to trade-off between the sharpness requirement and the coverage rate due to the fact that one can trivially sacrifice one for the other, so we omit the direct comparison. However, we can still observe that our NRC algorithm behaves reasonably well, with a stable coverage rate performance and relatively sharp interval length amongst the rest benchmarks.
4.2 Bike-sharing time series dataset
This experiment explores the potential of our algorithm to quantify the uncertainty of the time series objective, where independent variables are aggregated a sliding window fashion from the past few days. The Bike-Sharing dataset provided in Faaaeeee- and Gama (0114) is used to visually evaluate our proposed algorithm. Due to the limitations in the representation ability of the feed-forward structure, we instead deploy LSTM networks for the time series data. We design LSTM-HNN,
8

---

Table 1: Experiments the standard UCI datasets. Each experiment (for combination of dataset and method) is repeated times order to estimate the range of value fluctuation of each metric. On each dataset, the minimum loss for each metric is marked in bold and red, while omit the comparisons on the Length and the Coverage due to the difficulty of trading off between sharper intervals and better coverage intervals. Full experiment details deferred to Appendix E
Daaand
Mere
Decpisscmble
ISR
XEAOR
CQR
NBC
NRC-ER
Bostnz
CheckSooee
AGCE
0.226:004
19852
027H0007
1s0.3
0.27-1.0
0:2500005
0.365600
8.73:6.1
0.19-0055
1.1:03
0.19:8045
LeighhCosenge
15::108))
9.211999)
70926)
1051/9499
9.7 ()
MACE
03442000
004P21000
0047:5002
AGEE:
02.0004
BJK ILBS
19-00
0 18-000
0.211005
0.15 0.05
8 19 19.
ChekSron
16:03
19:0X
255.0
2.102
172.86
1.4 8,2
2.18
LeephCCorrags
(S198955)
1999999)
MACE
116-002
01K1-004
0015-1000
00064-0000
B 0,948
Enargy
CheckScore
AGCE
021:301
14:33
0.10::0004
06428000
021:0.00
0.600.11
021:0044
0.23022
0.18:0000
0,10:0002
0.22:006
0 0.24
0.5220.2
LonathCoo
924100%)
3.95)))))
62630949)
S59.1/999)
MACE
01341001
0077
0013
001120000
00412002
0JI13 6,864
0 L00907
Kasšen
CheckSrore
AGCE
0065-9000
007564
0:51 -0.000
O4TZ-4000
0000-000001
0
0005-00000
0U775-0000
0037-00000
DC07-0 03
10002-00000
LeeghCoornpe
MACE
0.125000
0811:99005
0.43500000
C477=8077
0:012:0:005
0001700/077
Navi
ChuidSore
AGCE
00013
0.15=005
8.0022
00000701000007
8004200000
QU001 100007
0001:0000
CUP98:3005
0.05288860
0001
0
BU002-00006
MACH
626-001
BBRS 000
0079-1000
G 18 M1T
AGET:
81-0M
-000
CheckScore
1220003
1.2-0.05
20630979)
Iean 5
12:02208)
MACE
02263992
00002000000
0102928.00
0816::9
2017250006
AOCE
0094:10.00
#007
LegghCosrryy
4-40 (99055
5-40 12965)
12-00000053
MACE
0.11-003
0088-3000
CheckSeoee
19:36
1.1004
1.5000
0322-009
0.1920007
0.18-683
0,460.2
(8150495)
10
(0:0915844)
3.2 18336)
17240 42725)
S.3119956
122:1199))
which is ordinary LSTM network followed by fully connected layer with output units; and LSTM-NRC, whose regression model shares the same LSTM setting as LSTM-HNN except for the ending layer having only one output unit (without variance/uncertaintyy unit). For more details, Appendix E.6
LSTM-HNN
Observations
#12-
Predictions
951% Interval
330-
LSTM-NRC
1000
Predictions
95.% Interval
622-
428-
160
125 150 175 200
Time Stamp
Figure 2: Experiments on the Bike-Sharing dataset. The demonstrated time stamps are consecutively taken from the test set, thus temporally separated from the training set. The top figure gives the predicted 95% confidence intervals of an LSTM-HNN model, and the bottom figure is produced from random projection variant of our NRC method (LSTM-NRC.. The target variable shows strong periodicity, which both two models capture. Nevertheless, our NRC method delivers sharper confidence intervals and is better adapted the underlying distribution shift. In particular, our method also has a better coverage rate on the extreme values at time stamps between 100 and 200.

---

Figure 2 visually compares the two models on the 95% confidence interval prediction task on the test [E for how the test set is split out). A first observation is that our LSTM-NRC is producing sharper and better calibrated 95% intervals. Also, in spite of the satisfying coverage rate of both the methods on time range 0 to 100 (time stamp 11 is exactly the terminating point of the training set, thus time range 0 to 100 roughly denotes the near future of the historical events), on time range 100 to 200 the behaviors of the two models diverge. After an evolution of 100 time stamps the target variable (rental count per hour) is under shift in distribution: the peak values get higher than historical records would foresee. Confidence intervals produced by LSTM-HNN fail to cover these "extremer" peaks, while our LSTM-NRC method could easily adapt it, taking advantage of a high-quality regression model. Note that both LSTM-HNN and LSTM-NRC share the same regression architecture, argue that the additional variance unit of LSTM- MNNN may hurt the mean prediction performance of the model.
set (see Appendix
4.3 Higher-dimensiona. datasets
One major drawback of the non-parametric type methods is its unsatisfying dependence the dimension of the data. Our paper attacks such an obstacle from two aspects simultancously: on the one hand, subtract regression model from the original data to gain smaller Lipschitz constant I.; on the other hand, we apply several dimensional reduction methods to reduce the data dimension manually. In this paper, we mainly implement three of them: random projection, covariate selection, and second-to-last layer feature extraction, of which the details can be found in Appendix We call the variants of our NRC algorithm NRC-RP, NRC-Cov, and NRC-Embed, respectively. We examine them against several high-dimensional datasets: the medical expenditure panel survey (MEPS) datasets (panel 19 2017 panel 20 2017 panel 21 2017], suggested in [Romano et al.)
E.3
2019
The datasets' details are also deferred to Appendix E.3
Dataset
Metric
DeepEnsemble
ISR
CQR
OQR
NRC-RP
NRC-Embed
MACE
0.2611000
0.12+0.05
0.05=0.01
0.055-0.02
0.024=0.01
0.0087=0.006
AGCE
0.261:0.04
0.1411004
0.079:10.01
0.091-0003
0.062:0.02
0.04:10.02
19
CheckScore
25-77
2.4
2.54+0.2
3,24-0.4
2.4--0.7
Length(Coverago
230=140(99%)
50++6(93%)
28+:2899%)
25+:2(85%)
35++8(00%)
24=005(90%%
MACE
0.25++0.02
0.16+ 0.08
0.046+00.01
0.057-+0.024
0.0052-0.000
0.0058-0.005
AGCE
0.25++0.02
0.18+00.88
0.061=0.01
0,08++0.03
0.033=0.01
0.041=00.07
meps__2
CheckScore 190.55((99%)
23=88
7.4=1
2.520.11
2.510.1
3:20.4
2.4:10088
Lungthi(uveerage)
69.38(94%)
28.1((89%)
27.14(85%)
32.15(90%%
27.13(91%)
MACE
0.25:10.02
0.13:20.1
0.05
0.0087-0.004
0.012110.02
0.0088:0.005
AGCE
0.26:1.002
0.15
0.063.10.008
0.075-0.02
0.045
0.041
meps_21 21
CheckScore
36+-20
23--10
2.5+-0.3
2.6+00.3
3.5+-0.7
2.7+1
LengthiCovernge)
160::40(999)
43=20(93%)
27::3(86%)
29-3(
394+300%%
24#:(118)) %)
Table 2: Experiments on MEPS Datasets
In Table 2 we summarize the dimension reduction variants of our NRC algorithm against several benchmarks. The result sees a substantial advantage of our method over the selected benchmark methods, while most of the unshown benchmarks either are too computationally expensive to run or do not deliver competitive results on the high-dimensional dataset.
5 Conclusion
In this paper, we propose simple nonparametric methods for regression calibration, and through its lens, we aim to gain better theoretical understanding of the problem. While numerical experiments show the advantage, we do not argue for universal superiority of our methods. Rather, we believe the contribution of our work is first to give a positive result for individual calibration and an awareness of its importance, and second to show the promise of returning to simple statistical estimation approaches than designing further complicated loss functions without theoretical understanding. Future research includes (i) How to incorporate the uncertainty output from a prediction model with our nonparametric methods; (ii) While we do not observe the curse of dimensionality for our methods numerically and we explain by the reduction of Lipschitzness constant, it deserves more investigation on how to adapt the method for high nigh-dimensional data.
10

---

References
Ahmed M Alaa, Zeshan Hussain, and David Sontag. Conformalized unconditional quantile regression. In International Conference on Artificial Intelligence and Statistics, pages 10690-10702. PMLR, 2023.
Jason M Altschuler, Victor-Emmanuel Brunel, and Alan Malek. Best arm identification for contami- nated bandits. J. Mach. Learn. Res., 20(91):1-395 2019.
Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané. Concrete problems in ai safety. arXiv preprint 1eXiv:6606006655.. 2016.
Gah-Yi Ban and Cynthia Rudin. The big data newsvendor: Practical insights from machine learning.
Operations Research, 67(1)900-108, 2019.
Peter L Bartlett, Olivier Bousquet, and Shahar Mendelson. Local rademacher complexities. The 2005. 2005.
Annals of Statistics,
Blair Bilodeau, Dylan J Foster, and Daniel M Roy. Minimax rates for conditional density estimation via empirical entropy. arXiv preprint arXiv:2 109. 10461, 2021.
Stéphane Boucheron and Maud Thomas. Concentration inequalities for order statistics. Electronic Communications Probability, 17(none):1-12, 2012.
Asaf Cassel, Shie Mannor, and Assaf Zeevi. A general approach to multi-armed bandits under risk criteria. In Conference on learning theory, pages 1295-1306. PMLR, 2018
Xiangli Chen, Mathew Monfort, Anqi Liu, and Brian D Ziebart. Robust covariate shift regression, In Artificial Intelligence and Statistics, pages 1270-1279. PMLR, 2016.
Youngseog Chung, Willie Neiswanger, Ian Char, and Jeff Schneider. Beyond pinball loss: Quantile methods for calibrated uncertainty quantification. Advances in Neural Information Processing Systems, 34: 19971-10848,, 2021.
Peng Cui, Wenbo Hu, and Jun Zhu. Calibrated reliable regression using maximum mean discrepancy. Advances in Neural Information Processing Systems, 13:17164-77175, 2020.
Andreas Damianou and Neil D Lawrence. Deep gaussian processes. Artificial intelligence and statistics, pages 207-215. PMLR, 2013.
Sanjoy Dasgupta. Experiments with random projection. arXiv preprint 1rXiv:1301.3849,200
Dheeru Dua and Casey Graff. Uci machine learning repository, 2017. httppttt://arrciie... uci. edu/ml
Aryeh Dvoretzky, Jack Kiefer, and Jacob Wolfowitz. Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. The Annals of Mathematical Statistics, pages 642-669, 1956.
Hadi Fanaee-T and Joao Gama. Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, 2:1 23113,, 2014.
Shai Feldman, Stephen Bates, and Yaniv Romano. Improving conditional coverage via orthogonal quantile regression. Advances in neural information processing systems, 34:2060-2017 2021.
Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. The limits of distribution-free conditional predictive inference. Information and Inference: A Journal of the IMA, 10(2):455-422, 2021.
Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning, pages 1050-1059, PMLR, 2016.
Isaac Gibbs, John J Cherian, and Emmanuel Candès. Conformal prediction with conditional guarantees. arXiv preprint arXiv:2305.126111 2023.

---

Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger. On calibration of modern neural
networks. In International conference on machine learning, pages 1321-1330. PMLR, 2017.
László Györfi, Michael Köhler, Adam Krzyžak, and Harro Walk. distribution-free theory of nonparametric regression, volume 1. Springer, 2002.
Kurt Homik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are universal approximators. Neural networks, 2(5):359-366, 1989.
Steven R Howard and Aaditya Ramdas. Sequential estimation of quantiles with applications to a/b testing and best-arm identification. Bernoulli, (8(3):1704-177 2022.
Alex Kendall and Yarin Gal. What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems, 30, 2017.
Ravi Kumar Kolla, LA Prashanth, Sanjay P Bhat, and Krishna Jagannathan. Concentration bounds for empirical conditional value-at-risk: The unbounded case. Operations Research Letters, 47(1): 16-20, 2019.
Vladimir Koltchinskii. Local rademacher complexities and oracle inequalities in risk minimization. The Annals of Statisticcs, pages 2593-2656, 2006.
Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon. Accurate uncertainties for deep learning using calibrated regression. In International conference on machine learning, pages 2796-2804. PMLR, 2018.
Ananya Kumar, Percy S Liang, and Tengyu Ma. Verified uncertainty calibration. Advances in Neural Information Processing Systems, 32, 2019.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems, 30, 2017.
Jing Lei and Larry Wasserman. Distribution-free prediction bands for non- con-parametric regression. Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 71-96, 2014.
Jing Lei, Max G' Sell, Alessandro Rinaldo, Ryan Tibshirani, and Larry Wasserman. Distribution-fres predictive inference for regression. Journal of the American Statistical Association, 1 113(523): 1094-1111, 2018.
Rachel Luo, Aadyot Bhatnagar, Yu Bai, Shengjia Zhao, Huan Wang, Caiming Xiong, Silvio Savarese, Stefano Ermon, Edward Schmerling, and Marco Pavone. Local calibration: metrics and recalibra- tion. In Uncertainty Artificial Intelligence, pages 1286- 1295. PMLR, 2022.
David JC MacKay. A practical bayesian framework for backpropagation networks. Neural computa- tion, 4(3):448-47,, 1992.
Pascal Massart. The tight constant in the Probability, pages 1269-1283 1990.
inequality. The annals of Andreas Maurer and Massimiliano Pontil. Empirical bernstein bounds and sample variance penaliza- tion. arXiv preprint arXiv:0907.3744 2009.
Medical expenditure panel survey panel 19, 2017. URL laleees__eatlll.......
ahrq
Medical expenditure panel survey panel 20, 2017. URL http:::ttptssps..teeetea
tttaa_sttttt/oooooadddddo hidd_dat__filees ctalless_eeellllee1 httpl/aeddddaaallppppppoooooool hat.i///ppppooufffuuubrrr-------
Medical expenditure panel survey panel 21, 2017. URL https:://eps..attt http://pppp..arrr..oo///////s/// data_stata/doomlood fit.aa.fleee.ddaaa...ppppcoooooo
Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19-23, 2002 Proceedings 13, pages 345-356. Springer, 2002.
12

---

Tim Pearce, Alexandra Brintrup, Mohamed Zaki, and Andy Neely. High-quality prediction intervals for deep learning: A distribution- free, ensembled approach. In International conference on machine learning, pages 4075-4084. PMLR, 2018.
Joaquin Quinonero-
Masashi Sugiyama, Anton Schwaighofer, and Neil D Lawrence. Dataset shift in machine learning. Mit Press, 2008
Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019.
Hugh Salimbeni and Mare Deisenroth. Doubly stochastic variational inference for deep gaussian processes. Advances in neural information processing systems, 30, 2017.
Hao Song, Tom Diethe, Meelis Kull, and Peter Flach. Distribution calibration for regression. In International Conference on Machine Learning, pages 5897-5906. PMLR, 2019.
Ingo Stainwart and Andreas Christmann. Estimating conditional quantiles with the help of the pinball loss. Bernoulli, 17(1):211-223 2011.
Balazs Szorenyi, Róbert Busa-Fekete, Paul Weng, and Eyke Hüllermeier. Qualitative multi-armer bandits: A quantile-base -ased approach. In International Conference on Machine Learning, pages 1660-1668. PMLR, 2015.
Ichiro Takeuchi, Quoe V Le, Timothy D Sears, and Alexander J Smola. Nonparametric quantile estimation. Journal of Machine Learning Research, 7(45):1231-126,, 2006.
Jayaraman J Thiagaraaan, Bindya Venkatesh, Prasanna Sattigeri, and Peer-Timo Bremer. Building calibrated deep models via uncertainty matching with auxiliary interval predictors. Proceedings of the AAAI Conference on Artificial Intelligence, pages 6005-6012, Apr. 2020.
Léonard Torossian, Aurélien Garivier, and Victor Picheny. Xaammed bandits: Optimizing quantiles, cvar and other risks. In Asian Conference on Machine Learning, pages 252-267. PMLR, 2019.
Vladimir Vovk. Conditional validity of inductive conformal predictors. In Asian conference on machine learning, pages 475-490. PMLR, 2012.
Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning random world, volume 29, Springer, 2005.
Junfeng Wen, Chun-Nam Yu, and Russell Greiner. Robust learning under uncertain test distributions: Relating covariate shift to model misspecification. In International Conference on Machine Learning, pages 631-639. PMLR, 2014.
Jonathan Wenger, Hedvig Kjellström. and Rudolph Triebel. Non-parametric calibration for classifica- tion. In International Conference on Artificial Intelligence and Statistics, pages 178-190, PMLR, 2020.
Jia Yuan Yu and Evdokia Nikolova. Sample complexity of risk-averse bandit-arm selection. In IJCAI, pages 2576-2582, 2013.
Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding deep learning (still) requires rethinking generalization. Communications of the ACM, 6403):107-115. 2021.
Mengyan Zhang and Cheng Soon Ong- Quantile bandits for best arms identification. In International Conference on Machine Learning, pages 12513-12523, PMLR, 2021.
Shengjia Zhao, Tengyu Ma, and Stefano Ermon. Individual calibration with randomized forecasting. In International Conference on Machine Learning, pages 1 11387-11397. PMLR, 2020.
13

---

A Related Literature
A.1 Estimating the whole distribution
The most straightforward thought is to estimate the whole distribution which accomplishes the task of mean regression and uncertainty quantification simuitaneously. Generally, there are two ways to fit a distribution: either fitting a parametric model or turning to nonparametric methods. Representative parametric methods Bayesian Neural Networks [MacKay 1992] and Deep Gaussian Process (DGP) Damianou and Lawrence 2013). of which the latter assumes the data is generated according to Gaussian process of which the parameters are further generated by another Gaussian process. The major draw backs of these Bayesian-type methods are their high computation cost during the training phase and statistical inconsistency under model misspecification. To resolve the computational issue, [Gal and Ghahramani (016) proposes a simple MC-Dropout method that captures the model uncertainty without changing the network structure. Despite the view of Bayesian, the task of distribution estimation also falls into the range of frequentists: by assuming the underlying distribution family (say, Gaussian), one can fulfill the task by estimating the moments or maximizing the likelihood. For example, the Heteroscedastic Neural Network (HNNN) [Kendall and Gal [017] gives an estimation of both mean and variance the final layer by assuming the Gaussian distribution. An ensemble method called Deep Ensemble [Lakshminaray yanshannaraaaaan et al. [017)] is constructed on the basis of HNN. (Cui al. 2020 designs an algorithm that finds the best Gaussian distribution to minimize the empirical Maximum Mean Discrepancy (MMD) from the observed data so as to ensure the marginal calibration property of quantiles. (Zhao et al., 2020] raises a randomized forecasting principle that minimizes the randomized calibration error with negative log-likelibbe likelihood regularization term to ensure the so-called sharpness. Although they claim the construction of (randomized) individual calibration (named PAIC) from the existence of mPAIC predictors, they do not provide any non-trivial mPAIC predictor except for the one that permanently predicts the uniform distribution.
The other approach to estimating the entire distribution is the nonparametric way without assuming the underlying distribution. For example, Kernel Density Estimation (KDE) utilizes kernel methods to estimate the conditional probability density function. Although the method enjoys certain theoretical consistency guarantees [Lei and Wasserman 2014 Lei et al., 2018 Bilodeau et al., [02T]_ it is often computationally expensive and not sample-efficient in high dimensional cases, which makes it impractical for uncertainty quantification. A similar idea is applied in [Song et al., 2019) via a post-hoc way. The method is to improve on some given distribution estimation via the Gaussian process and Beta link function. However, unlike the KDE literature, the method in [Song et al. 2019 lacks theoretical guarantee of its consistency.
A.2 Estimating quantiles
Due to the excessive and sometimes unnecessary difficulties in estimating the distribution, some researchers suggest estimating the quantiles only, since a satisfying quantile estimation suffices in many downstream tasks. The quantile prediction is now regarded as new task, where one natural idea is to combine the average calibration error on quantiles with a sharpness regularization term. [Pearce] al., 2018 appends a calibration error (with respect to coverage and sharpness) to the loss functions to obtain the calibrated prediction intervals at pre-specified confidence level. [Thiagarajan et al.) 2020] induces an extra interval predictor with a similar combined loss. But their methods require the pre-determined confidence or quantile level and the whole model needs to be retrained if one wants another level of prediction. Such methods in [Pearce et al. 20188] Thiagarajan et al. 2020 are only evaluated empirically without theoretical guarantee for even the marginal distribution. [Chung ctt al. (022] proposes an average calibration loss combined with a similar sharpness regularization term, which achieves average (but not individually, with these models that aim to give a quantile prediction out of scratch, [Kuleshov et al.) 2018 develops a post-hoc way of calibrating (on average) a pre-trained quantile prediction model via isotonic regression methods. The aim is to choose a specific level of quantile prediction that reaches the desired marginal calibration of which the consistency relies solely the pre-trained model in contrast to our model that applies a model-agnostic post-hoe calibration.
will discuss later) calibration. Along
Another line of research is to use Kernelized Support Vector Machines (KSVM) to predict the conditional quantiles in one step, which fixes the inconsistency issue (Takeuchi et al. 2006 (Stinwarrt and Christmann 2011 Essentially, methods from both papers search for quantile prediction
14

---

function over the Reproducing Kernel Hilbert Space (RKHS), say H. Theoretically, Takeuchi et al. [0066] derive a performance guarantee for the conditional quantile estimator based on the assumption of finite H-norm, which is in parallel with our I Lipschitz class (Assumption | (a)). 'he Lipschitz function class and the bounded H- H-norm function class overlap, but neither one contained in the other. (Stainwart and Christmann [2011) analyze the same algorithm but adopt a different approach. Instead of imposing a bounded H-norm, they impose a decay rate for the eigenvalues of the kernel integral operator, alongside some other assumptions on the underlying distribution. Those works focus on the RKHS space while we focus on the Lipschitz space. For the Lipschitz function class, to the best of our knowledge, our work is the first result to provide an individual guarantee. From the practical point of view, both papers rely on solving a kernelized learning problem that cannot scale well with the sample size or dimension, while our algorithm is simple implement and of low computational cost.
A.3 Conformal prediction
very closely related field is conformal prediction, where the goal is to give coverage set that contains the true label with a certain probability. Theoretically, the goal of conformal prediction is no harder than our quantile prediction task since correct quantiles induce desired coverage sets, while empirically different measurements are considered: the conformal prediction cares about giving the sharpest covering set rather than specifying which quantiles are estimated. For the regression problem, there is an "impossible triangle" summarized in a concurrent work (Gibbs et al. 2023] :
(i) a conditional/indiiidual coverage (as our Definition (3) (ii) no assumptions on the underlying distribution (iii) finite-sample guarantee and asymptotic consistency. Such triangle is shown by
Vovk [2012). Foygel Barber et al. 2021] to be impossible to achieve simultaneously. Some works in the conformal prediction literature aim to reconcile the impossibility: |Romano al. [2019] consider
the marginal/average case (as defined in Definition D. which greatly relaxes (i). Gibbs et al. (2023)
violates (i) in a milder but different way, focusing on a midway between marginal coverage and
conditional coverage defined by the linear function class.
The most related paper to ours may be Feldman et al. [202T]: both theirs and ours attempt to address the minimum cost of keeping the requirement (i). Before diving into the detailed comparison, we summarize the main difference first: Feldman et al. [2021] relaxes the remaining two, while paper keeps (ii)) and relaxes (ii) in a much milder way. Feldman et al. [202T] achieves a theoretical guarantee that the true conditional quantile is the minimizer of the regularized population loss under realizability assumption (that is, the true conditional quantiles lie in the function class considered) and an assumption that the conditional distribution is continuous with respect to features. However, its
theoretical results have several limitations. First, it achieves only asymptotic consistency compared to our finite-sample guarantee in Theorem I which means that [eldman requirement (iii). Besides, it only achieves a necessary condition that the true conditional quantile is the minimizer of the population loss but not vice versa. Second, Feldman et al. [202[] violates (ii)
al. [2021] violates
to much deeper content than ours. Note that the key assumptions made in Feldman et al. [221]] assumption and the continuity assumption. The former means that
are the seto-approximation-nnoe there should be no model misspecification at all, which is not a weak one in the statistical learning theory and limits the guarantee from a desirable "distribution-fred property. On the contrary, our algorithm can have finite-sample guarantee with only three very general assumptions summarized in Assumption T Furthermore, our Lipschitz assumption is made only to achieve a finite-same sample guarantee (of which the necessity is supported by Theorem (). while it can still be relaxed to a weaker one than Feldman et al. [202T] to reach the same asymptotic result. To see that, our result only requires the conditional quantile to be continuous, while Feldman et al. (202T] requires continuity for the whole conditional distribution. Third, the regularization term in Feldman et al. [202T) is related to zero-one indicator, which is discontinuous for gradient descent. The authors twist their algorithm by replacing it with a smooth approximation to overcome the issue without any formal guarantee on the twisted one. In contrast, our theory accompanies the exact same algorithm that we propose.
Our paper can also serve as starting point for understanding the empirical success of the "split" procedure that is common among the conformal prediction algorithms. As summarized before, most conformal prediction methods theoretically abandon requirement (i) with only a marginal coverage guarantee. But tricky point is that if one only needs the marginal coverage, it suffices to merely use the empirical quantile without the splitting step at all. Why does the learner bother to "split"? Empirically, the splitting procedure helps obtain a more homogeneous coverage across the whole
15

---

feature space, while our paper is the first theoretical justification up to our knowledge; to reduce the Lipschitz constant and smooth the target. Such an intuition is justified by Theorem I and Theorem
2 from both the upper and the lower bounds.
B Failure of Existing Algorithms in Individual Calibration
In this subsection, we will give brief discussion of the existing methods in regression calibration literature. We give a close inspection of the existing calibration algorithms, proving that they cannot achieve a desirable result (or even avoid naive and meaningless ones).
Some methods assume that the underlying distribution is Gaussian [Cui et al. 2020 Zhao et al. 2020]. The most straightforward drawback one can think of is that these Gaussian based methods are not consistent once the Gaussian assumption is violated. However, we will prove with some counter-examples that they cannot prevent the naive or meaningless prediction even if the underlying model Gaussian.
Cui et al. [2020) use a two-step procedure to estimate the conditional distribution, where at each step the underlying model is neural network predicting the mean μ(z) and variance o² (z). The first step is simply fitting the observed {(X;, Yo) into the model, getting an initial guess of and σ.. The second step utilizes the maximum mean discrepancy (MMD) principle. To be more specific, the second step iterates the following steps until convergence: first generate N new samples (Y:),N from the estimated Gaussian model N
(X.)), then update the model estimation (ūt, 8² to minimize the empirical MMD + Σx ((Y)) + E o(Y.)]]. But the problem is: the corresponding relation between each {(X;, Yo)} pair is wiped out from the model by calculating the empirical MMD. In other words, the procedure produces the same result even if one permute estimated conditional distributions {())),), (X.))) arbitrarily. This can lead to a failure in predicting the true conditional distribution even under the Gaussian condition, especially when they use multi-layer neural networks which are well-known universal approximators that can fit any continuous function HHornkk et al. 1989
Zhao et al. 2020 try to give a randomized prediction of the edf (in their paper denoted by hfz, r-).), where r is some Unif[o, 1] random variable). The idea behind this is that if one gets a good cdf estimation, then F. (Y) should be close to the uniform distribution. So the target is to minimize the empirical absolute difference between hfz, (y) and (denoted by I,). If L, < c with probability at least 1- 6, then they call it (c, -roobbbly approximately individual calibrated (PAIC). But minimizing L, alone cannot avoid being trapped by naïve and meaningless calibration: M/,, rl) = r. So they combine the loss with another maximum likelihood type loss L2 (negative logarithm of likelihood, by assuming the underlying distribution of Y is Gaussian). But combining that loss cannot guarantee a meaningful calibration. To illustrate this, one just needs to consider the continuous distribution case, where each X; Z; is only observed once almost surely. A global minimizer for the empirical loss L1 = - 1 (a, be N(u - (r:)) o²), where N' is the Gaussian and (--) is the cdf of standard Gaussian. can be easily checked that this Gaussian distribution's rithh quantile is exactly y:). The maximum likelihood type loss L2 is now 1/ log(o²) + "C__-biia
(a ry) - r.| is obtained if one predicts h(2,,r,](.) to
_/2 log(o²) + constant. Minimizing it means that the o² term will asymptotically go to zero, leading an overfit and meaningless result of hIr; ra](.) + ky.- Again, this overfitting problem can be exacerbated by the fact that they train the model via neural network, which is commonly believed be highly overparametrized that can fit into any data nowadays. We must also point out that although Zhao et al. (2020) claims that they prove the existence of an (e, 8 predictor, the existence is only guaranteed by the existence of another monotonically PAIC (where monotonicity requires hiz, r((x) is monotonically increasing with respect to r) predictor (see their Theorem 1). But they do not provide the existence of any non-trivial mPAIC predictor other than the trivial solution h(z (v) = r. Therefore, we think their argument least incomplete, if not suspicious.
For those existing methods that do not assume Gaussian distributions, unfortunately, they are still inconsistent from an individual calibration perspective (Kuleshov al. 2018 Chung et al. 202T].
Kuleshov al. (2018] develop an isotonic regression way to calibrate the quantiles in the regression setting, which is inspired by the widely used isotonic regression in the classification setting, hence is recalibration method based on some initial estimation of all quantiles (say (Ộ), (7).
16

---

Then the authors propose an isotonic regression method that learns a monotonically increasing function ((Fepecced)) = Tobserved- where Tobserved 1/4 Ex=1 1(x; I Qo The final output is Q = (r)- The major drawback of their method is that it heavily relies on the initial calibration Q. Hence the final marginal (or individual) recalibration result is only valid when the function class {Q- : E [)] contains marginal (or individual) calibrated prediction. Despite the marginal case, such a realizability condition is too strong to obtain the individual calibration set, leading to the failure of their recalibration methods in the individual case.
Chung et al. [ propose combined loss in their Section 3.2 with respect to both average calibra-
tion and sharpness. The first term minimizes the calibration error: L1 1 Σ 1(Y, S Q+)
7) L
n E I0, - -x,())) > Q (X.)) + 1 1[y - - -
(V) (x.)). Note that T1 is zero if the prediction Q- is marginally calibrated. Σ (X,) Q1 (X,), called sharpness regularization term. But that combined loss L = (1- A)L + XL2 will fail even for calibrating such a simple in- stance: Unif 1], u ~ 1ni and Y = ux. Any -Th( > 0.5) quantile should be (Fa),) = (T - 0.5),, while minimizing their proposed loss for E [0, λo] ends up with Q, = a for those x < T 0.5 and @ = 0 for those x 2 т - 0.5. For larger À E [Ao. 1]. minimizing their proposed loss will lead to an all-zero meaningless prediction Q = 0. The detailed calculation can be found in Appendix.
The second term L2
Finally, we point out another issue that using the training residuals as errors to estimate the error distribution will lead to a biased (usually over-confident) estimation. For example, the Section 3.1 of Chung et al. [002] proposes a combined method of first training a regression model f on {(x,, Y.)] and then train regression model on res, = YY, - (Xx)). But they use the same set of data to do both steps of training, leading to statistically inconsistent estimator. The inconsistency becomes especially severe for modern highly overparametrized machine learning models that can fit into even random data (Zhang et al.) [221].
Remark 1. Many researchers [Chung et al. 2021 Zhao et al., 2020) argue that combining a sharpness regularization term with the average calibration loss will help to force the calibration model to stay away from marginal calibration. However, in this section, we show that those combining methods fail to provide meaningful individually calibrated quantile prediction. In conclusion, a sharpness regularization term may help achieve a marginal calibration, while does no good for the goal of individual calibration.
C Discussions on Negative//mposibility Results for Individual Calibration
Before elaborating on those impossibility results, we adopt the lower bound result from the nonpara- metric (mean) estimation literature Györf et al. [2002]. showing that any estimator may suffer from an arbitrarily slow convergence rate without additional assumptions on the underlying estimation, even in a noiseless setting. The original result established for the conditional mean regression, and the same result holds for the conditional quantile estimator in our case.
Theorem 6 (Theorem 3.1 in Györf et al. 2002], rephrased). For any sequence of conditional quantile estimators (Qr... } and any arbitrurily slow convergence rate {an) with lim Un = 0 and a 02 2 >> 3P, s.t. X ~ Unif(o, 1), UIX 5o(X) (hence Q+ x) = =g(x), VT 0),
where g(X) E {-1,+11, and
E[1êra(X) -xxx)) 21.
lim sup D0
an
Such negative result further justifies our Assumption Tin way that one cannot achieve any mean- ingful convergence rate guarantee without making any assumptions on the underlying distribution.
Discussions on Zhao et al. [2000: The negative result obtained Zhao et al. (2020) (their Proposition 1) is for the target of estimating the conditional distribution's 's edf. The authors argue that the quality of an estimation fr (.) should be measured via how close the distribution (Y) to the uniform distribution Unif(0, 1]. The measurement of closeness is defined by the Wasserstein-1 distance.
17

---

Their counter- exanter-examle is constructed in a way that Pyjx- is a singleton (say, a delta distribution on the point g(z), where ((x) is chosen adversarially by the nature). They prove that any learner cannot distinguish between the true distribution and the adversarialyy chosen distribution from the observed data, thus any algorithm that claims a small closeness between (Y) and Unif(0, 1] cannot guarantee that the observed data are not sampled from the adversarial distribution. Since the Wasserstein-1 1 distance between a delta distribution and Unif(0, 1] is at least 4. the closeness claim is false for at least the adversarial distribution. We need to note that the measurement they take is a little suspicious: if the conditional distribution itself is delta distribution, then the closeness of Fz(Y) to the uniform distribution cannot be regarded as a good measurement of distribution calibration, since even the oracle will also suffer an at least 1 Wasserstein distance. Although a certain negative result when Y|X= 6g(X) is also made in our Theorem 6 we need to note that the closeness of (Y) to Unif (0, 1] can only be a good calibration measurement if the conditional distribution is continuous. On the contrary, for the general continuous conditional distributions (with only those mild assumptions), our methods can achieve a minimax optimal rate individual calibration guarantee with respect to any quantile (Theorem D-
Discussions on Lei and Wasserman (2014). Vovk [012]: Vovk 2012 's impossibility result (their
Proposition 4) is a version of Lemma 1 in Lei and Wasserman [014].. Thus we will only discuss Lemma in Lei and Wasserman [2014) here. The impossibility result established for those distributions which have continuous marginal distributions Px. They focus on a conditional coverage
condition, meaning that one should obtain (1 11-11 p) coverage set Č C(Xist) such that PYciss E
C(Xiest | X =) 2 1 p for any x E X almost surely. Lei and Wasserman 2014 prove that for any distribution P and its any continuous point To. for any Jinite n. such an almost sure conditional coverage must result in an infinite expected length of covering set with respect to Lebesgue measure. This negative result seems to be disappointing for individual calibration at first glance. However, an almost sure guarantce in a finite sample case is usually too strong in practice. Our results are established for two weaker but more practical settings; either the strong consistency in the asymptotic case as n - oo or mean squared error guarantee at the finite sample case. similar spirit of relaxing the requirement in order to get a positive result can also be found in Lei et al. (2018)) where the authors get an asymptotic result that the covering band will converge to the oracle under certain assumptions. Note that their requirement is to get the sharpest covering band, meaning that they need to estimate the probability density function, which leads to their more restricted assumptions and different analysis, compared to ours.
D Comparison with Other Results
In this subsection, make comparison of the technical contributions with other related literature such as concentration of order statistics and Kernel Density Estimation (KDE).
Current quantile concentration analyses are all based the i...d. assumption, which is not applicable to the continuous case conditional quantile estimation problem, where the learner almost surely does not observe identical features. But to complete our discussions, we give a brief introduction to the existing quantile concentration methods. There are two lines of research, which we call vertical and horizontal correspondingly. Vertical concentration guarantee is based on curving the convergence rate of the Smirnov process (the process of the empirical edf converges to the true edf) and applying the Dvoretzky- Kiefer- Wolto-wolf inequality (Dvoretzky et al. 1956 Massart 1990] to get sharp coverage along the vertical line (in regards to edf function value) that is, the get a sharp 5 >>>>> 0 such
that Qr < Q < Q with high probability (Szorenyi et al. 2015 Altschuler et al. 2019 But turning such a vertical guarantee that guarantees the coverage probability into horizontal one that guarantees the estimation error is highly non-trivial. By assuming the cdf to be strictly increasing
and continuous, Szorenyi et al. [2015), Torossian et al. [2019] Kolla et al. [0019] Howard and
Ramdas [2022] get horizontal concentration of order O(L/vn) Cassel et al. [2018] assume
Lipschitz edf (hence an upper bound for pdf) and get a rate of 0(1) / Vn). Altschuler al. 2019
directly assumes a function (which they call R) of the edf to transform the vertical concentration
into a horizontal one, Another line of quantile concentration is to directly get horizontal guarantee.
Yu and Nikolova [2013) use Chebyshev's 's inequality to prove the horizontal concentration while remaining coarse-graing -oarse-grained compared to Hoeffding type guarantee. Zhang and Ong [202T) utilize
an alternative assumption and analysis of Boucheron and Thomas [012 that relates to certain
18

---

type of distribution with non-decreasing hazard rate (or monotone hazard rate, MHR). Although the MHR assumption common in survival analysis, such an assumption is no longer appropriate for many other distributions, such as Gaussian. In our proof, we make a combination of the finite sample theory of the M-estimator and the local strong convexity of the conditional expected pinball loss to derive a new proof for the concentration of empirical conditional quantile around the true conditional quantile. We do not assume the upper bound or lower bound on the pdf, the strict increasing edf, or the non-decreasing hazard rate. We also generalize the i.i.d setting of quantile concentration to the independent but not identical setting, which is the first one to our knowledge in the quantile concentration literature.
Another important line of research is the kernel density estimation (KDE) literature focuses on giving a full characterization of the conditional distribution. Such a goal is aligned with the need for the full characterization of pdf some downstream tasks. For example, if one wants to find the sharpest covering set (or covering band, if one further assumes the conditional distribution is univariate), then one must identify the probability density function so that thresholding can thus be constructed [[ei and Wasserman 2014, Lei et al. 2018). But for simpler and probably wider class of downstream tasks, such as decision-making based on quantiles (newsvendor problem) and evaluating and minimizing the quantiles (robust optimization), an accurate estimation of the conditional quantiles should suffice.
From the point of developing practical algorithms, KDE methods are often computationally intractable for conditional quantile estimations. To get quantile estimation, one has to integrate the full probability density function. The whole procedure is highly costly even if one discretized the integration. Another point is that theoretical guarantees in the KDE literature are often established with respect to the mean integrated squared error (MISE), while it does not simultaneously guarantee convergence with respect to some quantile. Consider simple instance where the estimated probability
density function p(y) differs from true probability density function p(y) by a factor of O MM
The MISE converges with a rate of O (+). However, the cumulative pdf estimation error (denoted by 3n caused by the non-exchann ex-eehagggaabbiity between the square and the integration.
lØ(t) p(t)ldt may diverge (since Sa 1 dt is a divergent integration). The divergence is [2]
Since the final goal is to get estimation of the conditional quantile, we can relax the assumptions in the KDE literature largely. For those nonparametric methods such as Lei and Wasserman (2014), they require the smoothness (to be more specific, Höldemess) of the underlying conditional pdf with respect to y and the Lipschitzness with respect to r in the infinite norm, which a far stronger assumption than our case. Instead, only assume that the conditional quantile (rather than the full pdf with respect to the sup norm) is Lipschitz. For those methods that establish the KDE convergence results in generalization bound via complexity arguments, our assumption is still much weaker. For example. Bilodeau et al. [2021] make a realizability assumption (that the true pdf must lie in the hypothesis class) to achieve generalization bound, while our assumption just requires the true conditional quantile lies in a bounded set (-M, M] C IR
In word, our analysis combines the idea of decomposing the error into bias and variance terms as in nonparametric regression methods and the spirit of the covering number arguments in the analysis of parametric models. Such combination is novel and non-trivial, which is critical for the relaxation of the assumptions.
E Numerical Experiment Details
In this part, we provide detailed experiment settings for Section
4
E.1
Evaluation metrics
Mean Absolute Calibration Error (MACE): for (expected) quantile level 7 and a quantile prediction model Q(), define "observed quantile level" by
n
10x; < Q().,
n
19

---

For pre-specified quantile levels 0 S < ... < TM S 1, MACE is calculated by:
M
1
MACE( Y)) Q)
T5
obs
M
In practice, .M 100 and pick 7's from 0.01 to 0.99 with equal step size.
Adversarial Group Calibration Error (AGCE): introduced in Zhao al. [220 AGCE is calculated in practice by first sub-sampling (possibly with replacement) BI. B, S {(X,, H and then retrieve the maximum MACE from the pool:
AGCE(( (B Q) == max Q)).
Check Score: check score is an empirical version the expected pinball loss. For pre- specified quantile levels (I I 71 < Tm S 1 (in our experiments r's are set in the same way as in MACE), the check score is defined as:
CheekSccre(((X,, (.), Q)
Σ Pry (Q(X,T,),))
with p.
Length: our length metric gives the average interval length of a predicted 90% confidence interval. The interval is given as the area between the lower quantile (at level Flow 0.05)
and the upper quantile (at level Tup 0.95):
Length(((X) Y.) Q)
= n Σ (Q((, - Q(X, Now)).
Coverage: encourage fairer comparison, always attach the empirical coverage rate to the length metric. Coverage rate evaluates the percentage of real data that covered by the given 90% confidence intervals:
Now) :) HQQ(, Tow) SYY Y.S.S Q(Xi,
We make brief explanation for the selection of the performance metrics as follows. The former three metrics are typical measurements in the literature of quantile regression and uncertainty calibration that measure the quality of the calibrations, while the latter two are classic metrics in the field of conformal prediction. Although the quantile calibration task very closely related the conformal prediction task, the goals are still different. Conformal prediction aims to provide as sharp as possible prediction sets that cover the true outcomes with desired rates. While precise quantile prediction guarantees the coverage rate, such a covering band may not be the sharpest. As is discussed in Appendix B the additional sharpness regularization term could harm the goal of precise quantile prediction. For example, people may select 0.05 and 0.95 quantiles to construct 0.9 coverage band, but the sharpest covering band probably would deviate from such quantiles. However, high- quality quantile prediction alone is of independent interest for many downstream tasks, such as the newsvendor problem. This points to difference between the objectives calibration and conformal prediction. As result, the popular performance metrics applied in conformal prediction such as the average interval length and the coverage rate are not direct measurements for the quantile calibration problem. Nevertheless, other mcasurements, for example, the independence between the coverage indicator and the band length in Feldman al. [2011) are only necessary condition of the quantile calibration and are therefore not considered here.
still present the Length and the Coverage metrics for completeness. Some
F.2 Synthetic Dataset
For an example of heteroscedastic Gaussian noises, consider underlying generation scheme of (X, Y) R2 from X Unif(O, 15) and YJX (Noo(X), (f))) with [x( sinnin 2m X) and (X) max(0.2X (.11). We generate Ahrain 40000 samples for training and
20

---

Fest 1000 for testing. For our NRC method, we randomly select = 36000 for training regression model and Terain n, 4000 samples for the quantile calibration step. All neural networks used in this example have hidden layers of size [100, 50). As shown in Figure [ our NRC method successfully captures the actual fluctuation of conditional variance, while the traditional benchmark model Deep Ensemble which relies on 5 heteroscedastic neural networks cannot capture the individual dependence on the feature of the noise variance.
E.3 Dimension reduction implementation
As given in Theorem for high dimensional data (i.e., large d) the upper bound guarantee may get too loose theoretical guarantee. Motivated by this dilemma, we propose applying dimension- reduction techniques when doing nonparametric approximation. We summarize them as the NRC-DR algorithms. There are loads of choices for dimension reduction, in this work we only implement three of them: random projection, covariate selection, and second-to-lass layer feature extraction.
we use the Gaussian random projection method, in which we project the d-dimensional inputs on the column space of a randomly generated do d matrix (assuming that do < d, otherwise original input will be retained), whose elements are independently generated from (O, 5).
Random Projection Dasgupta 2013
Covariate Selection: we select the most relevant features from the input vector. In this experiment, we select dg out of d features that have the largest absolute value of Pearson correlation coefficient with the target variable.
Second-to-Last Layer Embedding: The second-to-last layer of trained neural network has been widely used for feature extraction when the input dimension is high. We add an additional layer of size da to the regression network f(-) in Algorithm 2 and directly use this layer as the feature embedding.
NRC algorithms implemented with each of these dimension-reductio techniques are called NRC-RP, NRC-Cov, and NRC-Embed respectively. In the UCI-8 experiments implement NRC-RP and NRC-Cov, with the full result presented in T while for high- nigh-dimensional datasets we choose the medical expenditure panel survey (MEPS) datasets [panel 19, 2017 panel 20 2017 panel 21 2017] as suggested in [Romano et al. 2019 to test out our NRC-RP and NRC-Embed variants. The sizes and feature dimensions of these datasets are listed in Table 3 Due to the relative abundance of training data, we set the hidden layers of neural networks to [64, 64]. the rest hyperparameters are pre-determined by grid search. For the results presented in Table we reduce the dimension to do- 30.
E.4 Benchmark methods
We explain the benchmark models and uncertainty quantification methods we competing against in Section 4 The first type of method is based on the Gaussian assumption. They assume that the conditional distribution of Y given X can be fully determined by the conditional mean and variance. These models are of the same output structure of two neurons: one for prediction of (XX), another for (X). At the implementation of the experiments, (X) share the same overall network structure, with two hidden layers of size (20, 20] for small datasets and [64, 64] for large or high dimensional datasets. Amongst them, MC-Dropout (MCDrop) Gal and Ghahramani (0066), Deep Ensembles (DeepEEnemmlle) Lakshminarayana
assume that (xX) and
et al. 2017
and Heteroscedastic Neural Network (HNN) Kendall and Gal [017 are trained by minimizing the negative log- ligeliieldhood (NIJ.) error. The Maximum Mean Discrepancy method (MMD) Cui et al. (2020) trains by first minimizing the NLL loss of an HNN model, and then re-train the model to minimize the MMD between the model prediction and the same training set at the second phase. The method introduced in section 3.2 of Chunn al. [202T] introduces a combined calibration loss (CCL) with an average calibration loss and a sharpness regularization term, and the model is trained by minimizing it. Apart from these ferd-rd structures, there is also Deep Gaussian Process model (DGP) Damianou and Lawrence (2013) which assumes that the data are generated from Gaussian process, and is trained by the doubly stochastic variational inference algorithm [Salimbeni and Deisenroth [0177] The DGP model is of high training cost, we thus only test it on small datasets.
(HNN))Kendal)
We also run experiments on host-hoc calibration methods, such as isotonic regression (ISR) suggested by Kuleshov et al. 20188. In our experiments, we evaluate the post-hoc ISR method by attaching it
21

---

to pre-trained HNN model. We also implement another post-hoe method named Model Agmostic Quantile Regression (MAQR) method given in section 3.1 of Chung al. [202T]. Their algorithm is similar to ours in terms of training regression model first, but different in two aspects: first, MAQR uses the same dataset to both regression training and quantile calibration; second, MAQR obtains quantile result at the second phase via training new quantile regression model, which in our implementation is a neural network. For the algorithms coming from the conformal prediction society, we test the Conformalized Quantile Regression (CQR) method introduced by Romano et al. [2019] and the Orthogonal Quantile Regression (OQR) approach from Feldman et al. [202]]] The CQR algorithm also takes training step followed by calibration step, with the training step learning the target quantiles by minimizing the pinball loss and the calibration step further calibrating the quantile predictions by conformal prediction. Based on the CQR paradigm, the OQR algorithm further improves conditional coverage by introducing an additional regularization term to the vanilla pinball loss used for quantile regression. The regularization term encourages the independence between interval length and the violation (falling out of the predicted interval) event of the response variable.
E,5 Experiment details on the 8 UCI datasets
We have run extensive experiments to compare our methods against multiple state-of-the-ar bench- mark models. The full results are listed in Table A brief review of all the benchmarks has been provided in Section E,4 and realizations of algorithms NRC-RP and NRC-Cov are explained in Appendix
short summary of the 8 UCI datasets is listed in 3 For all the network structures used in Section we fix the size of the hidden layers to [20, 20] (and [10] for the DGP model). For dimension reduction algorithms, the target dimension that we reduce to is set to 4. The rest hyperparameters (including the learning rate, minibatch size, kemnll width, etc) are pre-determined by running grid search. The learning rate is searched within [10-2.5 5 x 10-3, 10-3), the minibatch size is searched within [10, 64, 128), and the rest hyperparameters exclusive to each model searched in the same fashion. For each possible combination of model and dataset, we save its optimal hyperparameters into separate configuration files.
During the main experiment, for each combination model and dataset, we load optimal hyperparameter configurations, we randomly split out 10% of the whole sample as the testing set. On the remaining set, for our NRC algorithms additionally separate out 30% for recalibration, and for the rest algorithms, additional partitioning is required. The rest data is then fed into the whole train-validation loop, and up early stopping scheme with patience count of 20. That is, if for 20 epochs decrease in loss is observed on the validation set, then the training is early stopped.
which we repeat the whole training process 5 times. Each time
Table 3: Descriptions of Datasets. Each row gives the number of examples contained in the corre- sponding dataset and the number of features in the input variable.
(a)
Dataset # Examples
# Features
(b)
Boston
506
13
Dataset
# Examples
# Features
Concrete
1030
Meps_19
15785
139
Energy
768
Meps_20
17541
139
Kin8nm
8192
8
Meps _2
15656
139
Naval
11934
17
Bike-Sharing
17379
11
Power
9568
Wine-Red
1599
11
Wine
4898
Auto-MPG 392
Yacht
308
22

---

Table 4: Full table of experiment results 8 UCI datasets
Durse
HNN
MCDsp
DeeyErocable
CCL
ISR
DGP
MMD
MAGE
0055-884
0.11:000
0.11
Boatn
ChookSexxe
AOCE
02
0,2280044
192022
0,275007
1203
0.342501
1.1204
0,272818
1,130,3
04150009
1.2254
(921199)
16611999)
65=100%)
(009))
T9=3059))
(0:19999)
MAK.E:
0
0.112000
005422022
0.1720022
MAA2
0220004
0.184903
0,22883
B 10
0.24
CheckSeree
22306
361012
19:03
23:08
25:08
17:00
3951
MACT:
B MIC
CONS 18.04
OUSS 1803
0 DOOS
Eacrgy
AGESI
1910001
1.2
1010
V149975)
72106945)
2,19926999
MACR
01s--om
CO77-18.0
01--000
DO1S-0
0.13-2101
ChalkScren
00014000
003 0-00
101000000
IUZ2 4,003
0.76 (CCPO8T5)
05-6000995)
0.75 027-16600077)
025-20000009))
0.7% 000011993
D 12 2000919)))
0U57 BITUM65)
Keol
ChelkScore
AGCS
00005-5000000
0,175000
4008
01560005
000644-00002
05
D0001 m935
0.22601
0 000075-80000
D:0daa0 009
0 1006600000 0006
0.15
D DOKG-O 8X37
02:000055
TenghWCieengg)
0011-0/00099))
002-00000
RO1S50 1003/1003)
(008000))
D005-0 (042576)
DO05-0
0.004
776)
MAAE
00454062
07 D7
DOTE
013
0215
Prour
CheokScoxe
AOCE
01:002
1.30009
0.0
N
0075e0000
1
0.1 05
1.5402
1
0.16200007
300 003
0.33-2655
1.78112
LewgleCweergee
3079:
(6-35)979)
34 1090))
9/90))
17
MAKE
0106
R057200000 100
ROIG20007
Wine
CheckScoxe
MACE
0262.01
01200
6112005
0.122000 RO
01
01200000
0.345902
22002)88))
2.2000006
2350
MALE
0.11 10101
AG422
0.29
Loverne
KI5200%)
102
18 1099)
10:0
53 2795))
2222095
Datzatt
Mosic
MAQRL
OQR
OUR
SRC-RF
SRC-HP
Ciww
MACI:
0851 L862
ORSS 1802
BJA4K 6882
BLOS 1110
ILDSS L0122
Redto
ChockSeres
AC4S:
1
11100
6.2 8.1
LO
1.100
1
1-2104
0.26 LG
1.1003
8.99 LURS
MACE:
0019-1665
0n99 0003
0.46-11.05
0008-1000
00ET-D00
AGE
1S10
035-0.05
011-005
0 LUB
04-0044
0240006
Langhitwwrggs)
(T121799)
7010 70RE99
21 0088555555)
19: 909555
21 15995)
21 11995
MACE
0063440022
0115-0003
01
0066000
0042000
004-0001
CholkScore
06603
031-02
18:005
0.60-02
0.991-03
LengheCowrnge)
62= 1545)
55-0
3760
4. 1589)
550 16999)
85% S0115)
MACE
0017=00007
Karlan
CheckScone
BUCHA0 001
0012:002
0075-0002
0009
6009-0001
RU62200 002
0 D2
DO 001
006T-002
0.0020000 003
LagghCCeengee
(773))
025
0000000)
028-000005550 (006
0.7700 0036009)
004319022
0.242022
D012500005
Nard
ChookScoxs
AOCE
0 0003
08512001
0,09
00002200001
0.3720J
0002200007
0852220005
000005500000
0002
BO6S
0U00 1003
CovennC
001550235555
(929))
(006119)
(X2/915)
MACE
0829-0022
016720007
OJHN 10865
RC1620005
Powex
Checkheree
AGASE
0161000
12:803
110000
111:008
122000
11010
1.32100
22:0
MACI:
0639 1862
BO06 DLDIT
COS DOCOS
BOIZLO 1MS
SJIEe B41
OBET LAGDS
AGESE
6.151005
03271 1802
DLOTE 10.00
ILIPS LO
2510-19955)
2.2 (022(99))
2.311220055)
2 (833005
MACK
011100
0106--005
0.087-1004
ILDWS-OOS
ChackScoen
0.521019
0119-0007
041-033
0444-01
ILIS-IUIS
0S5-102
D44-117
535 1077))
17-004(7))
161000045)
53-199933
10--1915)
4.5. 1975))
55-16815)
E.6 Experiment details on the Bike-Sharing dataset
The Bike-Sharing dataset is a time series dataset that records the number of bike rentals as well as the value of 11 related features at each hour. As briefly described in Section 4 for this experiment we design the LSTTM-NNN and LSTM-NRC models, each having 2 layers, 40 hidden neurons for each layer, and a sliding window of size 5. For the LSTM-NRC model, in the recalibration step we are required to calculate (as part of the whole kernel matrix calculation) the distance between to input variables, and we retrieve the input vector by concatenating all the features within the history window (into a 55-dimension vector). As the dimension is high, we adopt and use random projection to map the concatenated feature vector to a 4 4-dimensional space.
redductionn techniques
Most experiment details, in terms of hyperparameter selection, early stopping criteria, etc, as similar to the settings described Appendix E.5 One difference lies the splitting of the dataset. Due to the temporal property, we are only allowed to predict "the future" from "the past", thus we split the starting 80% of time stamps for training (and further split for recalibration, if required), and the ending 20% of time stamps for testing.
23

---

E.7 Robustness to Covariate Shift
Covariate shift [Quinonero Caainonero-Candlla et al. 2008 describes a frequently observed phenomenon that the distribution of the covaniates (also called the independent variable, feature variable, etc) that differ between the training and testing datasets. In the presence of covariate shift, calibration model attained from the training set that only achieves average calibration (i.e., naive predictions of quantiles) may no longer suit well for the testing dataset, while an individually calibrated model maintains its optimality in terms of all the metrics mentioned above.
Following the experiment design in Chen both real datasets and half-synthetic ones (all of which are available on Dua and Graff (2017)) The experiment design almost the same as in Section except for difference the processing of testing sets. The experiment results on subset of benchmarks are presented in Table 5 which further validates the robustness of our algorithm covariate shift, especially from the dominating performance on the half-synthetic datasets. A few remarks on the datasets: the Wine dataset and Auto-MPG dataset both induce natural interpretation of covariate shift, with the training set and the testing set having different colors in the experiment on the Wine dataset, and different origin cities on the Auto-MPG dataset; and on the rest two datasets we synthetically simulate a covariate shift for the testing set, as recommended in Chen et al. (2016). To be specific, we introduce a covariate shift into the testing set by following steps:
[2016] and Wen al. [014]] experiment on
1. Randomly split out 10% of the whole dataset as the "testing pool".
2. On the training set X7, Calculate empirical mean and variance j XT, E Cov(Xr).
3. Sample X (á, Σ).
4. Resample (with replacement) 1000 examples from the testing pool, following density (Xess 0.35.).
Table 5: Experiments covariate shift datasets. The first three rows real datasets, while the last two are half-synthetic ones. For the Auto-MPG dataset, models are always trained on data samples from city 1, and tested on dataset of either city 2 city denoted by "MPG2" and "MPG3" respectively. Our algorithm has an outstanding performance, especially so on the two half-synthetic datasets.
Dataset
Metrie
HNN
MCDrop
DeepEnsemble
ISR
DGP
NRC
MACF
0.32+0.1
0.22+0.07
0.32+0.06
0.31++0.09
0.26-0006
0,2++0.03
Wine
AGCE
0.33++..1
0.23++0.07
0.34+00.6
0.33++0008
0.28=0005
0.23=00.3
CheckScore
0.51=002
0.29+4.06
0.43+0.11
0.48+00.2
0.31=0.06
0.31=0.09
MACE 0.065+0.03
0.12+0.01
0.07+0.03
0.059++0.02
0.17+00.1
0.069-0.008
MPG12
AGCE
0.21-0.05
0.22-10.44
0.1910.03
0.21+10.06
0.30=0.05
0.29
CheckScore
1.2:10.2
1.S.10.2
1.1-0.04
1.110.06
1.7:1.22
1.2:10.05
MACE 0.1+0.02
0.16+0.02
0.11+0.02
0.1+-0.03
0.23++0.03
0.08=0.03
MPG13
AGCE 0.240.009
0.2400.3
0.2-0.02
0,2410.06
0.37=0.07
0.26:10.1
CheckScore
1.1:00.44
1.6:0..3
1.2-0.08
1+:0055
1.4:10007
0.9
MACE
0.110.04
0.1710.09
0.12:10.33
0.059-10001
0.19
0.05
Concrete
AGCE
0.13=1.04
0.18110088
0.1540.04
0.093-10.22
0.21
0.07
CheckScore
4.11:00.2
4,8+00.8
4.3+00.2
4+-0.3
0,37
2.5+
MACE
0.24:10.88
0.33
0,2-0.01
0.11100.2
0.08
0.05 0.02
Boston
AGCE
0,26=0.07
0.34+0.11
0,21+0.01
0.14+0.02
0.1+0.03
0.075+0.02
CheckScore
2.1+0.7
3.2+1.2
1,5
1.1+0.1
0.71+0.33
0.61
24

---

Some Basic Probability Results
In this section, we present several well-known probability theory results.
Lemma 1 (Bernstein's 's Inequality). Let :n be N independent random variables. Suppose & < Mo almost surely. Then Vt 0,
N B(E,)) 21 S exp E N
2
N
Var(E:) + 1 Mot
Lemma 2 (Hoeffding's Inequality). Let En be N independent random variables. Suppose a, S E S b, almost surely. Then Vt 0,
N
21
P
Σ
(E.
ElE,)) 2
S exp
Σ i=1 (b, a.)²
21²
P
(B[C.) 21
2t exp
Σ Σi=1 (b. a.)2
Lemma (Poisson Approximation). Let Sn Ei= En.ks where for each n the random variables Sn,k, k E [ni] are mutuallly independent, each taking value in the set of non-negative integers. Suppose that Pn,k = P( = and 6n = ( 2) are such that as + 00,
(a) Σx_1 Pr.k + X < 00;
(b) maxke(m] (pm,k} 0:
(c) Ei=1 En, + O.
Then, Sn Dg (oisson(A) of a Poisson distribution with parameter as n - co.
G Proofs in Section 2
Proof of Proposition 7 From the fact that Y = = f(X). replace the Y term in each probability with U f(X), proving that the LHS and RHS probabilities are identical.
Proof of Proposition 2
By Lemma we know that
(x(x) E arg min (u, U)| z] WER
Since (+1,, = Ir (wn + c, u2 c), Ve E R, we have
(x) arg min min DE 14 f(x). =]], RCR
which completes the proof.
H Proofs in Section 3.1
We start the analysis fixed point X. We first analyze the case where n(x) 0, iec. Di [n], - x|| < We denote those ||X;, - x| < by (Xi, [m(z) We denote Feix,, by FR. We denote the average distribution of F& by = Exlr) .( Fi. and denote its -thh quantile by (F). We abbreviate the estimation QSNQ obtained by Algorithm to be Q in this section. To bound the difference ( QU .r), we factorize it into two terms, as analogy to the bias-variance factorization nonparametric analysis:
(=)
Q - = =)| - [,)
(x)
where the former term is called the bias term and the latter the variance term. We now deal with them one by one.
25

---

Lemma 4. Under Assumption (a) with Lipschitz constant L, we have
(Q-(F) - QU( = =>) - Lh.
Proof of Lemma 4
By Assumption (a), vk E [n(x)),
1Q+ X Q-(UIX = x)) & Lh,
which implies that
Fx (1O =x) + +b))* -U Lh) 21,
(0(0 ST.
Hence,
n(w)
(10+(UX==))tIM)) -(Q-(UIX==)ttt x) Σ =
= +b)) 2T.
n(x)
n(x)
2) x))
Σ n(x)
=) Fi (t) Lh) ST.
which justifies that
- Q = Lh.
Cemmaaccntronn the bias term. For the variance term, we need a more careful analysis. We start with inspecting the properties of the expected pinball loss function. We will show that under Assumption (b) and (c), the empirical risk minimization solution obtained at the second step of Algorithm I Iis actually close to the true quantile of by two steps:
Step 1.
The empirical risk minimizer Q (ar) concentrates around its corresponding population risk's unique minimizer for sufficiently large n n(x) C'p.r. M), where the minimizer denoted by u".
Step 2. The unique population risk minimizer u² in Step is identical to the population risk minimizer of [u (u, U.
We prove the latter Step 2 first due to its simplicity. The following Lemma 5 proves Step
2
Lemma 5. We define the population risk stated in Step by [=(x) == EU E> miz 1r Uk)
Suppose Assumption hoold for some L.P. M. Then
(a) I is twicely differentiable and convex;
(b) I, is identical By=Fl--u, U] up constant;
(c) For sufficiently small h < z.z 3ro 2 5. s.l. the minimizer of L, is unique (denoted by u*),
u² e [-M, M). and I, is p-strongly convex within a ball of radius To around its minimizer u*.
To prove Lemma 5 we introduce the following Lemma o whose proof is postponed to Appendix !
Lemma 6. For any distribution with edf F(u) = P(U S u), the expected pinball loss with respect to F is derivvavve with respect to u, and
θ E Eu-rl U)] = F(u
0 E. (- U)] = F(u) TT
Furthermore, for those F(u = F(u), the derivative exists.
Proof of Lemma 5 By Lemma 6 we have
1 n(w)
1 n(x)
v,I(u)= n(x)
V.E.AA
[(, n(x)
Σ
Fu(u)
k=
26

---

Since we assume that (X, U) follows a continuous distribution, and the conditional pdf exists, then
1 x(())
(1) (1x) Pa(u).
where pk is the pdf of Fk. vk E [n(x)], which proves part (a).
For part (b), one only needs to notice that F has the edf the form En(2) m(x) Fe. which implies that Fu (u, U)] has the same derivative as by.
For part (c), by Assumption ! (c), we have (xx) balls B(Q, = Xx)),) (.)))))) k E [n(x)], such that (u(u)= 2 p. Vu (Q, x.)).)).
At the same time, from the proof of Lemma we know that
Lh $ 1
[Q+(F) - Q-UX-------------- x
r,
which means that
max max(Qr = I
-XX)) -
I min QU = -ix)
(F)
S max Q (U = )in))
1
S &
min Xx,)) + 2
Such an inequality implies that for at least a ball of radius 5,
n(x)
Flu, U) T(x)) k=1
E (k
22.
AUE B(Q+(F),
which proves the uniqueness u² and the local strong convexity.
Since min, Q- = Xix) S Q+(F) < max Q+ = X and Q+ = Xi,) E [-M, MI, have
u' = E
Now we go back to prove Step To get the optimal convergence rate, our arguments share the same spirits of Bartlett et al. (2005), Koltchinskij (2006) of getting fast convergence rates for the M-estimators of low variance. But for simplicity. we derive our proof of uniform convergence via Bernstein's 's inequality and covering number arguments, which is more similar to the way of Maurer and Pontil 2009].
We are now interested in
(u) (u, U.2) - (Q, A direct computation shows that
Yue (-,,-, Un ~ ~ Xixs Ke [n(x)].
Lemma 7. Vu [-M, M]. have
(1)) I max(T, 7} (1u [ QF)
1 Var(Ex(u) <
Proof of Lemma
Note that the pinball loss can be expressed in an alternative way:
4uu, U)
- ()) (U- u)),
((Q,(F),U) = max
From the fact that - -1--00100000011111--------------
(F)
and that IT - u) -U - Q-())) ( (T(x)Q-(FF), we prove the first statement via Lemma 15
For the second claim, one just need to notice that any almost surely bounded by (a, b) random variables have at most variance of 1 (b a)
27

---

Combining Lemma 1 and Lemma
we directly prove the following
Lemma 8 (Pointwise Convergence). 6o 0, 0,f()) 1 the following holds
0, u E I-M, M) with probability at least
m(2)
n(a)
[[Ex(u)]
Σ Eu(
I
2
[ - QF)( (FFF)llog log 1
((x))
3n(xx)
1=1
1
+ V 3n(x)
[u (-) log
8
Proof. Direct corollary from Lemma 1
and Lemma 7
To reach a uniform convergence guarantee from pointwise convergence, we utilize the covering number arguments.
Definition (e- (e-Covering Numbers), Ve > O, a function class F is said to have an e-covering number of
(c, 11,0) = Fool: Fo C CF.F EF, E Fo, s.t. |l follso I e)-
We compute the function class formed up by the functions (u, (F). Lemma (Covering Number). Assume Assumption / (b) holds with parameter M. Define
M.M = (u, (F)) u [-M, Then VE > 0,
N (c. 11=0) diam((- M, M)) 2M
Proof of Lemma 9 By a similar argument to that Lemma 7 we have
(u,,-) - (FQ. (w2, )) L, (Q. )) s |u u2).
where the left-hand-side the function infinity form by itself.
Combining the above fact with another that
jui u2! S sup
1ug - u/ll -iiam([-M,
of wj. M.M]
complete the proof.
Since there will be many coefficients dependence, we abbreviate any constant C' that depends polynomially some other factors (c,,'., CGg) by C(ei,... C.y) Now provide the uniform convergence result by combining the pointwise convergence and the covering number.
Lemma 10 (Uniform Convergence). Ve, 8 > O, assume n(x) is sufficiently large such that nı(x) 2
log( c5
2M we have with probability at least - 8,
n(w)
1
2
2M
E [E)
(x(
S 4e +
|u Q (F)
log
(x)
T(
3n(x)
2
2M
[u - (F)) log
3n(z)
cd
holds for Vu [-M,
Specifically, if we select n(z]] m(a)s
we have vn(x)= 2 C2 = 6 6llo(mmx (2M, 10})
((x))
n(x)
(4 (4+ 8M log(2n(x)M))
B[ca(u)]
n(x) Σ {u(
I n(x)
(xx) k=1
4
Iu Q+(F)) 2 log(2n())M))
V 3n(x)
holds for Vu M..
28

---

Ve > 0, we can select an e-cover of the function class L denoted by Proof of Lemma 10
= -
Note that
(u( = L (u, Uia) - (F), I Vu E I,
which means that Bu,= E S, s.t.
(Ex(u) - E (2y)1 = [L( - [U,,)(U [|][4] - L| S.
Now vo > 0. we apply Lemma 8 on {Ex (u,) 15 E S with probability guarantee of at least
N(c.C M, we have Vuj S,
m
B[Gx(u))]
m Σ (n(z)
m
S
2
11 Q-((F) log
N(e, M 11ba)
m
3m
2 [1 Q. (F) log
[N(e, [L ll (la0) (lmp))
+
3m
Combining above inequality with the fact that - (u;)) S and the way we select the cover Lemma | such that Iu - -úyl S e, we complete the proof.
- Ex
To verify the specific conclusion, one just need to notice that our choice of n(x) ensures that n (x) (2 for n(z] n(x)M
Lemma 11 (Generalization Bound for Q,) If take (x(x) = arg min, Exel i Unr) as we do in Algorithm 1 then for sufficiently large m 2 C2 -,, .M), we have with probability at
least 1 n(x)af
Calp M, M, =0 - (111)) (n(x)
- S
n(x)
where C33 == 8 +
8M 3
(og(n(x)M) + 16log(n(x))M)
3p
More specifically, with probability at least n(x)
- Calp-1, M, log((())M) = (n(z) (V((To
Vn(z)
where Ci = /Cs 2
Proof of Lemma 11 (x)x) must fall into the neighborhood B(Q+(F). 5) stated Lemma 5 (c). fact, one just needs to verify that
We first claim that for sufficiently large n(2) > Co. the empirical minimizer
8M log (2((2))) + Vn(a)
16 V
M log(2n(z)M)) s p-r²
4+
nxx)
3
3
8
which will definitely hold for sufficiently large n(x) > C2(p 1 7,
From the fact that minimizing Exer) mea L,, Ui,) is equivalent to minimizing it minus some constant, we know that
n(z) Σ
n Ea (x)))) s O.
n(x)
A=1
29

---

By Lemma and Lemma 10
we have
0 s (=)(((x)) 1 (Q-
n(x)
Σ
ElEx (2)(z))
n(x)
kml
$
n(x)
1
+
SM
+ log log(22 (2((2n())A))
#1
|
3n(x)
16
(4-(x) Q-( (og(())))
<
4
SM
log (2((2
+
16
2
(1)(((()))) log(2n(x)
(())
3n(x)
where the first inequality follows from Lemma last from (c)
(a), the second inequality from Lemma 10 and the
Solve the quadratic inequality
B n(xx)
2²
w,
n(x))
where a = 4+ SM 3
log
(2g(((()))
we have
3p
8
B
4a
a(a)
1
(VB+Va).
w
Vn(z)
where w exactly
Hence
1
16M log 1og(2n(z)M) +
64 loloo222 1og(2n(<z)MM)
((Q-()))-T(((-(33) n(x)
(23 + 2a) 8
nxx)
3p
The latter conclusion follows from the local strong convexity of 4 again, as is shown in Lemma (c).
Lemma proves our S finally. We can combine Step T and Step2rin now to get an upper bound that is critical for the following analysis:
Lemma 12 (Bounding (2)xx) Q+ = =x)/² Suppose Assumption holds. For sufficiently small h < r ZL. have
(2,(z) - QUX x)) S Lh C23M
Ca Vn(x)
+ +1fn(x) C2) M Bounds in Lemma 77 fail] +
where C2 C2(p-1,1-1, M), C2 Ca(p M, log(n(z) in Lemma 11 and the failing
probability Lemma
is no larger than m(a)M provided n(x) 2 Cz.
Proof of Lemma 12 A direct corollary from Lemma and Lemma 11
Proof of the pointwise consistency part of Theorem Since the Ô
term and the failing
Vn(z)
probability term in Lemmal 12 tend n((x) - oc, only need to verify that select
such that as n - oo, we have
(a) h + 0:
30

---

(b) AC > 0,F(n(z) 2 C) 1.
We will show that both conditions are true for any wo where the pdf p(x) is positive and continuous at x = Tg. we select h = Cşn कमッ
The first argument automatically holds as + co. We turn to prove the second one.
Since p(x) is continuous p(ro), we can find sufficiently large No. such that
p(ao)
p(a) 2
2 Va (am,, h).
Then the probability of 1(x; E B h)} can be lower bounded by
P(X; B(xo, h)) p(au) vha = Cen पे.,
where V the volume of the unit ball X C R C. = p(aa) VCg. AC > 0.j E N_, we can define sequence of parameters {Ay} such that for any Poisson(Aj).
2
P(n 2 2 j j=1,2....
where we assume without loss of generality that Xÿ is non-decreasing.
For each j, we can define an auxiliary sequence independent random variables {v for
I -= [g]. 1,2,.....3 ........ such that each random variable Vi,q. Bernoulli( A ).
By Lemma
3 can find sequence of { 1 such that
P Σ Vlass 20) 2
p(p P(n 2 C) - 2--21-1 -
2 4q qj. = 1,2,....
where we select q, without loss of generality that A %
is non non- con-increasing.
As {1fX, E B(zo, K-) is also an independent Bernoulli random variable sequence with param-
eter Cgn कें we can select a non-decreasing sequence of {N₂} such that Vn Nj.
A n 2 l.
CBn
Then Vn 2 Ny, we have
n(x0) = 1 B M) 2 2
2 P
Σ Unny 2
2 21- j'
1=1
where the first inequality follows from the result of stochastically dominant of 1fX, E B(rohh) over Vi,n.j.
Ve 0, we select sufficiently large j 3. then Vn. 2 max (No, N₃).
F(n(zo) 2 C) 2 1-6,
which completes the proof.
To prove the mean squared error part of Theorem variable that will be proved later.
we state a lemma on the Binomial random
Lemma 13. For any Ç ~ Binomial(n, p). we have
[ 03} < (n+1IF'
(a) E
(b) Vr E N.r np, we have PP( <r)
(n-r)p
(mp-rfr-
31

---

Now we return to the proof of the mean squared error part of Theorem
Proof of the mean squared error part of Theorem have
By direct computation from Lemma 12 we
E
(0ie)))))
Q- = )
-xxe
E Exu=
BBin (Xrsst) - QU Xresa) [9]]
<
Xxx
Ex.m
2
4C4(0-1, M, log(((XXes))MM) 4M²
+
n( n(Xiea)) ) ((n(Xws) 2 C2}
n(X
+ 4.M2111 ( < Cab
Cilp M,M log(M(M ((Xees)
I 42 Ex...
Exi
1 (nn(Xws)) 2 C2] 4M1 (n(us) <
where the first inequality comes from the fact that (a + +++++++ d) < 4(a² + + + dP), and the second from that (XXres) S n. Note that we replace the original C4 that depends on the logarithmic term (Xeea) with a new C₁ that does not depend logarithmically on n( (((eu) but logarithmically on the entire n. From the definition of C4 in Lemma 11 we recall that the dependence on the logarithmic term is polynomial.
By Lemma 13 (a), we have
Exi..
(n((ieu 2 C2)3
Exi..
1
1 1fn(Xicn 01}
n(Xiea)
n ((Xua)
2
(n 1)) E B(Xistt h))
We try to bound Ex ((X))
Since X' C [0, 1]d, we can generate an cover of X, denoted by S (21) where
C
Ca S hd
where C' is some intrinsic constant with respect the space Rd equipped with metrie || ||-
Then
1
= 1 P(B(,,)) dp
Ex
P(X. E B(X,
Ca
1fx 1 ()) dp
j=1 x P(B(a,h)
Ca
11 1) P(B(I,, h))
Σ
1=
C
= Ca
hd
where is the probability measure with respect to X, and the first inequality follows from that S is 2 coover of X. the second from the fact that B(z,, 4)) B(x, h) for those || 2z S 4. and the inequality from the fact that f 1f E B(x,, 4) = P(B(z, h)).
now try to bound Expi Exi... [1sn(X c21l]
32

---

(n,(((zzhh) 2 2| then by Lemma 13 (b). we have
n C|C2]
4
P(n((a) C2)
(nP(B(z, h)) - [C2])² nP(B(z,h)
If EnP(B((, < 2[]] then for sufficiently large n 2 4[C2] - 1, by direct calculation we have,
P(B(a, h))
S 2]]] n n-1
2
$
[1
4
C2]
3
[C2]
4
[C2] n 1 1C2
[C2]
P(Bx, h))
n-
n2
n n-3
which implies that
ch (1- P(B((())))) (P(B(r, (
22 2 [c2] 0 (1 P(B(e,))))
4
We [[C]]]]
where C₁ is the number of combinations of (n,1). Hence
P(n(a)) < C2) s P((n(x) $ (C2l)
1C2]
Σ C² (1- ((B(s, n-1 ( (P(B(z,)))
t=
8(C2)²
(1-P( (i))
3
8(C2)²
exp( Ppo B h)))
3
Since exp(-a) S max; (b (exppxb)) b)) ea we have
C7
P(n(x) < C2) < nP(B(e,h))
where C7 8(Ca)) ) 3e
Therefore, we can conclude that (w.l.o.g, we assume C7 2 4)
ma C
P(n(x) < C nP(B(a,h)) nP(((z,h))
Hence, we can finally give the bound by
E (Xuss)
Q = [2
Capp-1 ,M, lo,.o)(MM) (n(Xicu) 2 C2) + (n( < Caa
< 14T.212 + Exmi
EXIO
n(Xusa
1 <TT + C,Ch-
+ C7Ch-d." 1
n
n
1
= =4L² + Cg x-- log(n
For the special case of = 0, we can select = (MM), hence reaching fast rate of convergence of order O
If T. > 0, by selecting h पtz T (dC)) गें we get the convergence result that
(F < L 2 (L. -1, M, Mog(n)) /
where the dependence of C²22 log(n) is polynomial. We conclude that
(Fx) - (+(x)) - (0 (n -
33

---

2 is almost 1, compared to other terms. Remark 2. A for large d cases, the left behind term (dC2) d-2
The major contribution to the error upper bound is made by the 1 - term. If L is replaced by kL, then n must be at least h to keep the upper bound to remain the same, which implies that the contribution of I is much more significant than
As for the lower bound, Theorem 2 follows immediately from Theorem 3.2 in Györfi et al. 2002 with Hölderness parameter (1, L). Note that class pL we construct has a conditional distribution identical to the Gaussian of the unit variance, which automatically fulfills our Assumption (b) and (c). The only additional requirement in pL. compared to their DD(..)) is that we require the existence for some x E [0, 1] their proof, they construct a sub-class based on the division of [0, 1] into smaller cubes, where the mean function μuxx) (1, L-Hölder (i.e. L-Lipschitz) on
of ((x))
each cube while being zero on the boundary of those cubes. Hence our additional requirement of
Er E [0, 1]d such that p(x)
does not affect the result.
Lemma 14 (Theorem 3.2 in Györfi et al. 2002]). Consider a conditional mean estimation problem (i.e. regression problem) with the estimator sequence { and true conditional mean p. For the class PL, the sequence
24 On La गीス
a lower minimax rate of convergence. In particular,
E Ellpp - Hl]3]
lim inf An (xYY)-PPPEP'
sup
Cs > 0,
L di2
for some constant Cg independent of L.
Proof of Theorem 2 Since for any Gaussian distribution with known variance, estimating its mean is
equivalent to estimating its any 7thh quantile, Theorem 2 is direct corollary from Lemma 14 once
one verifies that pt. satisfies Assumption 1
Verifying Assumption I (a): Since the - conditional quantile is of the type (ux)) + ((T) and p(X) L-Lipschitz, the conditional quantile also of course L-Lipschitz.
Verifying Assumption (b): Because μ(()) is I--Lipschitz and X' bounded, its image (x() is bounded. Furthermore, as we know 0 E (,)t), we can bound ju(t') by [-Lvdd Lvd,, leading to the boundedness of the conditional quantile set (x(² + ((r).
Verifying Assumption 1 (c): The pdf around pu(x) + ((()) is of the same shape as that of standard Gaussian at ((T), which is definitely lower bounded from zero in neighborhood.
I Proofs in Section 3.2
Proof of Theorem 3
By the definition of mutual independence, we have A E F(A), BE BEF(),,
IP(U E BIZ E =(A))P(X E AlZ z(A)) = IP(U B, E A A x
Hence,
E A/Z =(A))
P(U BIZ 2(A))=
P(X E AlZ E x(A))
P(U' E A, 2,
E x(A))
E =(A))
A
P(X A, x(A)) A)
P(
IPU E B),
leading to the conclusion that
FuIx= Fuiz=z(a)-
34

---

Thus, if we have any consistency guarantee mean squared error guarantee on
QUIZ =(x)) (UIZ 2(x))),
then also get the same result for
(UI = =(x)) - Q- =)|-
Notice that the dimension of Z is do & d, we can conclude the proof.
Proof of Theorem By definition of quantiles, we know
P(U & QU())) Z) FUIZ(Q-(UIZ)) 2 T.
Since we assume that Feiz is absolutely continuous with respect to the Lebesgue measure, the Radon-Nikodym derivative exists, which is the conditional probability density function. Thus Fojz must be locally Lipschitz. Then for any fixed bounded neighborhood containing (U)Z), Fojz Lipschitz on I with some Lipschitz constant L(I). Therefore Ve and small enough such that e), C,
Fu)z(Q, 2 Fuiz(Q,(UUZ)) L(T)e.
By the definition of quantiles, we have
Fu]z(Q,(U\Z) - <T.
Hence
Fu]z(Q-(U1Z)) < L(T)c +
Taking have
(u)z(Q+(U1Z)) T. Combining with the fact that (U)2())(U)Z)) 2 r, we prove that
PPU S (-(UIZ) 2) -Fu))(Q-(U1Z))= =
Proof of Theorem 5 arguments.
We only prove the middle inequality, since the remaining two hold from similar
By a similar argument as Proposition one knows ((d2), then the best single point decision they can make with respect to conditional expectation of pinball loss is
(U(l))(da) E arg min (M, U))
WER
Since (dd) E F(Z(4)),, if another learner makes decision based the first d, components of 7(da), say, (d())) then they must suffer a pinball loss smaller than that (dda) learner, which is
(Q- U)) z (II((Q- (
By taking the expectation with respect to (d2)) the tower property of the conditional expectation tells us that
e ((((0- w)
= E (Q+ ),U)) [da]]
E
El1-(-- 2(6)) U) [](2))
((Q,001222 (4)) (.)
E
Noticing that Y - + f(x) and (a (e (= b), we have
4uu+ + f(X), = 1. (u, (u,U) U);
which verifies the proof finally.
Remark 3. One may wonder if one conduct the proof of Theorem 5 in simpler way via Jensen's inequality since 4, (..)) is convex with respect to its first component. But unfortunately, such an argument does not work for the quantile case, since in general,
((Q((U))) * (U).
35

---

J Proofs of Auxiliary Lemmas in Appendix H
Proof of 6 Denote the probability measurement by P. From the definition of pinball loss and expectation, we have
(T 1)( - ujdP + +(U u)dP
00
(T- 1)(U - wjdPP + /
T(U - uddP
For small but positive A > 0. we have
Bu (u U)] - Br
000
(T- - 1)(U - u + A)dP+
T(U u- A]dP
-n-
(T- 1)(U w)dP 1 T(U uddP
1)(U + A)dP
T(U + A)dP
1)(U w)dP ( T(U uddP
1)(U u A)dP+ I" A
T(U a u+ A)dP
(T 1)dP + TdP
(U -u- A)PP
= A(T - F(u + J. Ju--A (U - u- A)dP.
Note that any edf càdlàg (right continuous with left limits). Thus, Ve 0, 36 s.t. AA <8.
IF(u A) F(u Sc.
Then we have
05
(U - u A)dP
a A
AdP
Ac.
Hence
A(T - A(TFF(u)) < Ey 11 - A, - -b < -()
Then
Bu A,
By 14, (u, U)]
T F-F(u) < lim inf A -++
A
Ev [4 (u - A, U)]
Eu (u, U)
< lim sup
A
-+
ST F(u ) +c..
36

---

Since e can be arbitrarily small, the limit exists, and
Ey (il(uu_,U)) - (u,L,(U F(u
a_ LLNFFLL = =- lim A-0+
A
The other side of semi-derivative can be derived in a similar way, which we omit for simplicity.
Lemma 15. For any a1, 02 bi, b2 E I, we have
max (01, b1) max b2)l S max (l01 02). (b b2
Proof of Lemma 15 If a; 2 by (or a, S b,) simultaneously for i 1, 2, the claim is straightforward.
Without loss of generality, we assume d1 b1 b2 02. Then
LHS a2 b1 S 02 a1 RHS.
Proof of Lemma 13
As for part (a), notice that
E
-1- p)-
1+c
Ci+t
1+1 -1--
(n + 1)p
/=0
n+1
$
R+1P (1 - -
(n + 1)p
1=
(p+1 p)) nh
(n + 1)p
(n + 1)p'
where CA2 is the number of combinations of (n,
2 for any h 2 1, we have Since S 1+1
2
E (
115
> 03
<
<
[+
(n + 1)p
For part (b), we prove its symmetric argument that Vr > np,
rg
P(5 2 (r np)
where g - p.
The probability ratio of two adjacent value is
P(C = r +
(n r)p (r + 1)9
P(C = r)
(n r)p
rq
r - n1p
rq
1.
Hence for any 2 r, we have
P(( = k) SP(( = =r) (1-- np
r
k-y
37

---

Summing all these k's from r. to oo, we have
20
P(C 2()2r)= Σ P(c=h)
00
R-r
up
$ P(s r)
1
rq
kmr
P(( r)
r-np
rg
P(C r)
rq
np
From similar argument of inspecting the probability ratio between [np] S k S
know that
P(S k) > PP(C r).
There are at least r [mp'] + +12r np such k's (including r itself), so we have 1 2 (r - pp) - r),
which leads the conclusion that
P(( S
rq
Vr > np.
(r np)
If have r npp, then we repeat the above arguments by replacing with n and reversing p and g, we can prove
(n r)p (mp r)2
P(CC S r) <
Vr np.
38


Reviews:
[{"review_id": "TeoAWVcpNq", "reviewer": "Reviewer_aPS1", "summary": "This paper proposes a new method for quantile regression and for calibrating prediction intervals in regression. The paper first proposes a simple quantile regression method and shows that this method estimates the true quantile curve at the minimax-optimal rate. For calibrating prediction intervals, the basic idea is to 1) decompose the conditional distribution of the response variable into a conditional mean + noise term, 2) estimate the conditional mean, and 3) apply the above quantile regression method to the distribution of residuals (i.e., respose minus estimated conditional mean). The approach is agnostic to the method used to fit the conditional mean itself, and the experiments demonstrate that the proposed method (as well as a supplementary method that also applies dimension reduction) performs well using feed-forward and recurrent neural nets, as well as random forests.", "questions": "1) One of the paper's key ideas is to decompose the prediction interval problem into a mean regression problem and a noise quantile prediction problem. Are there any cases where this decomposition would fail, or would be expected to perform worse than a method that learns the regression quantiles directly?\n\n2) Relatedly, how does the performance of the mean estimator affect the performance of the prediction interval (e.g., if the mean estimator under- or overfits)?", "limitations": "The paper is fairly clear about its limitations, such as the gap between the theoretical performance (which suffers from the curse of dimensionality) and the strong empirical results.", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is quite easy to follow. The proposed method seems both simple and effective, and makes very weak assumptions. Figure 1 makes the advantage of individual calibration very clear, and the theoretical guarantees are hence both useful and impressive in light of existing results on the impossibility of individual calibration. The experiments are also fairly thorough and well-presented.", "weaknesses": "1) Table 1: The \"Ours Best?\" column seems misleading, because it is taking the best of two different methods (NRC and NRC-DR). In several cases, only one of the proposed methods performs best, while the other method performs worse than competitors. However, in practice, one must typically pick one method before knowing which of the two will perform better. So, this column is not really informative of how the proposed methods would perform in practice. I suggest either removing this column. Perhaps a vertical rule could be added to distinguish the current paper's methods from previous methods.\n\n2) The introduction is quite long, and it is not clear to me whether all of this information should be presented so early in the paper.\nFor example, the discussion on Page 2 about \"individual calibration\" is unclear since \"individual calibration\" isn't defined until Page 3. Much of the content could also be refactored into a specific \"Related Work\" section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the feedback and the raised questions. We believe clarifying these questions improves the positioning of our work.\n\nDecomposition/two-step procedure:\n\nAs noted by the reviewer, our paper adopted a two-step approach which first predicts the mean and then calibrates the quantile of the noise. Empirically, we do find that the two-step procedure is better than the one-step procedure that directly predicts the quantile of $Y|X$ (see the attachment for the new experiment). Theoretically, this has a nice justification too. Specifically, from the minimax lower bound in Theorem 2, the problem’s hardness depends on the Lipschitz constant $L$, the sample size $n$, and the feature dimension $d$. The two-step procedure subtracts an initial regression model $\\hat{f}(X)$ from the label $Y$: in the following, we argue why such subtraction can effectively reduce the Lipshitz constant. \n\nGenerally speaking, the conditional expectation function $E[Y|X]$ is highly related and very likely to wax and wane together with the quantile function $Q_{\\tau}[Y|X]$. In this light, subtracting the conditional expectation function $E[Y|X]$ can very likely smooth out the quantile function $Q_{\\tau}[Y|X]$, resulting in a smaller Lipschitz constant. In metaphor, this is quite like the method of controlled variates in Monte Carlo simulation which introduces a controlled random variable that is correlated with the target random variable so as to reduce the variance of the estimator. Here, the conditional expectation function $E[Y|X]$ works as the “controlled variate” for the original “target variate” of $Q_{\\tau}[Y|X]$ (the Lipschitz constant $L$ here corresponding to the variance to be reduced in Monte Carlo simulation). Of course, there might exist data distributions where $E[Y|X]$ and $Q_{\\tau}[Y|X]$ are completely unrelated, and in that case, the two-step procedure doesn’t guarantee an improvement over the one-step procedure. Also, in reality, we will only use the fitted conditional expectation, $\\hat{E}[Y|X]$ (from the mean prediction model). Regardless of overfitting or underfitting, as long as the function $\\hat{E}[Y|X]$ behaves relatedly with the quantile function $Q_{\\tau}[Y|X]$, it can effectively reduce the Lipschitz constant in logic. In the extreme case, the Lipschitz constant L=0 for the error quantile function $Q_{\\tau}[Y-E[Y|X]|X]$ corresponds to the case of homoscedasticity, which is still reasonable for some data distribution; but meanwhile, we can’t really expect the original $Q_{\\tau}[Y|X]$ to have $L=0$, i.e., to be a constant function. \n\nPaper writing and introduction: \n\nWe thank the reviewer for the kind advice. We agree with the comments on the column of “Ours best” and will adjust it accordingly in the next version of our paper. Also, for the introduction, we apologize for the usage of concepts before formal definitions and will revise it. \n\nIn terms of the length of our introduction, it is partly due to the nature of the studied problem that it lies on the intersection of calibration, quantile regression, and conformal prediction, spreading over several fields of machine learning and statistics. The ML literature has a related but different taste than the statistics literature. The calibration literature, which originates from the community of machine learning, is often based on modern machine learning models such as deep neural networks. It usually proposes models at the price of sacrificing some theoretical rigorousness. Concurrently, the statistics literature develops concepts that always have theoretical guarantees. But sometimes it is not so straightforward: some recent advances in conformal prediction claim their excellence in dealing with heteroscedasticity, while still lacking a finite-sample theoretical guarantee on reaching such an individual coverage. For more discussion on this, we refer to our response to Reviewer bxwF. We skip some quantile regression and conformal prediction literature reviews that are less relevant but still result in a long introduction. Moreover, papers from different communities usually don’t cite each other much, though they are working on a similar or even the same problem. In this light, we hope our work can make at least some contribution to improving the awareness and synchronizing the language/understanding for researchers from different communities on this same problem. \n\nWe hope our response addresses the raised questions. If there are any follow-up questions/concerns, we will get back to you timely in the following discussion week. "}, {"review_id": "aFKdVfraUn", "reviewer": "Reviewer_oVNP", "summary": "The paper considers the uncertainty quantification problem for regression models. First, they proposed an algorithm for simple nonparametric quantile estimator. Then, they further proposed the nonparametric regression calibration algorithm. They also provide theoretical analysis of the proposed algorithms and implications.", "questions": "1)\tFor Algorithm 1, do other kernels work for the proposed algorithm? Are there any requirements for the kernels?\n2)\tFor Algorithm 2, could the split proportion be different than half and half? The proportion would influence the results. Are there any experimental results to check the effect of the proportion?\n", "limitations": "NA", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1)\tThe paper is well-organized and well-written. \n2)\tThe paper provides several theoretical results and implications of the theory. \n3)\tThe paper includes extensive experiments. ", "weaknesses": "1)\tFor Algorithm 1, do other kernels work for the proposed algorithm? Are there any requirements for the kernels?\n2)\tFor Algorithm 2, could the split proportion be different than half and half? The proportion would influence the results. Are there any experimental results to check the effect of the proportion?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the questions and comments. \n\nThe choice of the kernel function: \n\nThis simple nonparametric method performs rather robustly with respect to the choice of the kernel function. Essentially, all standard choices of kernels specify a localized weighting regime for error quantile estimation and enjoy similar or even identical theoretical guarantees. Numerically, we presented the naive kernel and the Gaussian kernel (RBF) in the paper but we also tried out the Laplace RBF kernel, normalized triangle kernel, and higher-order kernels on different tasks. And they all give similar performances. We didn’t stress much on the point in the submitted version, but we will include more information about this aspect in the next version of our paper. \n\nData splitting scheme:\n\nWe remark that the split between the mean estimation data and the error quantile estimation data is quite flexible and it doesn’t have to be 50-50. We are sorry for causing the impression of a half-and-half split. The split ratio in fact trades off between learning a good mean prediction model and achieving a decent estimation of quantiles. In the numerical experiments, we take a grid search for ratios between 2:8 to 8:2, and the best proportion is picked based on the performance of the validation set. For some datasets, the best proportion happens to be half-and-half, but this is not always the case. Interestingly, many works on split conformal prediction do recommend such a 50-50 ratio but didn’t provide much theoretical justification. We will update our paper with the point and look forward to see more future investigations on this aspect. \n\nWe hope our response addresses the raised questions. If there are any follow-up questions/concerns, we will get back to you timely in the following discussion week. \n"}, {"review_id": "gkZ20dp5dM", "reviewer": "Reviewer_2q3v", "summary": "This paper studies uncertainty quantification for the regression problem. In particular, \nit considers the estimation of conditional quantiles (of the residuals) via the kernel method. The convergence \nrate of the proposed estimator is established, along with a matching lower bound. The proposed \nmethod is evaluated on multiple datasets and compared with other candidate methods. ", "questions": "1. As mentioned above, I wonder how this work \ncompares with the line of works of (conditional) quantile regression. \n2. It might be helpful to also show in the simulations \nthe results if one directly estimates the conditional quantiles  \nof $Y$.\n3. I wonder if the proposed method can be used within the \nframework of conformal inference and achieve distribution-free \nmarginal calibration as well.", "limitations": " The authors have adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper considers an interesting problem; the examples \nfor showing the unexpected results of existing methods are motivating;\nthe \nsolution provided has solid theoretical properties \nand show satisfactory empirical performance in \nnumerical experiments.", "weaknesses": "I was wondering about the position of this paper in the line of works \nof conditional quantile regression (e.g., Takeuchi et al. (2006); Steinwart et al. (2011)).\nA discussion in this direction will be appreciated.\n\nReferences:\n\nTakeuchi, Ichiro, et al. \"Nonparametric quantile estimation.\" (2006).\n\nSteinwart, Ingo, and Andreas Christmann. \"Estimating conditional quantiles with the help of the pinball loss.\" (2011): 211-225.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our work and raising inspiring questions. \n\nDirect estimation:\n\nWe include more experiments using that directly estimate the conditional quantiles by optimizing pinball loss as the attachment. It does show the advantage of the two-step procedure. We also provide more theoretical explanations for the advantage of this decomposition or two-step procedure in our response to Reviewer aPS1, and we are grateful if you have time and are interested in reading the explanations therein. \n\nMarginal calibration for conformal prediction: \n\nIn this work, we focus mainly on the individual calibration objective; and our main motivation when initiating the project is in fact to provide a positive result for the pessimism of individual calibration mentioned by the existing works. For marginal calibration, there are many existing methods in the literature of conformal prediction that enjoy both successful empirical performance and a nice theoretical guarantee (see also our response to Reviewer bxwF). For our method, if one aims for a marginal calibration guarantee, then a simple way to modify Algorithm 1 would be to remove the kernel weighting and assign a uniform weight to all the samples. The marginal calibration guarantee can then be derived under fewer assumptions (without Assumptions 1 (a) and (c)).\n\nPositioning against (conditional) quantile regression:\n\nFirst, we’d like to thank the reviewer for bringing up these two papers. We believe we have done an exhaustive search over the literature view part of the papers cited in our paper, but we still miss this line of works on conditional quantile regression. As we noted in our paper and our responses to other reviewers, it seems to us that several communities are working on the same problem but aren’t aware of/don’t cite each other. In this light, we have read the papers with great interest and will definitely include them in the next version of our paper.\n\nBoth papers utilize the kernelized method and directly predict the quantile function (one-step). Essentially, methods from both papers search for a quantile prediction function over the RKHS function space $\\mathcal{H}$. Apart from the comparison between the one-step method and our two-step method in our response to review aPS1, we discuss the positioning of our results against these two works in two additional aspects: \n\n- Theoretically, Takeuchi et al. (2006) derived a performance guarantee for the conditional quantile estimator. Their derivation is based on a Rademacher complexity bound over the function class with a bounded $\\mathcal{H}$-norm, and to ensure a consistency result, it requires the true conditional quantile function to have a bounded $\\mathcal{H}$-norm. This is in parallel with our Lipschitz class (Assumption 1 (a)). The Lipschitz function class and the bounded $\\mathcal{H}$-norm function class overlap but don’t contain each other. Steinwart and Christmann (2011) analyze the same algorithm but adopt a different approach. Instead of imposing a bounded $\\mathcal{H}$-norm, they impose a decay rate for the eigenvalues of kernel integral operator, alongside some other assumptions on the underlying distribution (on a minor note, their paper missed a condition in their theorem statement that the chosen kernel should lead to an RKHS dense in $L_1$ space). In short, both papers, together with ours, strive to fight against the impossibility of individual calibration by imposing conditions on the true quantile function so that it belongs to a certain function class. These two works focus on RKHS space while we focus on Lipschitz space; generally, we find the conditions imposed not quite comparable to each other. For the Lipschitz function class, to the best of our knowledge, our work is the first result to provide an individual guarantee. \n\n- Empirically/computationally, our algorithm features an analytical solution and thus is very simple and efficient to implement, while both papers rely on solving a kernelized learning problem that doesn’t scale well with sample size or dimension. For Steinwart and Christmann (2011), it is more of a theoretical work and doesn’t provide any numerical experiments. From a practical viewpoint, our method/framework is more about (i) justifying the widely-adopted two-step calibration procedure and (ii) providing a simple solution and positive result for individual calibration. Though our work and these two papers aim for the same individual guarantee, we don’t quite see our framework as a competitor to them. Because at the end of the day, one can adopt the kernelized method for the mean prediction part (in replacement of a linear model or NN) and/or the error quantile prediction part (in replacement of the simple nonparametric estimator) of our framework. \n\nWe hope our response addresses the raised questions. If there are any follow-up questions/concerns, we will get back to you timely in the following discussion week. "}, {"review_id": "YuwajLY8BB", "reviewer": "Reviewer_bxwF", "summary": "This paper addressed the uncertainty quantification problem in regression models, specifically focusing on individual calibration to characterize prediction model quantiles. To overcome these limitations, they proposed simple nonparametric calibration methods that are both computationally efficient and statistically consistent. Their approach provides insights into individual calibration possibilities and establishes upper and lower bounds for calibration error. They attempted to advance existing theoretical analyses by combining nonparametric and parametric techniques, offering new perspectives on regression calibration regarding the curse of dimensionality and reconciling previous findings on individual calibration impossibility. ", "questions": "I do not have further questions, all the concerns are summarised in the \"weakness\" and \"limitation\" sections.", "limitations": " The proposed method suffers from the curse of dimensionality as I mentioned before and the method relies on very strong assumptions. The authors have not implemented their proposed method on high dimensional datasets. ", "rating": 4, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "The paper is overall well written and the guarantees given are technically sound. The authors have demonstrated the effectiveness of their methods through several experiments. The ideas developed are novel to the best of my knowledge but their practicality is questionable. There are many advancements in the field of conformal inference that give similar guarantees with minimal assumptions. I think the paper would benefit from comparing their developed method in varied settings with the existing approaches in conformal inference to prove its efficacy.", "weaknesses": "In this paper, the authors tackle the uncertainty quantification problem for regression models, focusing on individual calibration. While the proposed nonparametric calibration methods and the accompanying analysis present some interesting ideas, I have serious concerns about the authors' familiarity with previous works in the field, as well as the lack of necessary citations to support their claims.\n\nOne significant issue is the absence of references to previous research on conditional conformal prediction, which is a relevant and well-established framework in uncertainty quantification. The authors should have acknowledged and discussed how their proposed methods relate to or differ from existing conditional conformal prediction approaches. Failure to do so undermines the paper's novelty and raises doubts about the authors' understanding of the current state-of-the-art in the field. There have been various works in this area namely:\n1. \"Conformalized Quantile Regression\" by \"Yaniv Romano, Evan Patterson, Emmanuel J. Candès\",\n2. \"Improving conditional coverage via orthogonal quantile regression\" by \"S Feldman, S Bates, Y Romano\"\n3. \"Class-Conditional Conformal Prediction With Many Classes \" by \"T Ding, AN Angelopoulos, S Bates, MI Jordan, RJ Tibshirani\"\n4. \"Conformal prediction with conditional guarantees\" by \"I. Gibbs, J. J. Cherian, E. J. Candès \"\n5. \"Knowing what you know: valid and validated confidence sets in multiclass and multilabel prediction\" by \"M Cauchois, S Gupta, JC Duchi\"\n\nThe proposed method also suffers from the curse of dimensionality and the authors have only included experiments with several UCI datasets which are relatively low-dimensional. Whereas on the other hand, advancement in the field of conformal inference has produced methods that seamlessly adapt to very high dimensional datasets like ImageNet, MNIST, etc. \n\nFurther, the theoretical guarantees in the paper rely on very strong assumptions, i.e., Assumption 1 of the paper, unlike previous works in this field. \n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the sincere comments and for pointing out our missing literature. In the past week, we spent quite much time working on these papers and related literature. We’d like to summarize our findings below. We hope this better clarifies our positioning, and we look forward to having more discussions with you and other reviewers in the coming week. \n\nPositioning of our work:\n\nThe reviewer raised a question on whether such algorithms (with the goal of individual calibration) are practical for high-dimensional setting. As noted in our paper, our work is more aimed at providing a positive result for individual calibration, given that the existing results on individual calibration for regression problems are largely negative. Let’s separate the regime into low-dimensional and high-dimensional. \n\n- Low-dimensional regime: after reading our work, we guess probably all the reviewers agree that individual calibration is to some extent achievable under the lower-dimensional regime; yet this is a fact that is not acknowledged by the existing literature. Importantly, we also justify the necessity of striving for such an individual calibration by the pinball loss from the downstream task such as inventory management or newsvendor problem (Proposition 2), which serves as a motivating example for many calibration papers. These two combined imply that being intimidated by the hardness of individual calibration is unnecessary and cost-ineffective. \n\n- High-dimensional regime: Our paper answers the question of “when individual calibration is achievable”, but we didn't say that individual calibration is always achievable. For high-dimensional data, we agree with the reviewer that individual calibration can be a goal too ambitious to target. Yet, on the positive end, as in Theorem 3 of our paper (the proof of which is technically very straightforward), if the high-dimensional covariates have low-dimensional statistics that induce some conditional independence, then individual calibration can be achieved with the low-dimensional statistics. Practically, the low-dimensional statistics can be obtained by training a neural network with a small number of neurons in the second-last layer. There has been a line of literature studying the ability of feature extraction for neural network models; Theorem 3 says that such development can be useful in achieving individual calibration in a high-dimensional regime. \n\n- The reviewer mentioned high-dimensional datasets such as MNIST and ImageNet. Both datasets are tested only as a classification task by the existing literature on conformal prediction/calibration, while our paper studies on regression problems. The eight datasets in the numerical experiments are chosen because they are the standard ones used in most (if not all) regression calibration papers. For demonstration purposes, we include a new dataset of Yacht which is of higher dimension (d=100, n=2000), and observe similar advantages of our algorithm compared to the benchmarks. The 100-dimensional vector might not be as large as the MNIST or ImageNet, but when working on an even higher dimension, the second-last layer of the neural network should have a dimensionality of a similar magnitude. \n\nLiterature on conformal prediction:\n\nWe appreciate the reviewer for bringing the papers to our attention. We believe, and as noted by other reviewers, our paper includes a much longer literature list than a normal conference paper; we didn’t miss these papers deliberately (Also, 2 papers mentioned by the reviewer are in fact posted after NeurIPS 2023 submission deadline). Also, as noted in our paper, we found researchers from different communities, conformal prediction, calibration, and quantile regression, working on similar or even identical problem setup but don’t acknowledge much each other’s work. For example, the existing calibration literature doesn’t cite much the papers on conformal prediction, and vice versa. To this end, we hope our work can make at least some minimal effort in syncing the understanding and awareness across communities, and in our paper, we did mention the perspective of conformal prediction throughout the presentation of our results. Hence, we read these conformal prediction papers with great interests and will include a discussion on this line of works in the next version of our paper. We defer a detailed discussion to the Author Rebuttal. \n\nWe also want to bring up one thing to the reviewer. The majority of work in conformal prediction follows a two-step procedure, also known as a split procedure, that first fits a mean prediction model and then calibrates the prediction error. Meanwhile, the majority of work in calibration literature follows a one-step procedure, that usually appends a calibration penalty to the mean prediction loss in the loss function. The two-step split procedure has the apparent advantage of avoiding overfitting and thus avoids too small confidence interval as the one-step procedure. Our paper provides a second justification for this two-step procedure (See our response to aPS1), which we believe complements the existing works on conformal prediction.\n\nLipschitzness assumption:\n\nThe Lipschitzness condition in our paper helps answer one key question of conformal prediction. It should be regarded as the minimum assumption to avoid violating requirements (i) and (iii) (see the Author Rebuttal for the requirements), where the term “minimum” is proved by our Theorem 2 and Theorem 6. Theorem 6 states that without such a Lipschitz condition, any estimator, even if asymptotically consistent, can suffer from an arbitrarily low convergence rate. In other words, it is an unavoidable cost to achieve a meaningful (iii). Theorem 2, on the other hand, proves that the minimax lower bound under a Lipschitz condition is matched (up to some poly-logarithmic factors) by our algorithm. "}]

Author rebuttal:
We thank the reviewers for taking time to read our papers and for all the helpful feedback. We look forward to more follow-up discussions in the coming week. 

We provide an individual response to each reviewer, and we'd like to use the extra space below to discuss the assumptions of our paper in response to the questions raised by Reviewer bxwF on assumptions and comparisons against conformal prediction literature. Sorry for the caused inconvenience and confusion.

The only assumption made in our paper is Assumption 1. Part (a) and part (b) are standard and critical in achieving individual calibration guarantees. Part (c) is used to obtain the finite-sample convergence rate. Essentially, a by-product of our analysis derives a concentration argument for quantile estimators. Part (c) is indeed milder than the existing literature on quantile bandits and nonparametric density estimation that also derives such quantile concentration results. We believe this part of the assumption compared to the literature is of independent interests and a contribution itself. The key is that we adopt a new analysis approach that combines both ideas for parametric and nonparametric analysis. We refer to more discussions in Line 210-219 of our paper. 

Now we proceed with a detailed discussion of the 5 papers mentioned by the reviewer. After careful examination, we don’t believe any of this work (or other works mentioned in our paper’s reference list and our response to the reviewers) provides a comparable result as ours under a comparable or weaker assumption. 

All 5 listed papers are developed in the area of conformal prediction, where the goal is to give a coverage set that contains the true label with a certain probability. Thus the goal of conformal prediction is no harder than our quantile prediction task since correct quantiles induce desired coverage sets. Among these 5 papers, Paper 3 and Paper 5 are for classification problems, and the remaining 3 papers are for regression problems as ours.

First, we rephrase the “impossible triangle” summarized in Paper 4: (i) a conditional/individual coverage (as our Definition 3) (ii) make no assumptions on the underlying distribution (iii) has finite-sample guarantee and asymptotic consistency. Such a triangle is shown by Vovk (2012) and Barber et al. (2021) to be impossible to achieve simultaneously. The listed works as well as ours follow different ways to reconcile the impossibility:

Paper 1 and Paper 5 focus on marginal coverage (as we define in our Definition 1), which greatly relaxes (i). Paper 3 and Paper 4 consider some milder violations of (i) but in different ways. Paper 3 considers a “clustered conditional coverage”, which is parallel to the group calibration defined in our Definition 2. Paper 4 considers another midway between marginal coverage and conditional coverage. Yet, these two papers are still far from conditional coverage (i). 

In short, all these 4 papers have relaxed requirement (i). But what’s the cost to reach such individual coverage theoretically and what’s the price of keeping the requirement (i)? These questions are not answered by these papers. Our paper answers them with minimal relaxations on the other two: we only relax the “no distribution assumption” requirement (ii) with our Assumption 1. In comparison, Paper 2 also doesn't relax requirement (i) and is the most related work to ours. However, it relaxes both (ii) and (iii), and it violates (ii) to a further extent than ours. To see it, the algorithm in Paper 2 is based on a regularized pinball loss/ interval score loss, where the regularizer forces the covering interval length to be independent of the coverage indicator. It achieves a theoretical guarantee that the true conditional quantile is the minimizer of the regularized population loss under (1) a realizability assumption (that is, the true conditional quantiles lie in the function class considered) and (2) an assumption that the conditional distribution of $Y$ on $X=x$ is continuous with respect to $x$. 

Apart from stronger assumptions, Paper 2's theoretical results have several limitations. First, it achieves only asymptotic consistency compared to our finite-sample guarantee in Theorem 1, which means that Paper 2 violates requirement (iii). Besides, it only achieves a necessary condition that the true conditional quantile is the minimizer of the population loss but not vice versa. Second, Paper 2 violates (ii) to a much deeper content than ours. Note that the key assumptions made in Paper 2 are the zero-approximation-error assumption and the continuity assumption. The former is not a weak one in the statistical learning theory, meaning that there should be no model misspecification at all, which limits the guarantee from a desirable “distribution-free” property. On the contrary, our algorithm can have a finite-sample guarantee with only three very general assumptions. Even for the continuity part of which the reviewer is unsatisfied, our assumption can still be relaxed to a weaker one than Paper 2 to reach the same asymptotic result. To see that, our result only requires the conditional quantile to be continuous, while Paper 2 requires continuity for the whole conditional distribution. Third, the regularization term in Paper 2 is related to a zero-one indicator, which is discontinuous for gradient descent. The authors twist their algorithm by replacing it with a smooth approximation to overcome the issue, without any formal guarantee on the twisted one. In contrast, our theory accompanies the exact same algorithm that we propose. 

In short, many existing works on calibration completely ignore works on conformal prediction. We didn’t do this in our paper; in contrast, we’d like to call for more mutual awareness between these communities. Also, while we position our work along with the calibration literature, our theoretical analysis complements the existing literature on conformal prediction.

Meta-review: The paper has received mostly favorable reviews, with the exception of reviewer bxwF, who has noted the paper is missing a comparison to work in conformal prediction, and expressed concern over the proposed method appearing to suffer from the curse of dimensionality. The authors have provided extensive responses addressing these two issues. After additional deliberation and discussion with the reviewers, I believe that - once properly extended and incorporated in the paper - the author's responses sufficiently address these issues. Consequently, I recommend acceptance subject to the following changes in the camera-ready paper:

* The literature review must include detailed discussion of all pertinent conformal prediction literature. Following commentary by reviewer aPS1, I encourage the authors to streamline the introduction by separating the detailed discussion of related work into its own section or subsection.

* The experiments must include the conformal prediction discussed during the review process as baselines.

* The experiments must report prediction interval sizes among other metrics, and discuss the pros and cons of various metrics along the lines mentioned during the review process.

* The authors should include additional discussion on the prospects of applying the proposed method to high-dimensional models such as neural networks, and perhaps add experiments with applying the model on second-to-last layer embedding from pre-trained neural networks.
Decision: Accept (poster)