Learning Prediction Intervals for Regression: Generalization andCalibration

Haoxian Chen
IEOR
Columbia University

Ziyi Huang
EE
Columbia University

Henry Lam
IEOR
Columbia University

Abstract
We study the generation of prediction intervals in regression for uncertainty quantification. This task can be formalized as
an empirical constrained optimization problem that minimizes the average interval
width while maintaining the coverage accuracy across data. We strengthen the existing
literature by studying two aspects of this empirical optimization. First is a general learning theory to characterize the optimalityfeasibility tradeo↵ that encompasses Lipschitz continuity and VC-subgraph classes,
which are exemplified in regression trees andneural networks. Second is a calibration machinery and the corresponding statistical theory to optimally select the regularization parameter that manages this tradeo↵, which bypasses the overfitting issues in previous approaches in coverage attainment. We empirically demonstrate the strengths of our interval generation and calibration algorithms interms of testing performances compared toexisting benchmarks.

1

Introduction

While most literature in machine learning focuses on
point prediction, uncertainty quantification plays, arguably, an equally important role in reliability assessment and risk-based decision-making. In regression, anatural approach to quantify uncertainty is the prediction interval (PI), namely an upper and lower limit fora given feature value X that covers the corresponding
outcome Y with high probability. The interval center
Proceedings of the 24th International Conference on Artificial Intelligence and Statistics (AISTATS) 2021, San Diego,
California, USA. PMLR: Volume 130. Copyright 2021 by
the author(s).

Huajie Qian
Alibaba Group

Haofeng Zhang
IEOR
Columbia University

represents the expected outcome, whereas the width
represents the uncertainty. A test point with a highexpected outcome, but wide PI width, means that theoutcome can still be low with a significant chance, thus
signifies a downside risk that should not be overlooked.
In this paper, we study the construction of PIs thatsatisfy an overall target prediction level across data,
known as the expected coverage rate (Rosenfeld et al.,2018). Compared to widely used conditional (on X)
coverage rate, this notion is advantageously more tangible to measure and easier to achieve. This meansthat a much wider class of models can be trained tobuild PIs, as less conditions are needed to impose on
the true relation and the model class to obtain satisfactory guarantees. In general, constructing a good PI in
this framework requires balancing a tradeo↵ between
the expected interval width and coverage maintenance,
which can be formalized as an empirical constrained
optimization. This viewpoint has been used in Khosravi et al. (2010) and Pearce et al. (2018) that focus onneural networks, Rosenfeld et al. (2018) that studies adual formulation, and Galván et al. (2017) that uses
multi-objective evolutionary optimization. It also relates to the learning of minimum volume sets (Polonik,
1997; Scott and Nowak, 2006) in which a similar tradeo↵ between volume and probability content appears.
Building on these works, our goal in this paper is to
study two key inter-related statistical aspects of thisempirical constrained optimization that enhances previous results both in theory and in practice:
Feasibility-Optimality Tradeo↵ for Interval
Models. We develop a learning theory for the PIs
constructed from empirical constrained optimizationthat statistically achieves both feasibility (coverage)
and optimality (interval width). Methodologically, webuild a general “sensitivity measure” that controls thistradeo↵, which in turn requires developing deviation
bounds for simultaneous empirical processes. Our theory in particular covers the Lipschitz continuous model
class (in parameter) and finite Vapnik–Chervonenkis(VC)-subgraph class, exemplified by a wide class ofLearning Prediction Intervals for Regression: Generalization and Calibrationneural networks and regression trees. Such type of
joint coverage-width learning guarantees appears thefirst in the literature. It expands the coverage-only results and the considered model classes in Rosenfeld
et al. (2018). It also generalizes Scott and Nowak(2006) as both our constraints and objectives possess
extra sophistication related to the shape requirementof the set as an interval, and also we characterize feasibility and optimality attainments explicitly insteadof the implicit metric in Scott and Nowak (2006).
Calibration Method and Performance Guarantees. We propose a general-purpose, ready-toimplement calibration methodology to guarantee overall PI coverage with prefixed confidence. This approach is guided by a novel utilization of the highdimensional Berry-Esseen theorem (Chernozhukov
et al., 2017). It is designed to combat the overfitting issue of interval models and perform accurately
on the test set. We demonstrate empirically how ourapproach either outperforms other methods in termsof achieving correct coverages or, for those methods
with comparable coverages, we attain shorter interval
widths. Moreover, our approach applies, with little
adjustment, to accurately construct multiple PIs atdi↵erent prediction levels simultaneously. This adds
extra flexibility for decision-makers to construct PIs
without needing to pre-select the prediction level inadvance.

2

Related WorkWe first review two most closely related methods, andthen move on to other works.
Conformal Learning (CL). First proposed in Vovk
et al. (2005), conformal learning (CL) is a class of
methods that leverage data exchangeability to constructs PIs with finite-sample and distribution-freecoverage guarantees. The original CL requires retraining for each possible test point and is therefore computationally prohibitive in general. Split/inductive CL(Papadopoulos, 2008; Lei et al., 2015, 2018) improves
the computational efficiency based on a holdout validation that avoids retraining, but at the cost of higher
variability and wider intervals due to less efficient data
use. Lying in between are variants based on more efficient cross-validation schemes, including leave-oneout (or the Jackknife; Barber et al. (2019); Alaa andvan der Schaar (2020); Steinberger and Leeb (2016)),K-fold (Vovk, 2015) and ensemble methods (Gupta
et al., 2019; Kim et al., 2020). Recently, quantile regression are combined with CL (Kivaranovic et al.,2020; Romano et al., 2019) to take into account theheterogeneity of uncertainties across feature values.
Despite its generality, the coverage guarantees fromCL are only marginal with respect to the training data(except split CL (Vovk, 2012)), whereas our proposed
calibration method provides a stronger high confidence
guarantee. Moreover, our approach explicitly optimizes the interval width, therefore typically generates
shorter PIs than CL.
Quantile Regression (QR). Quantile regression(QR) estimates the conditional quantiles of Y thatcan be used to construct PIs. Classical QR methods require strong assumptions (e.g., linearity or otherparametric forms; Chapter 4 in Koenker and Hallock (2001)). Approaches that relax these assumptions include quantile regression forests (QRF) (Meinshausen, 2006) and kernel support vector machine(SVM) (Steinwart et al., 2011). However, little is
known about their finite-sample coverage performance
because of estimation errors in the quantiles. Recently, Kivaranovic et al. (2020) proposes calibrating
the weight parameter in the pinball loss on a holdout data set to enhance PI coverage. This calibration
scheme, however, does not address overfitting on the
holdout set, thus could fall short of providing correctcoverages on test data as our experiments show.
Other Approaches. PI construction has been substantially studied in classical statistics. To understandthis construction, the error of a (point) prediction cantypically be decomposed into two components: model
uncertainty, which comes from the variability in the
training data or method, and outcome uncertainty,
which comes from the noise of Y conditional on X.
The classical literature often assumes well-defined andsimple forms on the relation between X and Y (e.g.,linear model, Gaussian; Seber and Lee 2012). In thiscase, model uncertainty reduces to parameter estimation errors. Outcome uncertainty, on the other hand,is intrinsic in the population distribution but not thetraining, i.e. it arises even if the model is perfectlytrained. An array of methods account for both sourcesof variability, which utilize approaches ranging fromasymptotic normality (Seber and Lee, 2012), deconvolution (Schmoyer, 1992), and resampling schemes suchas the bootstrap (Stine, 1985) and jackknife (Steinberger and Leeb, 2016).
To overcome the strong assumptions in classical statistical models, several model-free approaches have beendeveloped. Nonparametric regression, such as spline or
kernel-based methods (Doksum and Koo, 2000; Olive,
2007), removes rigid model assumptions but at the expense of strong dimension dependence. Gaussian processes or kriging-based methods (Sacks et al., 1989),
popular in the areas of metamodeling and computer
emulation, regard outcomes as a response surface andperform Gaussian posterior updates. In particular,
stochastic kriging (SK) model (Ankenman et al., 2010)

