Conformal Prediction using Conditional Histograms
Authors: Matteo Sesia
Yaniv Romano
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: Conformal prediction
quantile regression
histograms
skewed data.

Abstract: This paper develops a conformal method to compute prediction intervals for non-parametric regression that can automatically adapt to skewed data. Leveraging black-box machine learning algorithms to estimate the conditional distribution of the outcome using histograms, it translates their output into the shortest prediction intervals with approximate conditional coverage. The resulting prediction intervals provably have marginal coverage in finite samples, while asymptotically achieving conditional coverage and optimal length if the black-box model is consistent. Numerical experiments with simulated and real data demonstrate improved performance compared to state-of-the-art alternatives, including conformalized quantile regression and other distributional conformal prediction approaches.

Conformal Prediction using Conditional Histograms
Matteo Sesia
Department of Data Sciences and Operations University of Southern California, USA sesi arssi@@masshall usc.edu
Yaniy Romano
Departments of Electrical and Computer Engineering and of Computer Science Technion, Israel yromano@technion..acc ac.
Abstract
This paper develops conformal method to compute prediction intervals for non- parametric regression that automatically adapt to skewed data. Leveraging black-box machine learning algorithms to estimate the conditional distribution the outcome using histograms, it translates their output into the shortest prediction intervals with approximate conditional coverage. The resulting prediction intervals provably have marginal coverage in finite samples, while asymptotically achiev- ing conditional coverage and optimal length if the black-box model is consistent. Numerical experiments with simulated and real data demonstrate improved perfor- mance compared to state-of-ther oh-trr-rrtttttt alternatives, including conformalized quantile regression and other distributional conformal prediction approaches.
Introduction
Problem statement and motivation
We consider the problem of predicting with confidence a response variable Y E R. given features X E IRP for a test point n++1, utilizing pairs of observations ((X(), y()) }- drawn exchangeably (e.g., i.i.d.) from some unknown distribution, and leveraging any machine-learning algorithm. Precisely, Va (0, 1), we seek a prediction interval (x,o((x+1)) following three criteria. First, Ön,a should have finite-sample marginal coverage at level 1 a,
for Yo+1 satisfying the
[Yn+1 Cn, (Xn+))
2 11-a
(1)
Second, On should approximately have conditional coverage at level a,
Yn+1 Cn.a (x) X = 2 1 a-a Vr E IR,,
(2)
meaning it should approximate this objective in practice, and ideally achieve it asymptotically under suitable conditions the limit of large sample sizes. Third, Ćn,a should be as narrow as possible.
We tackle this challenge with conformal inference 24 35 which allows one to convert the output of any black-box machine learning algorithm into prediction intervals with provable marginal cov- erage (1). The key idea of this framework is to compute conformity score for each observation, measuring the discrepancy, according to some metric, between the true value of and that predicted
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

by the black-box model. The model fitted on the training data is then applied to hold-out calibration samples, producing a collection of conformity scores. As all data points are exchangeable, the empirical distribution of the calibration scores can be leveraged to make predictive inferences about the conformity score of new test point. Finally, inverting the function defining the conformity scores yields a prediction set for the test Y. This framework can accommodate almost any choice of conformity scores, and in fact many different implementations have already been proposed to address our problem. However, it remains unclear how to implement concrete method from this broad family that can lead to the most informative possible prediction intervals. Our contribution here is to develop practical solution, following the three criteria defined above, that performs better compared to existing alternatives and is asymptotically optimal under certain assumptions.
It is worth emphasizing that constructing short prediction interval with guaranteed coverage is a reasonable approach to quantify and communicate predictive uncertainty in regression problems, although it is of course not the only one. To name an alternative, one could compute a non-convex prediction set with analogous coverage [20). which might be more appropriate some situations, but is also more easily confusing. For example, it could be informative for physician to know that the future blood pressure of a patient with certain characteristics is predicted to be within the range [120,129] mmHg. However, it would not be more helpful to report instead the following non-convex region: [120, 120.012] U [120.015, 120.05 U [121, 122.7 7] U [123.1, 127.2] U [127.8, 129] mmHg. Indeed, in the second case it would not be clear (a) whether the multi modal nature of that prediction is significant or a spurious consequence of overfitting, and (b) how the physician would act upon that prediction any differently than if it had been [120,129] mmHg. Therefore, we focus on prediction intervals in this paper because they are generally easier to interpret than arbitrary regions, and they are also less likely to convey false sense of confidence.
1.2 Preview of conformal histogram regression
Imagine an oracle with access to Prix, the distribution of Y conditional on X, which leverages such information to construct optimal prediction intervals as follows. For simplicity, suppose Pr|x has a continuous density fly xx with respect to the Lebesgue measure, although this could be relaxed with more involved notation. Then, the oracle interval for Yn+1 X = would be:
Coracle(x)) = [1omacl (x),, uprracee (x)]
(3)
where, for any (0, 1]. poracle (x) and u gracle (x:) are defined as:
{ {u -u-: 1." fly [[ylx)dy= 2
[oraall (z), u (eracle(r)) := arg min (aa)RRR.1S
(4)
This is the shortest interval with conditional coverage 2 the solution to is not unique (e.g.. if f(- iz) is piece-wise constant), the oracle picks any solution at random. Of course, this is not practical method because f is unknown. Therefore, we will approximate 4 by fitting a black-box model on the training data, and then use conformal prediction to construct interval accounting for any possible estimation errors. Specifically, we replace f in (4) with a histogram approximation, hence why we call our method conformal histogram regression, CHR. The output interval is then
( [x)),, (+(())
(5)
where (x) and 0pp (I) approximate the analogous oracle quantitics in 4.. The value of f in (5) will be determined by suitable conformity scores evaluated on the hold-out data, and it may be larger than
1 - a if the model for f is not very accurate. However, if the fitted histogram is close to the true
Pylx, the interval in 55) will resemble that of the oracle
Figure previews an application to toy data, comparing CHR to conformalized quantile regression (CQR) [30] see Section 4.2 for more details. CHR finds the shortest interval such the corresponding area under the histogram is above T, for any 7 E (0, 1]. and then calibrates 7 to guarantee marginal coverage above - ox; this extracts more information from the model compared to CQR. For example, CHR adapts automatically to the skewness of Y X, returning intervals delimited by the 0%-90% quantiles in this example, which are shorter than the symmetric ones (5%-95% sought by CQR.
2

---

20
0.20
CHR lower
lower
CHR upper
COR upper
15
0.15
10
(a)
(b)
0.05
0.00
0.4
Figure 1: CHR prediction intervals in example with variable, compared those obtained with CQR [30] Both methods guarantee 90% marginal coverage and are based on the same deep quantile model. (a) Histogram estimate of Pylx for point with X z 0.2. The CHR interval corresponds to the shaded part of the histogram, whose area is approximately 0.9, as marked by the solid vertical lines. The dashed lines denote the CQR interval. (b) Prediction bands for the two methods, function of X. CHR: empirical marginal coverage 0.9, estimated conditional coverage 0.9, and average length 3.2. The corresponding quantities for CQR are: 0.9, 0.9, and 5.2, respectively.
1.3
Related work
This work is inspired by the conformity scores introduced by [1 for multi-class classification, the underlying idea of which can be repurposed here. Nonetheless, the extension to our problem involves several innovations. This paper connects [31] other conformal methods for continuous
responses 24
which sought objectives similar to ours by leveraging quantile regression [17 non-parametric density estimation [TO] 21,, sometimes considering multi-modal prediction sets instead of intervals [20]. Our approach also exploits black-box models for distributional
25
35
23
30
32
36
estimation; however, we introduce more efficient conformity scores.
We seek the shortest intervals with marginal coverage while approximating as well as possible conditional coverage, although the latter is impossible guarantee in finite samples [16 34]. The performances of prior approaches have been measured terms of these criteria, yet others have not sought them directly. Indeed, if the black-box model is consistent for Pyjx. our method becomes asymptotically equivalent to the oracle - under some technical assumptions. This property does not hold for other existing methods because they tend produce symmetric intervals, with fixed lower and upper miscoverage rates (the probabilities of the outcome being either below or above the output interval, respectively), which may be sub-optimal if the data have unknown skewness.
2 Methods
The proposed method consists of four main components: the estimation and binning of conditional model for the outcome, the construction of a nested sequence approximate oracle intervals based on the above, the computation of suitable conformity scores, and their conformal calibration.
2.1 Estimating conditional histograms
We partition the domain of Y into bins by), for some sequence by <... bm- With little loss of generality, assume Y is bounded: -C by < Y < bn, = C, for some C' 0. Then, we solve a discrete version of the problem stated in the introduction: we seek the smallest possible contiguous subset of bins with 1 - a predictive coverage. If m is large and the bins are narrow, this problem is not very different from the original one, although it is more amenable solution.
For simplicity, we present our method from split-conformal perspective [4] 35): extensions to other hold-out approaches 7 22 35 will be intuitive. Let trrain Deal C (1...., n) denote any partition of the data into training and calibration subsets, respectively. Dirain is used to train a black-box model for the conditional probabilities that Y is within any of the above bins: Vj (1, (1....,m).
(z) = P|Y E [b,-1,b,) X x².
6

---

There exist many tools to approximate Prix and obtain estimates (z) of Ty (x), including quantile regression 27 28 331, Bayesian additive regression trees [T]. or any other non corapaamtttric condi- tional density estimator [14] 19 26 Our method can directly be applied with any of these models, but we found multiple quantile regression to work particularly well 17 23
30 [2], and therefore we will focus on it in this paper. Referring to Supplementary Section S1.1 for implementation details and information about the computational cost of the learning algorithm (which is comparable to that required by CQR [30], we thus take these black-box estimates ft,22) (2) as fixed henceforth.
Note that estimating conditional distributions is more challenging the number of variables is larger. However, this is a fundamental difficulty of high- nigh-dimensional regression, not a particular limitation of the proposed CHR. Although our method utilizes conditional histograms learnt from the data, its performance is not directly measured in terms of how closely these resemble the true Prix Instead, as we shall see, CHR only needs to detect the possible skewness of Y X and estimate reasonably well some lower and upper quantiles of this conditional distribution. Therefore, its estimation task is not much more difficult than that of CQR [30], as skewness is relatively easy to detect.
2.2 Constructing a nested sequence of approximate oracle intervals
For any partition B- (bn..... ... b......,........... ) of the domain of Y. let (TI..... Tun) ) be unit-sum sequence, depending on x E IRP; this may be seen as a histogram approximation of Py (6. For simplicity, assume all histogram bins have equal width, although this is unnecessary. Then, define the following
.,))
bi-valued function S taking as input x E IR,, n, T E (0, 1). and two intervals s S (1,...,m)) m].
u
st,
arg min (L.)) (I....))) Piga
-u-1 11 Σ (T)()) S E ( S+ j=l
(7)
Above, it is implicitly understood we choose the value of (l, minimizing Σ)= among the feasible ones with minimal Iu 1- 2], if the optimal solution would not otherwise be unique. Therefore, we can assume without loss of generality the solution to (7) is unique; if that is not the case, we can break the ties at random by adding a little noise to I. The problem in can be solved at computational cost linear in the number of bins, and it is equivalent to the standard programming challenge of finding the smallest positive subarray whose sum is above a given threshold. Note that we will sometimes refer to intervals on the grid determined by B as either contiguous subsets of ..] (c.g.. 5))
) or as pairs of lower and upper endpoints (c.g., u)).
If S- = 0 and S+ = {l...., m.. the expression in 7 computes the shortest possible interval with total mass above according (zx) Further, if π. is the mass in the j-th bin according to the true Py then S(x, n, 0, f1..... m). 1 - a) is the discretized version of the oracle interval ( 4 In general, the optimization in 7 involves the additional nesting constraint that the output S must satisfy S S S C St, which will be needed to guarantee our method has valid marginal coverage (77. Intuitively,
helpful to work with a nested sequence because this ensures the prediction intervals are monotone increasing in T, essentially reducing the calibration problem to that of selecting the appropriate value of 7 that yields the desired marginal coverage. Note that the inequality in 77) involving may not be binding at the optimal solution due to the discrete nature of the optimization problem. However, the above oracle can be casily modified by introducing some suitable randomization in order to obtain valid prediction intervals that are even tighter on average, as explained in Supplementary Section $1.2
As á may be an inaccurate estimate of Py|x, we cannot simply plug it into the oracle in and expect valid coverage. However, for any approximate conditional histogram ât, we can define a nested sequence [17] of (randomized) sub-intervals of B, for different values of ranging from 0 to 1. Then, we calibrate 7 to obtain the desired
a marginal coverage. Precisely, consider an increasing scalar sequence 7. = t/T, for 1. E (0.... .., T} with some T E N, and define a corresponding growing sequence of subsets S, C (1..... m) as follows. First, fix any starting point i E (0...., . T) and define Sj by applying (7) without the nesting constraints (with S = @ and S+ {1, (I....,,... m)):
m), S = S(a, n, 0, (1...,,,,, Té);
Note the explicit dependence on x and π of the left-hand-side above is omitted for simplicity, although it is important to keep in mind that S; does of course depend on these quantities. Figure 2 (second row) visualizes the construction of S in toy example with 7 0.9.

---

(a) Estimated conditional histogram
(b) Nested prediction intervals
Ti+1 = 0.97
0.10
+ (wider)
0.05
0.00
0.10
0.9
S (initial)
0.05
0.00
0.10
T 1
0.8
S 1 (narrower)
0.05
0.00
-2
Figure 2: Schematics for the construction of a nested sequence of approximate oracle prediction intervals (8) TTOp. (a) Conditional histogram approximation of the distribution of Y X, based black-box model. The shaded areas delimited by the dashed vertical lines denote the shortest intervals with the desired mass (r) under the histograms, subject the nesting constraints. (b) Sequence of prediction intervals. The initial interval Sj is not subject to any nesting constraints. The wider (above), or narrower (below), intervals must contain St (above), or be contained in it (below).
Having computed the initial interval S. for t f, we recursively extend the definition to the wider intervals indexed by t = f- 1,..., T as follows:
S == S(x, n, St-1, {l...., m)
See the top row of Figure 2 for schematic of this step. Similarly, the narrower intervals S, indexed
by t = f - 1, 1, .....0 0 are defined recursively as:
S = S(x, π, 0, TT))
(10)
See the bottom row of Figure for schematic of this step. As a result of this construction, the sequence of intervals 1T is nested regardless of the starting point in ( as previewed in Figure? 2 However, different choices of may lead different sequences, any of which allows us to obtain provable marginal coverage, as discussed next. As our goal is to approximate the oracle in
accurately, the most intuitive choice is pick E such that Ti 11-60 a. A more involved randomized version of this construction, inspired by the more powerful randomized oracle, discussed in Supplementary Section [11.2] Note that the randomized version of the nested prediction intervals will be the one applied throughout this paper and, with a slight overload of notation, we will simply refer to it as {S Note also that we will highlight the dependence of this sequence on and i by writing it as S, Further, as we work henceforth with the randomized versions of these prediction intervals (described in Supplementary Section (1.2). we will refer them as St e, T)) where & is a uniform random variable in [0. 1]. independent of everything else.
2.3 Computing conformity scores and calibrating prediction intervals
Given any sequence of nested sets (, s, R), we define the following conformity score function E:
E(x, y. E, x := min {t E fo...., T) E S,(r,e, m)}-
(11)
In words, this computes the smallest index such that S, E, 7) contains y. as in 117 31I Equiva-
lently, one can think of these scores as indicating the smallest value of the nominal coverage r in (
necessary to ensure the observed Y is contained in the prediction interval. Our method evaluates [I]
on all calibration samples (X,, Y:) using the क learnt on the training data; for each i E Deal we generate E, ~ (ni 1) and store
E, E(X, Y,,,,
5

---

Then, we compute prediction intervals for Yx+1 by looking at the nested sequence in 8 TO corresponding the new X and selecting the interval indexed by the a quantile (roughly) of (E)).. The procedure is outlined in Algorithm I Note that the only computationally expensive component of this method the estimation of the conditional histograms (see Supplementary Section [T.J for details); the construction of the nested prediction intervals and the evaluation of the conformity scores have negligible cost because the optimization problem in is an easy one.
It may be helpful to point out that, if т provides an accurate representation of the true conditional distribution of Y X, then the above conformity scores are uniformly distributed [3T] In that ideal case, no calibration is needed and indeed our method simply reduces to applying to construct prediction intervals with 90% coverage. In practice, however, т can only be a possibly inaccurate estimate of Pyix (hence why we will refer to it as म from now on), which means that the distribution of the conformity scores may not be uniform and the conformal calibration is necessary to obtain valid coverage.
with 7 = 0.9
The next result states that the output of our method has valid marginal coverage, regardless of the accuracy of न. The proof relies on the sequence S, being nested; from there, coverage follows from the results of 17 31; see Supplementary Section ST.I
Algorithm 1: CHR with split-conformal calibration
Input: data {(x, Y.) Xn+1, partition B the domain of into m equal-sized bins, level a E (0, 1), resolution T' for the conformity scores, starting index E for recursive definition of conformity scores, black-box algorithm for estimating conditional distributions. Randomly split the training data into two subsets, Dtrain Deal
3 Sample Eạ ~ Uniform (0, 1) for all i E {1...., n + }}. independently of everything else.
4 Using the data in Dirain, train any estimate of the mass of Y X for each bin in B (6) see Supplementary Section [ for concrete approach based quantile regression.
5 Compute E, = EX,,, Y,, Ei, T) for each i E Deal, with the function E defined in [1
Compute i ((F) Dea as the [(1 a1 (Deml) ) ][th smallest value in (E4)
Select the f-th element from (S,(X En+1, a) 1T defined in
10p:
On Cn (X,+1)))))))))) S;
Output: A prediction interval Crc (X,,((Xn+1) (Xn+1) for Yn+1.
1}, are exchangeable, then the Theorem 1 (Marginal coverage). If (X fir, E {1,...,,, n.n + oufput of Algorithm satisfies:
(YYtI One (Xn+1) 2 -
(12)
Note that Theorem I provides only lower bound; nearly matching upper bound on the marginal coverage can be generally established for split-cou colffronformll inference if the conformity scores are almost-surely distinct [24 30 35 Although the CHR scores (TT) are discrete, our experiments will show the coverage is tight as long as the resolution T is not too small.
3 Asymptotic analysis
We prove here that the prediction intervals computed by CHR (Algorithm 1 are asymptotically equivalent, as + 20, to those of the oracle from 3 44 if the model î is consistent for Pryx and few other technical conditions are met. In particular, we analyze a slightly modified version of Algorithm in which there is no randomization; this is theoretically more amenable and equivalent in spirit, although it may yield wider intervals in finite samples. Our theory relies on the additional Assumptions 5 explained below and stated formally in Supplementary Section S2
1. The samples are i.i.d., which stronger than exchangeability; this the key our concen- tration results.
2. The block-box model estimates Prix consistently, sense analogous to that in This assumption is crucial and may be practically difficult to validate in practice, but

---

can be justified by existing consistency results available for some models under suitable conditions. such as random forests [27]. Further, the resolution m of the partition of the Y domain should grow with at a certain rate, and the resolution T of the scores E4 in (11) should grow as Tn = n.
3. The true Pylx continuous and with bounded density within finite domain. This assumption is technical and could be relaxed with more work.
4. The true Py|x is unimodal. This assumption also technical and could be relaxed with more work.
5. The estimated conditional histogram से preserves the boundedness and unimodality of Pyix; this assumption may be unnecessary but it is convenient and quite innocuous at this point given Assumptions 4
For simplicity, we assume the number of observations is 2n, the test point is (X2n+1; Y2a+ 1), ) and Dirain Deal although different relative sample sizes would yield the same results.
Theorem 2 (Asymptotic conditional coverage and optimality). Va E (0, 1], let One Ca (X2++1) denote the prediction interval for Yzu+1 computed by Algorithm arr level a without randomization.
Under Assumptions 1 Ose (X2++1) is asymptotically equivalent, as + x., to Coracle (X2++)) 1), particular, the following two properties hold.
the output of the orucle
(i) Asymptotic oracle length. For some sequences Yn. 00 (0add En. - +0 - 00.
IP O (X2++1)1 S |C (X2n+1)l + + (Carall((22n1))) 2 1 5m.
(ii) Asymptotic conditional coverage. For some sequences En + 0 and Sn + 0 as + oc,
P P Y E One (X2411))
X2n+1
2 a en 2 Ca.
Theorem is similar to results in 24 and 32 about the efficiency of earlier approaches to conformal
regression, including CQR [32 However, the increased flexibility of our method is reflected by the oracle in Theorem 2 which is stronger than those in [24] 32). In fact, the oracle in (24) does not have conditional coverage, and that [32 produces wider prediction intervals with constant lower and upper miscoverage rates. Other conformal methods based on parammaettric density are not as efficient as CHR, in the sense that Theorem 2 does not hold for them.
estimation 10
4
Numerical experiments
4.1 Software implementation
along with code to reproduce the following numerical experiments. This software divides the domain of Y into desired number of bins with equal sizes, depending on the range of values observed in the training data; we use 100 bins for the synthetic data and 1000 for the real data. Then, we estimate the conditional histograms using different black-box quantile regression models [27 33,, with grid of quantiles ranging from 1% to 99%; see Supplementary Section SI.I Our software also supports Bayesian additive regression trees [T] and could easily accommodate other alternatives. For simplicity, we apply CHR and other benchmark methods by assigning equal numbers of samples to the training and calibration sets; this ensures all comparisons are fair, although different options may lead to even shorter intervals (32) See [8] for rigorous discussion of how this choice affects the variability of the coverage conditional on the calibration data, which is an issue we do not explore here. See Supplementary Section [S3 for details about how the models are trained, and information about the necessary computational resources.
Python implementation of CHR is available online at athttps://gitthmm http:p/:///gittt...oo////s/s///e
4.2 Synthetic data
We simulate synthetic data set with a one-dimensional feature X and continuous response Y, from the same distribution previewed Figure T which similar to that utilized in 30 to present CQR. Our method is applied to 2000 independent observations from this distribution, using the first 1000 of
7

---

them for training a deep quantile regression model, and the remaining ones for calibration. Figure visualizes the resulting prediction bands for independent test data, comparing them the analogous quantities output by CQR. Both methods are based on the same neural network and guarantee 90% marginal coverage, but ours leads to narrower intervals. Indeed, the advantage of CHR is that it can extract information from all conditional quantiles estimated by the base model and then automatically adapt to the estimated data distribution. By contrast, CQR [30] can only leverage a pre-specified lower and upper quantile (e.g., 5% and 95% in this example), and is therefore not adaptive to skewness.
Figure 3 (a) summarizes the performance of CHR over 100 experiments based on independent data sets, as function of the sample size. We evaluate the marginal coverage, approximate the worst-slab conditional coverage [9] as in [31], and compute the average interval width. We consider two benchmarks in addition to CQR [3] distributional conformal prediction (DCP) [IO) and DistSplit (21). To facilitate the comparisons, all methods have the same base model. (We also applied DistSplit as implemented by [21] with different base model, but the version presented here performs better.) These results show CHR leads to the shortest prediction intervals, while simuitaneously achieving the highest conditional coverage. Compatibly with Theorem 2 the output of CHR becomes roughly equivalent to that the omniscient oracle as the sample size grows; the latter can be implemented exactly here because we know the true data generating process.
Supplementary Figure compares the performance of CHR these experiments to that of naive uncalibrated 90% prediction intervals based on the same deep neural network regression model, and obtained by simply plugging into the oracle prediction intervals do not generally have the desired marginal coverage; in this case, they tend to be too narrow if the sample size is small and too wide if the sample size large. Although the lack of coverage for the uncalibrated intervals is not very pronounced here because this black-box model relatively accurate even with n = 100, such naive approach can yield arbitrarily low coverage in general, especially if the learning task is more difficult (e.g., for high-dimensional X), and is thus not reliable.
with T 0.9. Unsurprisingy, the naive
Marginal
Conditional
Width
0.950
Method
Coverage
0.925
CHR
Width
(a)
0.900
5
CQR
0.875
DCP
DistSplit
0.850
100
300
1000
3000
100
300
1000
3000
100
300
1000
Sample size
Sample size
Marginal
Conditional
Width
0.950
Method
Coverage
0.925
CHR
Wiath
5
(b)
0.900
CQR
0.875
4-
DCP
DistSplit
0.850
31
0
2
Skewness
Skewness
Figure 3: Performance of our method (CHR) and benchmarks on synthetic data distributed as in Figure The dashed lines and curves correspond an omniscient oracle. The vertical error bars span two standard errors from the mean. (a) Performance vs. sample size. (b) Performance vs. average skewness of the conditional distribution of the response, with sample size of 1000. The maximum skewness (near 3) matches that of the data in (a).
Figure (b) shows analogous results from experiments in which we fix the sample size to 1000 and vary instead the skewness of the data distribution. Precisely, we flip a biased coin for each data point and transform Y into -Y if it lands heads, varying the coin bias as control parameter. At one end of this spectrum, we recover the same skewed data distribution as in Figure S (a); at the other end, we have symmetric Py x. Our results are reported function of the expected skewness, defined
8

---

as EY 1x(x))% where j(X) and ((X) are the mean and standard deviation of Y X, respectively. These experiments show all methods are equivalent terms of interval length if Pyxx is symmetric (skewness close to 0), while CHR can be much more powerful if Pyjx x is skewed.
4.3 Real data
We apply CHR to the following seven public-domain data sets also considered in 30 physicochem-
ical properties of protein tertiary structure (bio) [6,, blog feedback (blog) [1] Facebook comment volume [2], variants one (fbl) and two (fb2), from the UCI Machine Learning Repository [13]; and medical expenditure panel survey number 19 (meps19) [3). number 20 (meps20) [4]. and number 21 (meps21) [5). from [3.. We refer to [30 for more details about these data. As in the previous section, we would like to compare CHR to CQR, DistSplit, and DCP. However, as DCP 10] is unstable on all but one of these data sets, sometimes leading very wide intervals, replace it instead with hybrid benchmark that we call DCP-CQR. This improves the stability of DCP by combining it with CQR [30], as explained in Supplementary Section S1.3] This limitation of DCP may be explained by noting the method needs to learn a reasonably accurate approximation of the full conditional distribution of Y X, and its performance particularly sensitive the estimation of the tails, which is most difficult; see Supplementary Section ST.3 for more details. By contrast, CHR is robust because only needs estimate a histogram with relatively few bins a much easier statistical task and then it specifically focuses on finding the shortest intervals containing high probability mass. We apply all methods. including our CHR, based on the same deep quantile regression model. Their performances are evaluated as in the previous section, averaging over 100 independent experiments per data set. each experiment, 2000 samples are used for training, 2000 for calibration, and the remaining ones for training. All features are standardized to have zero mean and unit variance. The nominal coverage rate is 90%.
bio
blog
fb1
fb2
meps19
meps20
meps21
Conditional coverage
1.00
0.95
0.90
0.85
0.80
0.75
bio
blog
fb1
fb2
meps 19
meps20
meps21
Wiath (relative)
2.5
2.0
1.5
1.0
Method
Method
CHR
CQR
DCP-CQR
DistSplit
Figure 4: Performance of our method (CHR) and benchmarks on several real data sets, using a deep neural network model. All methods provably have 90% marginal coverage. The box plots show the distribution of results over 100 random test sets, each containing 2000 observations.
Figure 4 shows the distribution of the conditional coverage and interval width corresponding to different methods, separately for each data set. To simplify the plots by using shared vertical axis, the widths of the prediction intervals are scaled, separately for each data set, so that the smallest one is always equal to one. Marginal coverage is omitted here because all methods provably control it; however, can found in Supplementary Table ST All methods perform well in terms of worst-slab conditional coverage. CHR outperforms the others in terms of statistical efficiency, it consistently leads to the shortest intervals. CQR and DistSplit are comparable to each other, while DCP-CQR sometimes outputs wider intervals. Supplementary Figure S2 shows that analogous results

---

are obtained if we utilized a random forest model instead of a deep neural network. Supplementary Table STT summarizes these results in more detail, including marginal coverage and the omitted performance of the original DCP. Finally, Supplementary Figure S3 compares the performance of CHR in these experiments with real data to that of naive uncalibrated 90% prediction intervals based on the same deep neural network regression model, as in Figure [ST These results show that the naive prediction intervals do not generally have the desired marginal coverage; in some cases they are too narrow, and in others they are too wide.
5 Conclusions
This paper developed CHR, a non non-parametric regression method based novel conformity scores leading to shorter prediction intervals with coverage, and enjoying stronger asymptotic efficiency, compared to the state-of-the-ar alternatives. Of course, real data are finite and our theory relies on assumptions which may be difficult to validate; nonetheless, it is a sanity check and provides an informative comparison. Further, the experiments confirm CHR performs well in practice.
The ability of CHR to automatically adapt to unknown skewness may prove useful practice, as empirical data often follow distributions with power-law tails [12). Indeed, the data sets analyzed in Section 4.3 tend to have highly skewed outcomes with many observations equal to zero. At the same time, limitation of CHR is that it does not rigorously control the lower and upper miscoverage rates, which may be important for some applications; if that is the case, the modified version of CQR proposed by [30] would be better choice. Note that the standard implementations of CQR and of the other benchmarks TO 21 considered this paper are not guaranteed to separately control the lower and upper miscoverage rates. In any case, users of our method could naturally obtain approximations of the lower and upper miscoverage rates for any prediction interval by looking at the underlying conditional histograms, although these estimates of course not calibrated in finite samples.
Algorithm ST in Supplementary Section extends CHR to accommodate cross-van vross-validation+ which is often more powerful, and computationally expensive, compared to the split-conformal approach presented in this paper. The strategy is the same as that followed by [3T] in the classification setting, although it requires extra step, in which the standard cross- validation+ prediction set 77 is replaced by its convex hull to guarantee the final output is an interval [17]. Supplementary Theorem establishes that Algorithm leads to marginal coverage above 2a, applying the more general theory from [7]. Finally, a possible directions for future research may involve the extension of our method to deal with multi-dimensional responses Y.
From a broader perspective, this paper falls within a rapidly growing body of works focusing on improving the interpretability and statistical reliability of machine learning algorithms. Prediction intervals with marginal coverage provide convenient way of communicating uncertainty about the accuracy of any machine learning black-box, which is important to increase their reliability, to ensure their fairness (29), and to facilitate their acceptance. Furthermore, conformal prediction intervals provide principled metric by which to compare different machine learning algorithms [18].
Acknowledgments and Disclosure of Funding
The authors are grateful Stephen Bates, Emmanuel Candès, and Wenguang Sun for providing insightful comments about an carlier version of this mamuscript. M.S. thanks the center for Advanced Research Computing at the University of Southern California for providing computing resources. Y.R. was supported by the Israel Science Foundation (grant 729/21). Y.R. also thanks the Career Advancement Fellowship, Technion, for providing research support.
References
[1] BlogFeedbacck data set. http:://arcchie........gdddahahs hwwp:hhhiiee.iic......eeee//i/ii Accessed: July, 2019.
[2] Facebook comment volume data set.
htttss///archin sereerics.uce htw.a/eeeeeieeeeeeeeee/aaaaattta
Cecebool-Coommeett-oooueeeeeeeee Accessed: July, 2019.
10

---

[3] Medical expenditure panel survey, panel 19. https://mops.ahrq.u htt.a//g////eeewe///aaaaaa
Accessed: July,
2019.
[4] Medical expenditure panel survey, panel 20. httpps://mess.arqqq htw.:/arrgggooeeeeeeeebbaabbaaaa
stats/domnlos http///oooodddddaadaaadeeeaaaaaa setessddeellll.ee httpr//dee.......oooooooubbbbbea Accessed: July, 2019.
[5]
Medical expenditure panel survey, panel 21.
Accessed: July,
2019.
[6]
Physicochemical properties of protein tertiary structure data set. https://archive ics. www..eiw..uiiiiuu /ttp./////////tttsss////// http:///wwccccciaaaorrooooooooor Structure] Accessed: July, 2019.
[7] R. F. Barber, E. Candès, A. Ramdas, and R. J. Tibshirani. Predictive inference with the jackknife+. Annals of Statistics, 49(1): (986:77755 2021.
[8] S. Bates, E. Candès, L. Lei, Y. Romano, and M. Sesia. Testing for outliers with conformal p-values. preprint wrXiv:2104.08277 2021.
[9] M. Cauchois, S. Gupta, and J. C. Duchi. Knowing what you know: valid and validated confidence sets in multiclass and multilabel prediction. Journal of Machine Learning Research, 22(81):1-422 2021.
[10] V. Chernozhukov, K. Wüthrich, and Zhu. Distributional conformal prediction. preprint at 199.::1909007889 2019.
[11] H. A. Chipman, E. I. George, R. E. McCulloch, et al. BART: Bayesian additive regression trees. Annals of Applied Statistics, 4(1):266-29) 2010.
[12] A. Clauset, C. R. Shalizi, and M. E. Newman. Power-law distributions in empirical data. SIAM review, 51(4):661-703 2009.
[13] J. W. Cohen, S. B. Cohen, and J. S. Banthin. The medical expenditure panel survey: national information resource to support healthcare cost research and inform policy and practice. Medical pages $44-S50, 2009,
[14] N. Dalmasso, T. Pospisil, A. B. Lee, R. Izbicki, P. E. Freeman, and A. I. Malz. Conditional density estimation tools in python and r with applications to photometric redshifts and likelihood- free cosmological inference. Astronomy and Computing, 30: 10036,, 2020.
[15] D. Dua and C. Graff. UCI machine learning repository, 2017.
[16] R. Foygel Barber, E. Candès, A. Ramdas, and R. Tibshirani. The limits of distribution-free conditional predictive inference. Information and Inference: A Journal of the IMA, 2019.
[17] C. Gupta, A. K. Kuchibhotla, and A. K. Ramdas. Nested conformal prediction and quantile out-of-bag ensemble methods. preprint at arXiv:1910.10562. 2019.
[18] M. J. Holland. Making learning more transparent using conformalized performance prediction. preprint at ar0iv:2007. 04486, 2020.
[19] R. Izbicki and A. B. Lee. Nonparametric conditional density estimation in high-dimennnn dimensional regression setting, Journal of Computational and Graphical Statistics, 25(4): 1297-136,, 2016.
[20] R. Izbicki, Shimizu, and B. Stern. CD-split: efficient conformal regions high dimensions. preprint at arXiv:2007.12778, 2020.
[21] R. Izbicki, G. T. Shimizu, and R. B. Stern. Flexible distribution -ree conditional predictive bands using density estimators. preprint at 190.0511.0 2019.
[22] B. Kim, C. Xu, and R. Barber. Predictive inference free with the jackknifex-afier--oottt preprint at arXiv:2002.0990 2020.

---

[23] D. Kivaranovic, K. D. Johnson, and H. Leeb. Adaptive, distribution-free prediction intervals for deep networks. In International Conference on Artificial Intelligence and Statisties, pages 4346-436 PMLR, 2020.
[24] J. Lei, M. G'Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman. Distribution-free predictive prediccc inference for regression. Journal of the American Statistical Association, 113(523):1094-11 1111,, 2018.
[25] J. Lei and L.. Wasserman. Distribution-free prediction bands for non-parametric regression. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):11-66, 2014.
[26] M. Magdon-Ismail and A. Atiya. Neural networks for density estimation. In Proceedings of the 11th International Conference Neural Information Processing Systems, pages 522-528, 1998.
[27] N. Meinshausen. Quantile regression forests. Journal of Machine Learning Research, 7:983- 999, 2006.
[28] S. J. Moon, J.-J. Jeon, S. H. Lee, and Y. Kim. Learning multiple quantiles with neural networks. Journal of Computational and Graphical Statistics, pages 1-11, 2021.
[29] Y. Romano, R. F. Barber, C. Sabatti, and Candès. With malice toward none: Assessing uncertainty via equalized coverage. Harvard Data Science Review, 2020.
[30] Y. Romano, E. Patterson, and E. Candès. Conformalized quantile regression. In Advances in Neural Information Processing Systems, pages 3538 3548, 2019.
[31] Y. Romano, M. Sesia, and E. Candès. Classification with valid and adaptive coverage. Advances in Neural Information Processing Systems, 33, 2020.
[32] M. Sesia and E. Candès. A comparison of some conformal quantile regression methods. Stat, 9(1), 2020.
[33] J. W. Taylor. A quantile regression neural network approach to estimating the conditional density of multiperiod returns. Journal of Forecasting. 1944)299-311, 2000.
[34] V. Vovk. Conditional validity of inductive conformal predictors. Asian conference on machine learning, pages 470, 2012.
[35] V. Vovk, A. Gammerman, and G. Shafer. Algorithmic learning in random world. Springer, 2005.
[36] Y. Yang and A. K. Kuchibhotla. Finite-sample efficient conformal prediction. preprint at arXiv:2104118777 13871, 2021.
12


Reviews:
[{"review_id": "WdLi3rQ1m1", "reviewer": "Reviewer_mmq9", "summary": "This paper proposes an extension to conformal prediction that adapts to skewed data, and can achieve better conditional coverage (and provably achieves conditional coverage asymptotically). Conformal prediction in general is a methodology for constructing confidence sets that output likely response candidates $\\widehat{\\mathcal{C}}(X) \\subseteq \\mathcal{Y}$ for and input $X$, rather than a single value. The goal is to ensure that $\\widehat{\\mathcal{C}}$ covers the true response variable, $Y$, with specifiably high probability. While many prior methods provably control *marginal* coverage, *conditional* coverage is a much harder (and more practically important) goal (albeit impossible in finite samples in the general case). This paper attacks this problem by developing a novel conformalization strategy that leverages calibrated estimates of the conditional density $Y \\mid X$ to obtain approximate conditional coverage in finite samples (empirically), and asymptotic conditional coverage (theoretically). ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "=== Strengths ===\n\n- The motivation of the paper is very strong; achieving the harder task of approximate conditional coverage is far more important than marginal coverage for most practical problems. \n\n- The writing for the most part is fairly clear, with few exceptions (see questions below).\n\n- The theoretical analysis is quite comprehensive, with more precise guarantees than several other asymptotically conditional conformal methods.\n\n- Empirically the proposed method performs well relative to baselines. \n\n=== Weaknesses ===\n\n- Using estimates of the conditional distribution have been proposed before, most relevantly perhaps in CD-Split (cited) and HPD-Split (also by Izbicki et. al., 2021). Dist-Split is compared to in this work, though not the other methods (although CD-Split is designed for multi-model targets which are assumed to not exist here, it reportedly seems to work empirically better than Dist-Split even across unimodal tasks). Due to the similarities in approaches, it would be helpful to give better intuition as to why the proposed histogram-based method works better. The comments on line 192 are also confusing in terms of what exactly they mean, as Dist-Split in Izbicki et. al., 2020 has similar guarantees (Thm. 2.5 and Corollary 2.6). More discussion on this would be helpful.\n\n- Achieving conditional coverage relies on the consistency of the histogram-based conditional estimator, which scales linearly with the number of bins (which goes to infinity). It seems that achieving conditional coverage as the number of data points grows is at odds with computational complexity, as opposed to a potentially more scalable density estimator (either way, both need to be calibrated). \n\n- The quality of the conditional estimator will likely degrade considerably in high dimensional problems. The tasks considered in this work are also fairly simple. It would be good to see results on more compelling tasks; i.e. that would illustrate the value of the proposed extension across tasks relevant to the NeurIPS community.\n\n=== Questions ===\n\n- Section 2.2 is a bit dense to follow at first in terms of motivation, and it's not entirely clear what the motivation is for constructing the nested sets based off of the initial $\\bar{t}$. Instead of conformalizing the index of the set per Eq. (12), what is the downside to conformalizing $\\tau$ directly? I.e., one would find the smallest $\\tau$ such that $(1 - \\alpha)\\frac{1 + |\\mathcal{D}^{\\mathrm{cal}}|}{|\\mathcal{D}^{\\mathrm{cal}}|}$ fraction of calibration instances satisfy $y \\in \\underset{(l, u), l \\leq u}{\\arg\\\\!\\min} \\\\{ l - u  \\colon \\sum_{j=l}^{u} \\pi_j(x) \\geq \\tau \\\\}$. It seems that this is due to the randomization strategy in Eq. (8), which should produce slightly more efficient intervals depending on how coarse the binning is?\n\n=== Minor ===\n\n- Typo, line 188 \"increased flexibility\"\n- Supplement, line 16, I believe should have $\\hat{c}$ instead of $c$ in the numerator.\n\n=== Justification for score ===\n\nIn general, this paper is a valuable contribution to the growing body of approaches for approximate conditional conformal inference. There are a couple of clarity/differentiation concerns that could be addressed in the writing, and the method is similar to several already existing approaches (e.g., extending Romano et. al., 2020 to a regression setting, including the binning strategy as compared to Izbicki et. al., 2020)---but overall the work is quite solid.\n\n=== Citations ===\n- Izbicki et. al., 2020. Flexible distribution-free conditional predictive bands using density estimators.\n- Izbicki et. al., 2021. CD-split and HPD-split: efficient conformal regions in high dimensions.\n- Romano et. al., 2020. Classification with Valid and Adaptive Coverage.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I don't see any potential negative societal impacts other than those inherited by the underlying wrapped algorithm this method gives confidence intervals for. Conditional coverage is overall more desirable than marginal coverage in terms of mitigating risks like fairness. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}, {"review_id": "IPyeeGZw25D", "reviewer": "Reviewer_P36z", "summary": "This paper introduces a new method to construct conformal prediction intervals based on histograms of the conditional distribution of an outcome variable. Given a histogram of the conditional distribution of the outcome variable the method finds the shortest interval whose associated mass is no less than the desired coverage probability of the prediction interval. This yields intervals that automatically adapt to the skewness of the data. The intervals have provably correct marginal coverage in finite samples and correct conditional coverage and optimal length in large samples. A simulation study and numerical experiments on several benchmark data sets corroborate the theoretical results and demonstrate the advantage over several competing methods.\n", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This a well-written and clearly organized paper that makes an interesting contribution to conformal prediction. The paper addresses the important problem of how to construct conformal confidence sets that have finite-sample marginal coverage and large-sample conditional coverage while being as narrow as possible. Typically, adaptivity and optimality of (conformal) confidence sets are only argued heuristically (by judiciously choosing tuning parameters, e.g. Lei and Wasserman, JRSSB 2014) and are often excluded from the theoretical analysis. The novelty of the proposed method is that adaptivity and optimality are guaranteed.\n\nTo me the strength of the paper lies in the simple and intuitive idea of the proposed method (setting aside the somewhat tedious actual implementation). I believe that the basic idea of directly minimizing the length of the prediction interval subject to a constraint on the coverage probability can be extended and adapted to related problems. Thus, the paper does not simply propose a solution to an important problem, but also provides food for thought (see below).\n\nAn obvious weakness of the paper are the strong assumptions for the theoretical analysis. I am fine with Assumptions 1-3, but I would have preferred if the authors had done without Assumption 4 and 5 on the unimodality of the true and estimated conditional density of the outcome variable. At the very least, the authors could have included empirical results from an outcome variable with non-unimodal conditional distribution.\n\nI wonder whether the authors have thought about the following extension of their approach:  Why \"only\" consider prediction intervals of shortest length and not more general prediction regions with smallest Lebesgue measure? It seems to me that opt. problem (7) can be easily modified in such a way that even non-unimodal conditional distributions of the outcome variable can be handled. Consider solving\n\n$\\mathcal{S}(x, \\pi, S^-, S^+, \\tau) :=  \\arg\\min_{M\\subseteq \\\\{1,\\ldots, m\\\\}} \\left\\\\{ |M| :  \\sum_{j \\in M} \\pi_j(x) \\geq \\tau, S^- \\subseteq M \\subseteq S^+ \\right\\\\} \\quad{}\\quad{} (*)$\n\nIf the histogram of the conditional density of the outcome variable is unimodal, then $\\mathcal{S}(x, \\pi, S^-, S^+, \\tau)$ is an interval. Otherwise it is a collection of convex sets. Therefore, (\\*) seems to be even more adaptive than the original opt. problem (7).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, the authors have adequately addressed the limitation of their assumptions and the lack of control of lower and upper miscoverage rates.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "F3gJsZUVtmf", "reviewer": "Reviewer_prSn", "summary": "This paper proposes a new nonconformity-score based on conditional density estimates P(Y|X). The main idea is to approximate the optimal prediction interval based on the conditional density estimate. The proposed non-conformity score is experimentally assessed using the conformalization technique of split conformal. The overall method provably achieves marginal coverage, performs comparably to other methods for conditional coverage, and achieves the shortest prediction intervals on all datasets. ", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "My overall opinion of this paper is quite positive. The method makes sense and the experiments indicate a significant improvement on the efficiency/length of the prediction intervals. Thus I think the paper is a useful contribution to the ML community. In the following 'sections' I elaborate further on things that I liked about the paper, along with some suggestions for further improvements.\n\nHowever, I have one methodological question that is almost 'begging to be answered'. The optimal prediction **set** given P(Y|X) is simply to take the highest predictions until 1-\\alpha coverage is reached (Appendix F [2]). In the case of viewing it from a 'histogram' perspective, one would take the largest histogram bars until 1-\\alpha is reached. Perhaps considering the optimal prediction **interval** is more interpretable in some applications, as the authors do. However, in some applications, I think sets could also be ok? Thus, I think it makes sense to also explore if taking prediction sets helps reduce the widths. If the authors found in experiments that it does not help, could the authors briefly comment on this? Also, what if one takes the convex-hull interval of the PS to make it a PI? \n\nOn a related note, lines 43-46 do not answer the following question. P(Y|X) is not known, but after an estimate of P(Y|X) is computed, why is it discretized into bins? As far as I can tell, the P(Y|X) models that the authors use essentially have real-valued output. Why not work with the full real-valued estimate and define nested sets for them? \n\n*Writing and clarity:*\nThe introduction was clean and easy-to-follow. By the end of page 2, I already had a broad idea of what to expect from the rest of the paper. \nI have a significant concern regarding the density of notation and ideas introduced in Sections 2.2 and 2.3. In lines 165–167, the authors say: “In particular, we analyze a slightly modified version of Algorithm 1 in which there is no randomization; this is theoretically more amenable and equivalent in spirit, although it may yield wider intervals in finite samples.” I fully agree. In light of this, is there a reason the authors chose to describe the full randomization strategy along with the main conformal algorithm in Section 2.2 (rather than in a future section or Appendix)? In my opinion, Section 2.2 has way too much notation and many nitty-gritties of conformal inference are introduced very densely…. I think this could make the current description of the method *quite* difficult to grasp for someone unfamiliar with conformal prediction. \n\n*Originality:*\nThe main idea of using P(Y|X) estimates to produce efficient conformal sets has appeared in different forms, and is not striking. In this sense, the initial prediction sets that the authors describe (Sec 1.2) are natural. However, they pose the issue of not being nested sets making it tricky to apply conformal prediction. I really liked the way the authors translated this setup into a nested problem by first fixing a single set (eq. (9)), and then considering supersets and subsets of this initial set (eq. (10)). The experimental section verifies that this is a useful description of nested sets for producing short prediction intervals. \n\n*Terminology:*\nIn the nonparametrics statistics literature, histogram classification and histogram regression have been used to refer to methods that partition the feature space X into a number of bins, and compute a single prediction for each bin [1, and many others]. To the best of my understanding, there is very little intersection between the ideas proposed by the present authors, and the histogram regression literature. Due to this, I strongly suggest changing the title of the paper and the name of the proposed method. Something as simple as “Conformal inference using conditional histograms” avoids the confusion, and avoids overstepping on another interesting body of work. \n\nMinor comments: \n- Space permitting, it would be great if Assumptions 1–5 are latex enumerated.\n- It may be useful to demonstrate that the prediction interval based on the unconformalized conditional density estimate (based on Sec 1.2) does not have valid coverage. \n\nReferences:\n \n[1] https://www.jstor.org/stable/2242583\n\n[2] https://arxiv.org/pdf/1910.10562.pdf", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "7M9urTnc4Pg", "reviewer": "Reviewer_U8cY", "summary": "This paper considers constructing prediction set using histogram estimation of conditional density. The results and methods are new and interesting. ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Prediction set based on estimated density seems a natural and fundamental problem. The current paper provides a novel method and theory regarding this problem. Paper writing is clear and results are significant.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The term \"conformal\" needs to be further explained. As to common readers, the current approach is a \"plug-in\" method since unknown density is replaced by its estimator. Authors should explain the explicit reason of using \"conformal\" to name this procedure.\n\nSecond, can $Y_i$ be multi-dimensional? \n\nThird, can $f$ be estimated using general kernel density estimator?\n\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": ["I don’t know"], "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "6Q0EddXbPXi", "reviewer": "Reviewer_ucSx", "summary": "The paper shows how to obtain *short* prediction intervals that approximately have the right conditional coverage. In order to do so, the method estimates the conditional distribution of a new label given its features, and then uses a conformal score based on the quantiles implied by such estimate to create the prediction intervals. The experiments show that the method outperforms other quantile-based conformal methods in terms of width while still approximately controlling conditional coverage. Some theoretical results that prove converge to the oracle (i.e., the shortest prediction interval) are also shown.", "questions": "", "limitations": "", "rating": 9, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well written and technically sound. Given the importance of accurate uncertainty quantification methods, the paper is also very relevant.\n\nAs far as I know, this is the first paper on conformal methods that attempt to find *short* prediction *intervals*; most conformal-quantile approaches are based on intervals that have the form $(q_\\alpha,q_{1-\\alpha})$, which are only the shortest intervals if the underlying distribution is symmetric. Thus, the paper has an important and novel goal.\n\nAs the authors mention, there has however been work on how to recover short prediction *regions* (which may not be intervals). Although such regions will be smaller than the oracle intervals if the underlying distribution is **not** unimodal, the fact that the CHR yields intervals is an advantage in some settings (for instance, they are easier to report). This advantage should be emphasized.\n\nRelated to this point, the conformity score (Eq. 12) is essentially (up to the randomization need due to the discrete nature of the estimated density) the HPD-split score introduced by [20] (https://arxiv.org/abs/2007.12778). Indeed, the score is (a monotonic transformation of) the area under a region with high density. The main difference is that CHR has an additional constrain that the output regions must be intervals. Indeed, Theorem 2 is very similar to Theorem 22 of [20], but apparently, the HPD-split method does not need to satisfy Assumption 4 (unimodality) to recover the oracle region. Thus, in principle, HPD-split could give even smaller regions (as they don't need to be intervals), while still approximately controlling conditional coverage. It would be interesting to see how it compares to CHR in unimodal settings.\n\nIt also seems to me that the conformity score in (12) could actually be computed using any conditional density estimator (although Eq. 7 would have to be slightly generalized to mimick Eq. 4). That is, there is no need for the method to be restricted to the histogram-based approach (although such an approach is interesting). In other words, the CHR score is agnostic to the density estimator being used, which is an attractive advantage.\n\nMinor comment: It is not clear to me how conditional coverage was computed/estimated on real datasets. It doesn't seem to be trivial, because it is a conditional probability per se.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations have been properly addressed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "-Pw0iD1I2A_", "reviewer": "Reviewer_QWwn", "summary": "This paper introduces a conformity score that aims to decrease average interval lengths and improve conditional coverage from a black-box estimate  $\\hat{P}(y \\mid x)$ for split conformal prediction. The method, conformal histogram regression (CHR), involves first binning the space of $Y$, resulting in a conditional histogram from which approximate oracle intervals can be computed. A nested sequence of these intervals is then created for a sequence of predictive miscoverage values $\\tau$, where $\\tau$ will be close to $\\alpha$ if $\\hat{P}(y \\mid x)$ is a good estimate. The value of $\\tau$ is selected through a conformity score, and the authors show their method obtains finite marginal coverage and asymptotic conditional coverage. The authors then demonstrate the method on a few examples. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths: \n\nThis paper is an interesting natural extension of [Romano et al., 2020] from discrete to continuous $Y$, and also connects to the nested conformal interpretation of [Gupta et al., 2019]. Although the conformity score has been previously introduced, the construction of the nested prediction sets through histograms is novel. The asymptotic results are also stronger than previous works, where both asymptotic oracle length and conditional coverage are attained. The empirical results demonstrate superior conditional coverage and average interval lengths to other methods, especially for skewed data. The paper is well-written and enjoyable to read.\n\n#####################################################################\n\nWeaknesses/Questions:\n\nI only have a few minor points:\n\n1.) For equation (7), does treating $|u-l|$ as the length require the bins to be equally spaced? I don't think this is stated.\n\n2.) It may be good to briefly mention the negligible computational cost of CHR (which is in the appendix) in the main paper to help motivate the method. A rough example of some run-times in the experiments may also be useful for readers looking to apply the method.\n\n3.) Just a few typographical/communication points:\n- I found Section 2.2 slightly difficult to read, as the notation gets a little heavy. This may not be necessary, but the authors could consider presenting the nested intervals without randomization (e.g. after Line 119), with the randomization in the Appendix, as it is not needed in Theorem 2. This would give more room for intuitive discussions, related to my next point.\n- It may be helpful to introduce some intuition on the conformity score in equation (12) and why we need the sets to be nested for readers unfamiliar with previous work, perhaps at the start of Section 2.3.\n- Line 113: $\\epsilon$ is mentioned here before it is defined\n- Line 188: 'increased' instead of 'increase'\n\n#####################################################################\n\nOverall:\n\nThis paper is an interesting extension of previous work, and the provided asymptotic justifications of attaining oracle width and conditional coverage is useful. The method is also general and can empirically provide better average widths and conditional coverage than other methods, particularly under skewed data, making it useful in practice. \n\n#####################################################################\n\nReferences:\n\nRomano, Y., Sesia, M., \\& Candes, E. (2020). Classification with Valid and Adaptive Coverage. Advances in Neural Information Processing Systems, 33, 3581-3591.\n\nGupta, C., Kuchibhotla, A. K., \\& Ramdas, A. K. (2019). Nested conformal prediction and quantile out-of-bag ensemble methods. arXiv preprint arXiv:1910.10562.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have described the limitations of their method - in particular their method does not control for upper and lower miscoverage, and they provide alternative recommendations.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper presents an original contribution for conformal prediction with optimized bandwidth. The paper constructs conformal confidence sets that have finite-sample marginal coverage and large-sample conditional coverage while being as narrow as possible. This, I believe, will be a valuable addition to the current conformal prediction / regression problems, as confirmed by reviewers.
Decision: Accept (Spotlight)