Haoxian Chen, Ziyi Huang, Henry Lam, Huajie Qian, Haofeng Zhang

constructs PIs that account for both model and output variabilities by using a prior correlation structurethat entails both. However, SK does not deliver afrequentist coverage guarantee nor convergence rate.
More recently, uncertainty quantification of neural networks regarding their model and output variabilities
are studied, via methods such as the delta method
and the bootstrap (Papadopoulos et al., 2001; Khosravi et al., 2011). Nonetheless, like in classical statistical models, these approaches can only capture variability due to data and training noises, but not thebias against the true relation.
Lastly, our PI construction follows the high-quality criterion in works including Khosravi et al. (2010, 2011);
Galván et al. (2017); Pearce et al. (2018); Rosenfeld
et al. (2018); Zhang et al. (2019); Zhu et al. (2019),
which propose various loss functions to capture thewidth-coverage tradeo↵. They are also related to thehighest density intervals in statistics (Box and Tiao,2011). Our investigations in this paper provide theoretical guarantees in using this criterion.

3

PI Learning as Empirical
Constrained Optimization

Pn
E⇡ˆX [U (X) L(X)] = n1 i=1 (U (Xi ) L(Xi )). The
adjustment t, which can be viewed as a penalty term
for the empirical coverage rate, is used to boost thegeneralized coverage performance for the optimal interval solved from dopt(t): If no adjustment is made in
the constraint (t = 0), then, because of noise, the true
coverage can be lower than 1 ↵ with significant probability even if the empirical coverage is above 1 ↵. A
positive t decreases the probability of such an event.
Choosing t too large, however, would eliminate moreintervals from the feasible set, leading to a deterioration in the obtained expected width (objective). Oneof our main investigations is to balance the coverage
and width performance by judiciously choosing t.
We point out that, while our focus is on training L
and U directly, our approach can also be applied naturally when we are given in advance a well-trained
point predictor f : X ! R (obtained by any means
independent of D). In this case we may seek two
non-negative functions i : X ! [0, 1) such that
[L(x), U (x)] = [f (x)1 (x), f (x) + 2 (x)]. Our subsequent development applies by constructing lower andupper bounds for the “translational” data set D̃ :={(Xi , Yi f (Xi ))}i=1,...,n .

We consider the general regression setting where X 2
X ⇢ Rd is the feature vector and Y 2 Y ⇢ R is the outcome. Given an i.i.d. data set D := {(Xi , Yi )}i=1,...,neach drawn from an unknown joint distribution ⇡, ourgoal is to find a PI defined as:
Definition 3.1. An interval [L(x), U (x)], where both
L, U : X ! R, is called a prediction interval (PI)with (overall) coverage rate 1 ↵ (0  ↵  1) if
P⇡ (Y 2 [L(X), U (X)])
1 ↵, where P⇡ denotes
the probability with respect to the joint distribution ⇡.

4

We aim to find two functions L and U , both in a hypothesis class H. To learn the optimal high-quality PI
that attains a given coverage rate 1 ↵, we consider
the following constrained optimization problem:

opt(t) :min

L,U 2H and LUE⇡X [U (X)

L(X)]

subject to P⇡ (Y 2 [L(X), U (X)])

(3.1)
1

↵

where E⇡X denotes the expectation with respect tothe marginal distribution of X. Given the data D, weapproximate (3.1) with the following empirical constrained optimization problem

Joint Coverage-Width Learning
Guarantees

We establish finite-sample generalization error bounds
for two major classes: 1) finite VC dimensions, and 2)
Lipschitz continuous in parameters, by building on aunified “sensitivity bound” on the oracle optimization.
Corresponding results on consistency are provided inAppendix A. To begin, let Rt⇤ (H) be the optimal valueofmin

L,U 2H and LUE⇡X [U (X)

L(X)]

subject to P⇡ (Y 2 [L(X), U (X)])

1 ↵+t
(4.1)
which is (3.1) but with a higher target prediction level,
and correspondingly R⇤ (H) be the optimal value of(3.1). We make the following assumptions on the hypothesis class H, and the conditional distribution of Y
given X:
Assumption 1. For every function h 2 H, we haveh + c 2 H for every constant c 2 R.

Assumption 2. The conditional distribution of Y
given X = x has a density p(·|x) for every x 2 X .
L,U 2H and LU
Moreover, for every x, y such that P⇡ (Y  y|X = x) 2
subject to P⇡ˆ (Y 2 [L(X), U (X)]) 1 ↵ + t (0, 1), we have p(y|x) > 0 and that p(·|x) is continuous
(3.2)
at y.
parameterized by t 2 [0, ↵], where E⇡ˆX , P⇡ˆ are expectation and probability with respect to the empirThe simple closedness property in Assumption 1 turns
ical distribution constructed from the data D, e.g.,out to allow sufficiently tight and tractable analysis fordopt(t) :min

E⇡ˆX [U (X)

L(X)]

Learning Prediction Intervals for Regression: Generalization and Calibration

many, potentially complicated, function classes under
the same roadmap. Many common classes (e.g., linear, piece-wise constant such as tree-based models, andneural networks with linear activation in the outputunit) satisfy Assumption 1. It is also straightforward
to enforce a class to satisfy Assumption 1 by simply
adding one extra parameter. Assumption 2 is a mild
non-degeneracy condition on the conditional distribution (e.g., Assumption 2 holds when p(·|x) is Gaussianor uniform over a certain interval for each x).
We have the following upper bound for Rt⇤ (H) R⇤ (H)for 0 < t < ↵:Theorem 4.1 (Sensitivity bound with respect to target prediction level). Suppose Assumptions 1-2 hold.
For every x 2 X and
2 (0, 1), let (x, ) :=inf p(y1 |x) + p(y2 |x) : P⇡ (y1  Y  y2 |X = x) 1
and:= sup{z > 0 : P⇡X ( (X, ) < z)  }.
Then, for every ↵ 2 (0, 1) and t 2 (0, ↵), we haveRt⇤ (H) R⇤ (H)  6t/((↵ t) (↵ t)/3 ).
We briefly explain how Theorem 4.1 helps develop generalization guarantees. Let Ht2 , Ĥt2 ⇢ H2 := H ⇥ H
be the feasible set of opt(t) and dopt(t) respectively.
If a uniform error bound t over the class H2 canbe established for the empirical coverage constraint in22(3.2), then Ht+t ⇢ Ĥt , therefore the shortest interval we can potentially learn from dopt(t) can only be
wider than the oracle optimum from (3.1) by at most⇤⇤Rt+t (H) R (H). The density lower bounds in Theorem 4.1 ensure a sufficient increase of the coverage
probability as interval width increases, which inversely
constrains the growth of interval width as the required
coverage increases.
To derive finite-sample convergence rate, we first assume the availability of certain deviation bounds forthe empirical coverage rate and interval width that appear in dopt(t):

Theorem 4.2 (Rate of convergence). Suppose Assumptions 1-2 hold, and the following deviation bounds hold for every ✏, t > 0:
P suph2H |E⇡ˆX [h(X)] E⇡X [h(X)]| > ✏  1 (n, ✏, H)
and P supL,U 2H and LU |P⇡ˆ (Y 2 [L(X), U (X)])
P⇡ (Y 2 [L(X), U (X)])| > t  2 (n, t, H). Thenfor every t 2 (0, ↵2 ) and ✏ > 0, with probability at least 1
1 (n, ✏, H)
2 (n, t, H), we have,for every optimal solution (L̂⇤t , Ût⇤ ) of dopt(t),
that P⇡ (Y 2 [L̂⇤t (X), Ût⇤ (X)])
1↵ andE⇡X [Ût⇤ (X) L̂⇤t (X)]  R⇤ (H) + (↵ 2t)12t+ 4✏.
(↵ 2t)/3
Theorem 4.2 translates the deviation bounds of two
empirical processes into the probability of jointly
achieving optimality and feasibility. With this, ourfocus is to derive deviation bounds for important hypothesis classes. Such bounds are well-understood in

the literature, e.g., Chapter 2.14 in Van der Vaart andWellner (1996) and Chapter 3 in Vapnik (2013). However, in our case, we need to go beyond the standard
theory to control two empirical processes simultaneously by choosing a single function class H. Next we
present two fairly general choices of H for which exponential deviation bounds can be obtained for bothprocesses.
To present the results, we introduce some terminologies. Let vc (S) be the VC dimension of a class S ofsets. The VC dimension vc (G) of a class G of functions from X to R is defined to be the VC dimensionof the set of subgraphs SG := {{(x, z) 2 X ⇥ R :z < g(x)} : g 2 G}. G is called VC-subgraph if
vc (G) < 1. For a vector, k · kp represents its lp -norm
for p
1. For a function ⇠ : X ! R, we denote byk⇠k 2 := inf{c > 0 : E⇡X [exp(⇠ 2 (X)/c2 )]  2} its subGaussian norm under the distribution ⇡X , and denote
1/pby k⇠kp := E⇡X [|⇠(X)|p ]its Lp norm.
Theorem 4.3 (Joint coverage-width guarantee for
VC-subgraph class). If the hypothesis class H is suchthat the augmented class H+ := {h + c : h 2 H, c 2
R} is VC-subgraph, and H(x) := suph2H |h(x)E⇡X [h(X)]| satisfies kHk 2 < 1, then the deviation
bounds in Theorem 4.2 satisfy
1 (n, ✏, H)  2 exp

⇣

⌘
n✏2
,
CkHk 2 vc (H+ )2

8<4n+1 exp( t2 n)⇣⌘Cvc(H)
2 (n, t, H) 
2en
:4 Cvc(H)exp( t2 n)

if n < Cvc(H)2
if n

Cvc(H)
2(4.2)

where C is a universal constant, and e is thebase of the natural logarithm. If Assumptions 1-2
also hold (in which case H = H+ ), then for evCvc(H)ery ⌘ 2 (0, 1), when n
and we set t :=2q
Cvc(H)
18
2en↵
log Cvc(H)  4 , with probability at
n log ⌘ +nleast 1 ⌘, we have, for every optimal solution (L̂⇤t , Ût⇤ )of dopt(t), that P⇡ (Y 2 [L̂⇤t (X), Ût⇤ (X)]) 1 ↵ andE⇡X [Ût⇤ (X) L̂⇤t (X)]  R⇤ (H)s
Cvc (H)
1
8
2en
24
+
log +
log
↵ ↵6
n
⌘
n
Cvc (H)svc (H+ )4+· 16CkHk2 2 log .
n⌘Theorem 4.3 reveals that, after ignoring logarithmic factors, the sample size n needed to learn agood PI with guaranteed coverage pfrom dopt(t) isof order ⌦(vc (H)), if t of order O vc (H) /n is
adopted.
p The corresponding optimality gap (in width)
is O vc (H+ ) /n . Appendix B provides further discussion regarding the use of H+ versus H in the bound.
Similar sample complexities for VC-major H have been

Haoxian Chen, Ziyi Huang, Henry Lam, Huajie Qian, Haofeng Zhang

proposed in Rosenfeld et al. (2018) for a formulation that can be viewed as the dual of (3.1), where
the coverage rate is maximized under a mean width
constraint. Comparing VC-major and VC-subgraphclasses, they both cover many hypothesis classes commonly used in practice, e.g., linear functions, regression trees, and neural networks that are to be discussed momentarily. Nonetheless, in general neither
VC-major nor VC-subgraph implies the other, andtherefore our VC-subgraph results are in parallel tothose in Rosenfeld et al. (2018). More importantly,
their results provide finite-sample errors only for thecoverage rate, but not the interval width, whereas
we characterize performances jointly in coverage andwidth, thanks to our novel sensitivity measure from
Theorem 4.1. Moreover, our results appear to bypass
a technical issue of Rosenfeld et al. (2018). A key result there is that for a VC-major class H, the induced
set of between-graphs {{(x, z) : L(x)  z  U (x)} :L, U 2 H and L  U } is a VC class. When the class
H is uniformly bounded from below, say 0, then thisis equivalent to H being VC-subgraph. However, acounter-example for this conclusion is constructed in
Theorem 2.1, statement f, in Dudley (1987). Our approach overcomes such technical ambiguities.
Our next considered hypothesis class is on Lipschitz
continuity with respect to the class parameter:
Theorem 4.4 (Joint coverage-width guarantee for the
Lipschitz class in parameter). Suppose H = {h(·, ✓) :✓ 2 ⇥} where the parameter space ⇥ is a bounded
set in Rl . If the functions are Lipschitz continuous in the parameter, i.e., |h(x, ✓1 )h(x, ✓2 )| L(x)k✓1 ✓2 k2 , for all ✓1 , ✓2 2 ⇥, x 2 X for some
L : X ! R such that kLk2 < 1, and H(x) :=sup✓2⇥ |h(x, ✓) E⇡X [h(X, ✓)]| satisfies kHk 2 < 1,then the deviation bounds in Theorem 4.2 satisfy

1 (n, ✏, H)  2 exp

⇣

✏2 n

⌘

2
C max{log diam(⇥)kLk, 1}kHk2 2 l
kHk2
⇣
⌘2l
t2 n
2C log log CH· max{, 1}
exp( 2t2 n)
2 (n, t, H)  CH
2Cl

where diam (⇥) := sup✓1 ,✓2 2⇥ k✓1 ✓2 k2 is the diameter of ⇥, DY |X := supx,y p(y|x) is the supremum of
the conditional density, CH := 12diam (⇥) DY |X kLk1 ,and C is a universal constant. If Assumptions 12 also hold,q then for every ⌘ 2 (0, 1), when we
12l↵
set t :=
n log ⌘ + n · 4C log(CH ) log log(CH )  4 ,with probability at least 1⌘, we have, for every
optimal solution (L̂⇤t , Ût⇤ ) of dopt(t), that P⇡ (Y 2,[L̂⇤t (X), Ût⇤ (X)])

1

↵ and

E⇡X [Ût⇤ (X) L̂⇤t (X)]  R⇤ (H)
s
24
1
2
l
+
log + · 4C log(CH ) log log(CH )
↵ ↵6
n
⌘
n
s
diam (⇥) kLk2
l
4+
· 16C max log, 1 kHk2 2 log .
n
kHk2
⌘
(4.3)

The Lipschitz class is beyond the scope of Rosenfeld
et al. (2018). Theorem 4.4 states that the required
sample size for this class to achieve a certain learning
accuracy is of order ⌦(l log(kLk1 )), which depends onthe dimension of the parameter space and the Lipschitz coefficient. Correspondingly,
the optimality gap
pon the interval width is O( l log(kLk2 )/n). Here thelogarithmic factor associated with the Lipschitz coefficient is significant because kLk2 can be exponential in the number of layers for deep neural networks(see more details below). Note that our Lipschitzness is in the class parameter ✓ rather than the input x. The latter has recently been used to regularize
neural networks to improve generalization gaps (e.g.,
Bartlett et al. (2017); Yoshida and Miyato (2017);
Gouk et al. (2018)) and robustness against adversarial attacks (e.g., Cisse et al. (2017); Hein and Andriushchenko (2017); Tsuzuku et al. (2018)), but canpotentially lead to a loss in expressiveness of the network (Huster et al., 2018; Anil et al., 2019) becauseof the size restriction on network weights. Our result
does not restrict the sensitivity of the network in the
inputs, and in turn its expressiveness.
To further distinguish the technical novelty of Theorem 4.3 and Theorem 4.4, note that established results on uniform convergence bounds (UCBs) assuming VC or Lipschitzness on the hypothesis class H can
only handle the objective (width) on H, but not theconstraint (coverage) on a di↵erent hypothesis classof indicator functions on the joint (x, y) space. Our
analysis explicitly controls this constraint complexity to achieve joint optimality-feasibility guarantees.
The proof of Theorem 4.3 involves bounding the VCdimension of the indicator class through intersectionsof VC set classes, while that of Theorem 4.4 requires
directly bounding the bracketing number and evaluating the entropy integral. Moreover, even for the objective, it appears that the explicit UCBs with potentiallyunbounded VC or Lipschitz classes considered here are
new to the best of our knowledge.
We showcase the application of Theorems 4.3 and 4.4
in two important classes: Regression trees and neuralnetworks. Appendix C presents an additional class of
linear hypothesis.
Regression Tree.

Suppose we build two binary

Learning Prediction Intervals for Regression: Generalization and Calibration

trees to construct L, U respectively. The tree has at
most S + 1 terminal nodes, and each non-terminal
node is split according to xi⇤  q or xi⇤ > q forsome i⇤ 2 {1, . . . , d} and q 2 R. In other words,
at most S splits are allowed. A regression tree constructed this way P
is therefore a piece-wise constant
S+1function: h(x) = s=1 cs Ix2Rs , where each cs0,
S+1[s=1 Rs = X , and each Rs is a hyper-rectangle in Rd
that takes the form
8
9xik,1  qik,1 for k = 1, . . . , S1 >
>
>
>
<
=
x
> qik,2 for k = 1, . . . , S2
Rs = x 2 X : ik,2.
each qik,1 , qik,2 2 [ 1, +1] >
>
>
>
:
;
0  S1 + S2  S

Let hypothesis class H be the collection of all such regression trees. We use Theorem 4.3. Note that theaugmented class H+ = H. Since the regression treetakes constant values on each rectangle, its subgraphin the space X ⇥ R is a union of hyper-rectangles, i.e.,{(x, z) 2 X ⇥ R : h(x) > z} = [S+1s=1 Rs [ ( 1, cs ) .
Note that each Rs [ ( 1, cs ) is an intersection of at
most S + 1 axis-parallel cuts in Rd+1 , and the set of
all axis-parallel cuts is shown to have a VC dimension O(log(d)) (Gey, 2018). vc (H) can therefore be
obtained by applying VC bounds for unions and intersections of VC classes of sets (Van Der Vaart andWellner, 2009):
Theorem 4.5 (Regression tree). The class H of regression trees described as above is the same classas its augmentation H+ , and is VC-subgraph withvc (H) = vc (H+ )  CS 2 (log(S))2 log(d) for some universal constant C. If the trees are constructed in sucha way that maxx h(x) minx h(x)  M < 1, thenTheorem 4.3 can be applied with kHk 2  C 0 M foranother universal constant C 0 .
We note that, although upper bounds of VC dimension have been available for classification trees with binary features, and bounds for classification trees with
continuous features appeared very recently (Leboeuf
et al., 2020), here we consider continuous-valued regression trees with continuous features whose VC dimensions have not been addressed by previous works.
Our proof of Theorem 4.5 is based on recently developed VC results for axis-parallel cuts (Gey, 2018).
Neural Network.
Consider the class H of feedforward neural networks with a fixed architecture andfixed activation functions, indexed by the weights andbiases. Suppose the network has S 1 hidden layers,
one input layer, and two output units. Denote by W
the total number of parameters for weights and biases,
and by U the total number of computation units (neurons). Let Os 2 Rns , s = 0, . . . , S be the output of
the s-th layer. Then Os = s (Ws Os 1 + bs ) whereWs 2 Rns ⇥ns 1 is a matrix of weights, bs 2 Rns is avector of biases, and s = ( s,1 , . . . , s,ns ) is a vector of activation functions. ns denotes the number ofneurons in the s-th layer. In particular O0 = x is the
input vector and OS = (L(x), U (x)) is the final output
vector. We aim to characterize the class of one output
unit L(x) since the class of U (x) is the same.
We utilize Theorem 4.4. This approach can advantageously handle activation functions beyond sigmoid
and piece-wise polynomial (An alternate approach,which we present in Appendix D, uses Theorem 4.3
and Pollard’s pseudo-dimension (Pollard, 2012) thatapplies to sigmoid and piece-wise polynomial). Assume that each activation function s,k is globally M Lipschitz so that for some constant M0 the growth
condition | s,k (z)|  M0 + M |z| holds for all z 2 R,
and that each weight or bias parameter is restricted
to the bounded interval [ B, B] for some B > 0. To
apply Theorem 4.4, we show the following Lipschitz
property by a backpropagation-like calculation:
Theorem 4.6 (Neural network). The neural network class H = {h(·, ✓) : ✓ 2 [ B, B]W } definedabove
p the Lipschitzness
p condition
p with L(x) =p satisfies
C S(BM W )S (kxk2 + M0 U + BM W ) where C
is a universal constant. Therefore Theorem
4.4 can bepapplied
with
l
=
W
,
diam
(⇥)
=
2B
W
,
and
p
p
p
p kLk2 
C S(BM W )S (kkXk2 k2 + M0 U + BM W ).
We make two remarks. First, the sample size requiredin Theorem 4.4 to achieve a certain learning accuracy
is of order l log(CH ) log log(CH ). Applying the Lipschitz constants from Theorem 4.6 to evaluate the CH
reveals a required sample size of order W S, up to logarithmic factors, for neural networks. Second, the size
restriction B on the weights and biases enters into theerror bounds in a logarithmic manner. Therefore B
is allowed to be (exponentially) large and exerts little
impact on the training of the network.

5

Data-Driven Coverage Calibration

In this section, we propose a general-purpose PI calibration method to balance coverage and width performances in practice. On a high level, our proposal
selects the margin t in (3.2) in a data-driven manner
to guarantee a coverage maintenance.
More precisely, we recall that standard practice in validation requires 1) training models multiple times eachwith di↵erent hyperparameters, and then 2) evaluating the trained models on a validation set to selectthe optimal one. Our proposal aims to select the optimal PI model in 2). In the following, we thus assume
multiple “candidate” models are already available.
Algorithm 1 shows our procedure, which simultaneously outputs K PIs, each at a given prediction level

Haoxian Chen, Ziyi Huang, Henry Lam, Huajie Qian, Haofeng Zhang

1 ↵k (k = 1, ..., K). It starts from a candidate setof PI models, called {PIj (x) = [Lj (x), Uj (x)] : j =1, . . . , m}. These models can be obtained from setting m di↵erent values at a “tradeo↵” parameter (e.g.,
the dual multiplier in a Lagrangian formulation of theempirical constrained optimization; see Appendix G
for a neural net example), but can also be a moregeneral collection of PI models. We then use a validation data set Dv := {(Xi0 , Yi0 ) : i = 1, . . . , nv }, independent of the PI training, to check the feasibilityˆof eachP
candidate PI using the criterion CR(PIj ) :=
nv(1/nv ) i=1 IYi0 2PIj (Xi0 ) 1 ↵k +✏j for some selectedmargins ✏j .
Algorithm 1: Normalized PI Calibration
Input: Candidate PIs{PIj = [Lj , Uj ] : j = 1, . . . , m}, target coverage
rates {1 ↵k 2 (0, 1) : k = 1, . . . , K}, calibration
data Dv = {(Xi0 , Yi0 ) : i = 1, . . . , nv }, andconfidence level 1
2 (0, 1).
Procedure:
1. For each PIj , j = 1, . . . , m, compute its
empirical coverage
Pnvrate on Dv ,
1
ˆ
0
0
CR(PI
j ) := nv
i=1 IYi 2PIj (Xi ) . Compute them⇥mˆsample covariancewith
Pnv matrix ⌃ 2 R
1ˆ
00⌃j1 ,j2 = nv i=1 IYi 2PIj1 (Xi )
ˆ
ˆ
CR(PI
CR(PI
j1 ) IYi0 2PIj2 (Xi0 )j2 ) .
2ˆ
2. Let ˆj = ⌃j,j , and compute q1 , the(1)-quantile of max1jm {Zj /ˆj : ˆj > 0}
where (Z1 , . . . , Zm ) is a multivariate Gaussianˆwith mean zero and covariance ⌃.
3. For each coverage rate 1 ↵k , k = 1, . . . , K,
computemizes the interval width, indeed attains feasibility (target coverage) with 1
confidence. In this “metaoptimization”, we pool the training and validation sets
together in the objective to improve the width performance. We have the following finite-sample guarantee:
Theorem 5.1. Let 1 ↵ := maxj=1,...,m CR(PIj ),1 ↵min := 1 mink=1,...,K ↵k , and ↵˜ := min{↵min , 1
maxk=1,...,K ↵k }. For every collection of interval models {PIj : 1  j  m}, every nv , and 2 (0, 12 ), thePIs output by Algorithm 1 satisfy
PDv (CR(PIj1⇤ ↵ ) 1 ↵k for all k = 1, . . . , K)k✓⇣ 7log (mnv ) ⌘ 16
1
C1
+
nv ↵
˜
◆
✏2
exp
C2 nv min ✏,(5.1)↵(1 ↵)
with ✏ = max ↵min↵C1 (↵(1↵)/nv +

1/2)
, 0 , where PDv denotes
the probability with respect to the calibration data, andC1 , C2 are universal constants.
log(nv ↵min )/n2v ) log(m/The most important implication of Theorem 5.1 is thatthe finite-sample deterioration in the confidence level
using our procedure depends only logarithmically on m
and is independent of K. These enable both the use ofmany candidate models and the output of many PIs atdi↵erent prediction levels. The latter implies that ouralgorithm can advantageously generate all PIs for arbitrarily many target levels simultaneously with a single
validation exercise, thus is computationally cheap andcomprises a strength. We provide further interpretation on the error terms of (5.1) in Appendix E. Besides
coverage attainment guarantee in Theorem 5.1, ournncalibration procedure also possesses guaranteed perv
n 1 ⇣X⌘X
j1⇤ ↵k = arg min|PIj (Xi )| +|PIj (Xi0 )| formance regarding the optimality of the width, pron + nv i=1
1jm
i=1
vided that only the calibration data are used to assessoqˆ
the width in Step 3 of Algorithm 1. More details can1jˆ: CR(PI1 ↵k + pj)
be found in Appendix E.
nv

Output: PIj1⇤ ↵ for k = 1, . . . , K.

In addition, we provide and compare an alternate calibration scheme, viewed as an “unnormalized” (as opposed to “normalized”) version of Algorithm 1 whenhandling the standard error ˆj , in Appendix F.

The key of our procedure is to tune ✏j based on a
uniform central limit theorem (CLT) that captures
the overall errors incurred in the empirical coverage
rates. Denoting by CR(PIj ) := P⇡ (Y 2 PIj (X)) thetrue coverage rate of PIj , this CLT implies that, setting q1= (1)-quantile of maxj Zj /ˆj for someproperly chosen Gaussian vector (Zj )j=1,...,m , we havepˆCR(PIj ) CR(PIˆj / nv for all j = 1, . . . , mj ) q1
uniformly with probability ⇡ 1. The uniformity
over j ensures the solution in Step 3, which opti-

We discuss Algorithm 1 in relation to a naive, butnatural approach that simply selects the model withthe shortest average interval width, among candidate
models with empirical coverage rates on the validation
dataset reaching the target levels (i.e., t = 0 in (3.2)).
This latter approach is also known as the PAV validation scheme in Kivaranovic et al. (2020) (which we call
NNVA in Section 6). Our algorithms improve PAV in
two aspects. First, our proposal guarantees with highprobability that the target level will be achieved by thetest coverage thanks to a corrective margin, while PAV

where |PIj (·)| := Uj (·) Lj (·) is the width, and{Xi }ni=1 is the training data set.
k

Learning Prediction Intervals for Regression: Generalization and Calibration

does not o↵er such a guarantee and tends to fall short.
Second, our proposal enjoys a higher statistical power
than PAV in that unlike PAV whose analysis is basedon concentration bounds, Algorithm 1 is analyzed via
an asymptotically tight joint CLT. These guarantees
build on recent high-dimensional Berry-Esseen bounds(Chernozhukov et al., 2017) (which notably does not
require functional complexity measures but only thegeometry of “hit sets”). Moreover, we will observe inthe experiments in Section 6 that our proposal empirically performs better than PAV.

6

Experiments

Datasets. We evaluate our approaches on both synthetic datasets and real-world benchmark datasets
through comparisons to the state-of-the-arts. The
real-world datasets (“Boston”, “Concrete”, “Wine”,“Energy” and “Yacht”) have been widely used in previous studies (Hernández-Lobato and Adams, 2015;
Gal and Ghahramani, 2016; Lakshminarayanan et al.,2017) for regression tasks. The generative distributions for the three synthetic datasets are:TT2(1) : f (x) = c 2x + 10 sin( c 8x ) + ||x||10 ✏, x ⇠ N (0, I10 ),||x||2
1 T
2T(2) : f (x) = 8 (c x) sin(c x) + 10 ✏, x ⇠ N (0, I7 ),2(3) : f (x) = 12 cT x · cos(cT x)2 + ||x||10 ✏, x ⇠ N (0, I9 ).

where ✏ ⇠ N (0, 1), and c is a constant in [ 2, 2]10 ,[ 2, 2]7 , [ 2, 2]9 respectively.
Experimental Setup. We conduct experiments under two scenarios: the single PI case, where one PI at
a single prediction level 1 ↵ are constructed, and thesimultaneous PI case, where K PIs at K di↵erent
prediction levels 1 ↵1 , · · · , 1 ↵K are constructed.
Each trial is repeated for N times to estimate the confidence of coverage attainment. We adopt neural networks as our PI models with the following loss function: l(x, y, L, U ) := (U (x) L(x))2 + (max{L(x)y, 0} + max{y U (x), 0})2 , where > 0 is a penalty
for miscoverage and (L(x), U (x)) is the output vector
containing the lower and upper bounds. By adjusting , PI models with di↵erent coverage levels can betrained, which are then fed into our calibration algorithms to obtain the final PI. We implement two calibration strategies: the normalized Gaussian PI calibration in Algorithm 1 (NNGN), and the unnormalized version in Algorithm 2 in Appendix F (NNGU). In
addition, we also test the calibration scheme (NNVA)
that directly compares the empirical coverage rates on
the validation dataset to the target levels, without theGaussian margin. Note that this is the PAV validation
scheme in Kivaranovic et al. (2020).
Baselines. We compare our NNGN, NNGU with thefollowing state-of-art approaches: quantile regression

forests (QRF) (Meinshausen, 2006), CV+ prediction
interval (CV+) (Barber et al., 2019), split conformalized quantile regression (SCQR) (Romano et al.,2019), quantile regression via SVM (SVMQR) (Steinwart and Thomann, 2017), and split conformal learning (SCL) (Lei et al., 2018). Implementation details
can be found in Appendix I.
Evaluation Metrics. For the single PI analysis, ourmodels are evaluated on both exceedance probability(EP ) and interval width (IW ). EP captures the success in achieving the target confidence level, while IW
indicates the average interval width. For the simultaneous PI analysis, we use multiple exceedance probability (M EP ) and multiple interval width (M IW ).
M EP measures the proportion of trials where all PIs
reach the target prediction levels simultaneously (i.e.,family-wise correctness). Formally:
EP :
IW :
M EP :M IW :PN
1
P L}
NP i=1P {CRiNn1
(Ui,j Li,j )i=1
NnP
TK j=1N
1{{CR
P Lk }}i,k
N Pi=1 P k=1
K
N Pn
1
(U
Li,j,k )
i,j,k
k=1
i=1
j=1
KN n

where n is the size of testing data, CRi (CRi,k ) is
the estimated coverage rate from the i-th repetition
and P L (P Lk ) is the target prediction level 1 ↵ (1↵k ). Throughout our experiments, the confidence level1
is set to 0.9 in the calibration algorithms. For
both cases, the best result is achieved by the modelwith the smallest IW/M IW value among those with
EP/M EP
0.90. If no model achieves EP/M EP
0.90, then the one with highest EP/M EP is the best.
Single PI Analysis. Table 1 reports the values ofEP and IW for PI generation at 95% prediction levelon 3 synthetic datasets and 5 real-world datasets. It
is shown that QRF, CV+ and our NNGU and NNGN
are the only four methods that achieve the required
confidence level in all synthetic cases, and among thefour our NNGN consistently generates PIs of shortest
width. Moreover, NNGU attains the target confidence
level on all real datasets as well. NNGN seems to fall
below the target confidence for some real datasets, but
is better than all other baseline methods except CV+
and SCL. Among the methods with high EP , the interval widths of our NNGN and NNGU are the smallest
in 6 out of 8 datasets. In contrast, the EP values by
NNVA are below the target confidence in all cases.
Numerically, the averaged EP of NNGN/NNGU is
1.4/1.7 times higher than the one of NNVA. This showsin particular that NNVA may fail to ensure a correctcoverage rate with high confidence on test data, which
necessitates our calibration approaches.
Simultaneous PIs Analysis. Table 2 reports thevalues of M EP and M IW for simultaneous PIs at 19

Haoxian Chen, Ziyi Huang, Henry Lam, Huajie Qian, Haofeng Zhang

Table 1: Single PI at the 95% prediction level. The best results are in bold.
Synthetic1 Synthetic2 Synthetic3
BostonConcrete
EnergyWine
Yacht
Methods
EP IW EP IW EP IW EP IW EP IW EP IW EP IW EP IW
QRF
1.00 3.122 1.00 3.667 1.00 3.609 0.46 2.085 0.72 1.995 0.78 0.664 0.00 2.671 0.22 1.110
CV+
1.00 0.302 1.00 2.039 1.00 3.523 0.96 1.538 0.88 1.413 0.98 0.500 1.00 4.359 0.92 0.339
SCQR
0.90 2.264 0.78 2.863 0.98 2.804 0.66 2.309 0.62 2.289 0.68 0.837 0.30 2.810 0.86 1.681
SVMQR 0.00 0.256 0.00 2.351 0.00 1.855 0.06 1.340 0.10 1.376 0.18 0.411 0.48 2.975 0.30 0.483
SCL
0.98 0.313 0.70 2.211 1.00 3.754 0.88 2.053 0.88 1.630 0.74 0.769 0.80 4.437 0.92 0.906
NNVA
0.00 0.221 0.74 1.630 0.04 1.991 0.58 1.837 0.28 1.607 0.44 0.234 0.00 1.817 0.80 0.154
Ours-NNGN 1.00 0.296 1.00 1.921 1.00 2.176 0.90 2.477 0.86 2.375 0.62 0.398 0.74 2.155 0.92 0.217
Ours-NNGU 1.00 0.296 1.00 2.557 1.00 3.155 0.96 2.692 0.96 2.643 1.00 0.561 0.98 2.648 1.00 0.299

Table 2: Simultaneous PIs at 19 target prediction levels. The best results are in bold.
Synthetic1 Synthetic2 Synthetic3
BostonConcrete
EnergyWine
Yacht
Methods MEP MIW MEP MIW MEP MIW MEP MIW MEP MIW MEP MIW MEP MIW MEP MIW
QRF
1.00 1.911 0.92 1.657 1.00 1.805 0.46 1.099 0.72 1.108 0.78 0.328 0.00 1.460 0.22 0.703
CV+
1.00 0.176 1.00 1.078 1.00 1.935 0.66 0.780 0.60 0.736 1.00 0.245 1.00 2.289 0.86 0.114
SCQR
0.74 1.356 0.74 2.218 0.56 1.808 0.14 1.622 0.16 1.575 0.24 0.760 0.00 2.493 0.20 1.554
SVMQR
0.00 0.115 0.00 1.311 0.00 1.215 0.00 0.582 0.00 0.619 0.00 0.225 0.04 1.575 0.00 0.086
SCL
0.24 0.175 1.00 1.004 0.96 1.940 0.78 1.039 0.66 0.896 0.74 0.373 0.70 2.470 0.80 0.325
NNVA
0.00 0.160 0.28 0.969 0.04 1.606 0.26 1.086 0.16 0.956 0.00 0.151 0.00 1.147 0.04 0.079
Ours-NNGN 1.00 0.170 1.00 1.050 1.00 1.727 0.76 1.244 0.72 1.092 0.90 0.177 0.72 1.267 0.96 0.120
Ours-NNGU 1.00 0.168 1.00 1.083 1.00 1.734 0.82 1.287 0.78 1.123 0.60 0.183 0.98 1.276 0.72 0.145

target prediction levels: 50%, 52.5%, 55%, 57.5% . . .,95%. Our calibration approaches NNGN and NNGU
are always the best in terms of achieving the tightest interval under high M EP , or otherwise achieves
the highest M EP among all. Among the 6 datasets
where the target confidence level 0.9 can be attained,
NNGN/NNGU yields the smallest width in 4/2 ofthem. In the remaining 2 datasets, NNGU achieves
the highest M EP . NNGU attains the target M EP or
the highest M EP in 6 out of 8 datasets. Compared tothe case of single-level target, the M EP performancegaps are more significant in multi-level PI constructions. This is because the coverage rates in baseline
algorithms get increasingly overfitted as more simultaneous target levels are compared against. Thanks to
the use of the uniform safety margin, our NNGN and
NNGU schemes are free of overfitting even in this case.
These results demonstrate that our methods can accurately construct multiple PIs at di↵erent predictionlevels simultaneously. Finally, compared to NNGU,
NNGN tends to generate shorter PIs, and we recommend NNGN as the preferred choice.

We develop a general learning theory to characterize
the optimality-feasibility tradeo↵ in this optimization,in particular joint guarantees on both a short expectedinterval width and an attainment of the target prediction level. We also propose a readily implementablecalibration procedure, constructed based on a highdimensional Berry-Esseen Theorem, to select the best
PI model among trained candidates, which o↵ers apractical approach to build simultaneous PIs at multiple target prediction levels with statistical validity. We
demonstrate the empirical strengths of our proposed
approach by applying it to neural-network-based PImodels with our proposed calibration procedure, andcomparing them with other baselines across synthetic
and real-data examples.

We provide more experimental results in Appendix I.

References

7A. M. Alaa and M. van der Schaar. Discriminative
jackknife: Quantifying uncertainty in deep learning
via higher-order influence functions. arXiv preprint
arXiv:2007.13481, 2020.

Conclusion

In this paper, we study the generation of PIs for regression that satisfy an expected coverage rate. This
problem can be cast into an empirical constrained optimization framework that minimizes the expected interval width subject to a coverage satisfaction constraint.

Acknowledgments
We gratefully acknowledge support from the National
Science Foundation under grants CAREER CMMI1834710 and IIS-1849280.

C. Anil, J. Lucas, and R. Grosse. Sorting out lipschitz
function approximation. In International Conference on Machine Learning, pages 291–301, 2019.
B. Ankenman, B. L. Nelson, and J. Staum. Stochas-Learning Prediction Intervals for Regression: Generalization and Calibration

tic kriging for simulation metamodeling. Operations
Research, 58(2):371–382, 2010.
R. F. Barber, E. J. Candes, A. Ramdas, and R. J.
Tibshirani. Predictive inference with the jackknife+.
arXiv preprint arXiv:1905.02928, 2019.
P. L. Bartlett, D. J. Foster, and M. J. Telgarsky.
Spectrally-normalized margin bounds for neural networks. In Advances in Neural Information Processing Systems, pages 6240–6249, 2017.
G. E. Box and G. C. Tiao. Bayesian inference in statistical analysis, volume 40. John Wiley & Sons,
2011.
V. Chernozhukov, D. Chetverikov, K. Kato, et al. Central limit theorems and bootstrap in high dimensions. The Annals of Probability, 45(4):2309–2352,
2017.
M. Cisse, P. Bojanowski, E. Grave, Y. Dauphin, and
N. Usunier. Parseval networks: Improving robustness to adversarial examples. In Proceedingsof the 34th International Conference on Machine
Learning-Volume 70, pages 854–863. JMLR. org,
2017.
K. Doksum and J.-Y. Koo. On spline estimatorsand prediction intervals in nonparametric regression. Computational Statistics & Data Analysis, 35(1):67–82, 2000.
R. Dudley. Universal donsker classes and metric entropy. The Annals of Probability, pages 1306–1326,
1987.

J. M. Hernández-Lobato and R. Adams. Probabilistic backpropagation for scalable learning of bayesianneural networks. In International Conference on
Machine Learning, pages 1861–1869, 2015.
T. Huster, C.-Y. J. Chiang, and R. Chadha. Limitations of the lipschitz constant as a defense against
adversarial examples. In Joint European Conference
on Machine Learning and Knowledge Discovery in
Databases, pages 16–29. Springer, 2018.
A. Khosravi, S. Nahavandi, D. Creighton, and A. F.
Atiya. Lower upper bound estimation method forconstruction of neural network-based prediction intervals. IEEE Transactions on Neural Networks, 22(3):337–346, 2010.
A. Khosravi, S. Nahavandi, D. Creighton, and A. F.
Atiya. Comprehensive review of neural networkbased prediction intervals and new advances. IEEE
Transactions on Neural Networks, 22(9):1341–1356,
2011.
B. Kim, C. Xu, and R. F. Barber. Predictive inferenceis free with the jackknife+-after-bootstrap. arXiv
preprint arXiv:2002.09025, 2020.
D. Kivaranovic, K. D. Johnson, and H. Leeb. Adaptive, distribution-free prediction intervals for deep
networks. In International Conference on Artificial
Intelligence and Statistics, pages 4346–4356, 2020.
R. Koenker and K. F. Hallock. Quantile regression.
Journal of Economic Perspectives, 15(4):143–156,
2001.

Y. Gal and Z. Ghahramani. Dropout as a bayesian
approximation: Representing model uncertainty indeep learning. In International Conference on Machine Learning, pages 1050–1059, 2016.

B. Lakshminarayanan, A. Pritzel, and C. Blundell.
Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in neural
information processing systems, pages 6402–6413,
2017.

I. M. Galván, J. M. Valls, A. Cervantes, and R. Aler.
Multi-objective evolutionary optimization of prediction intervals for solar energy forecasting with neuralnetworks. Information Sciences, 418:363–382, 2017.

J.-S. Leboeuf, F. LeBlanc, and M. Marchand. Decision trees as partitioning machines to characterize their generalization properties. arXiv preprint
arXiv:2010.07374, 2020.

S. Gey.
Vapnik–chervonenkis dimension of axisparallel cuts. Communications in Statistics-Theory
and Methods, 47(9):2291–2296, 2018.

J. Lei, A. Rinaldo, and L. Wasserman. A conformal
prediction approach to explore functional data. Annals of Mathematics and Artificial Intelligence, 74(1-2):29–43, 2015.

H. Gouk, E. Frank, B. Pfahringer, and M. Cree. Regularisation of neural networks by enforcing lipschitz
continuity. arXiv preprint arXiv:1804.04368, 2018.
C. Gupta, A. K. Kuchibhotla, and A. K. Ramdas.
Nested conformal prediction and quantile out-of-bag ensemble methods. arXiv preprint
arXiv:1910.10562, 2019.
M. Hein and M. Andriushchenko. Formal guarantees
on the robustness of a classifier against adversarial
manipulation. In Advances in Neural Information
Processing Systems, pages 2266–2276, 2017.

J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, andL. Wasserman. Distribution-free predictive inference
for regression. Journal of the American Statistical
Association, 113(523):1094–1111, 2018.
N. Meinshausen. Quantile regression forests. Journal of Machine Learning Research, 7(Jun):983–999,
2006.
D. J. Olive. Prediction intervals for regression models.
Computational Statistics and Data Analysis, 51(6):3115–3122, 2007.

Haoxian Chen, Ziyi Huang, Henry Lam, Huajie Qian, Haofeng Zhang

G. Papadopoulos, P. J. Edwards, and A. F. Murray.
Confidence estimation methods for neural networks:
A practical comparison. IEEE Transactions on Neural Networks, 12(6):1278–1287, 2001.
H. Papadopoulos. Inductive conformal prediction:
Theory and application to neural networks. In Tools
in artificial intelligence. Citeseer, 2008.
T. Pearce, M. Zaki, A. Brintrup, and A. Neely. Highquality prediction intervals for deep learning: A
distribution-free, ensembled approach. In International Conference on Machine Learning, PMLR:
Volume 80, 2018.
D. Pollard.
Convergence of stochastic processes.
Springer Science and Business Media, 2012.
W. Polonik. Minimum volume sets and generalized
quantile processes. Stochastic Processes and Their
Applications, 69(1):1–24, 1997.
Y. Romano, E. Patterson, and E. Candes. Conformalized quantile regression. In Advances in Neural
Information Processing Systems, pages 3543–3553,
2019.
N. Rosenfeld, Y. Mansour, and E. Yom-Tov. Discriminative learning of prediction intervals. In International Conference on Artificial Intelligence andStatistics, pages 347–355, 2018.
J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn.
Design and analysis of computer experiments. Statistical Science, pages 409–423, 1989.
R. L. Schmoyer. Asymptotically valid prediction intervals for linear models. Technometrics, 34(4):399–
408, 1992.
C. D. Scott and R. D. Nowak. Learning minimum
volume sets. Journal of Machine Learning Research,
7(Apr):665–704, 2006.
G. A. Seber and A. J. Lee. Linear regression analysis,volume 329. John Wiley & Sons, 2012.
L. Steinberger and H. Leeb. Leave-one-out prediction
intervals in linear regression models with many variables. arXiv preprint arXiv:1602.05801, 2016.
I. Steinwart and P. Thomann. liquidsvm: A fast andversatile svm package, 2017.
I. Steinwart, A. Christmann, et al. Estimating conditional quantiles with the help of the pinball loss.
Bernoulli, 17(1):211–225, 2011.
R. A. Stine. Bootstrap prediction intervals for regression. Journal of the American Statistical Association, 80(392):1026–1031, 1985.
Y. Tsuzuku, I. Sato, and M. Sugiyama. Lipschitzmargin training: Scalable certification of perturbation invariance for deep neural networks. In Advances in Neural Information Processing Systems,
pages 6541–6550, 2018.

A. Van Der Vaart and J. A. Wellner. A note on bounds
for vc dimensions. Institute of Mathematical Statistics Collections, 5:103, 2009.
A. W. Van der Vaart and J. A. Wellner. Weak Convergence and Empirical Processes with Applications
to Statistics. Springer, 1996.
V. Vapnik. The nature of statistical learning theory.
Springer Science and Business Media, 2013.
V. Vovk. Conditional validity of inductive conformal
predictors. In Asian Conference on Machine Learning, pages 475–490, 2012.
V. Vovk. Cross-conformal predictors. Annals of
Mathematics and Artificial Intelligence, 74(1-2):9–28, 2015.
V. Vovk, A. Gammerman, and G. Shafer. Algorithmic learning in a random world. Springer Science &Business Media, 2005.
Y. Yoshida and T. Miyato. Spectral norm regularization for improving the generalizability of deep learning. arXiv preprint arXiv:1705.10941, 2017.
H. Zhang, J. Zimmerman, D. Nettleton, and D. J.
Nordman. Random forest prediction intervals. The
American Statistician, pages 1–15, 2019.
L. Zhu, J. Lu, and Y. Chen. HDI-forest: highest density interval regression forest. In Proceedings of the
28th International Joint Conference on Artificial Intelligence, pages 4468–4474. AAAI Press, 2019.