Year: 2020

Abstract:
Deep learning has achieved state-of-the-art performance to generate high-quality prediction intervals (PIs) for uncertainty quantification in regression tasks. The high-quality criterion requires PIs to be as narrow as possible, whilst maintaining a pre-specified level of data (marginal) coverage. However, most existing works for high-quality PIs lack accurate information on conditional coverage, which may cause unreliable predictions if it is significantly smaller than the marginal coverage. To address this problem, we propose a novel end-to-end framework which could output high-quality PIs and simultaneously provide their conditional coverage estimation. In doing so, we design a new loss function that is both easyto-implement and theoretically justified via an exponential concentration bound. Our evaluation on real-world benchmark datasets and synthetic examples shows that our approach not only outperforms the state-of-the-arts on high-quality PIs in terms of average PI width, but also accurately estimates conditional coverage information that is useful in assessing model uncertainty.

text: Deep learning has achieved state-of-the-art performance to generate high-quality prediction intervals (PIs) for uncertainty quantification in regression tasks. The high-quality criterion requires PIs to be as narrow as possible, whilst maintaining a pre-specified level of data (marginal) coverage. However, most existing works for high-quality PIs lack accurate information on conditional coverage, which may cause unreliable predictions if it is significantly smaller than the marginal coverage. To address this problem, we propose a novel end-to-end framework which could output high-quality PIs and simultaneously provide their conditional coverage estimation. In doing so, we design a new loss function that is both easyto-implement and theoretically justified via an exponential concentration bound. Our evaluation on real-world benchmark datasets and synthetic examples shows that our approach not only outperforms the state-of-the-arts on high-quality PIs in terms of average PI width, but also accurately estimates conditional coverage information that is useful in assessing model uncertainty.

heading: 1 INTRODUCTION
text: Prediction interval (PI) is poised to play an increasingly prominent role in uncertainty quantification for regression tasks (Khosravi et al., 2010; 2011; Galván et al., 2017; Rosenfeld et al., 2018; Tagasovska & Lopez-Paz, 2018; 2019; Romano et al., 2019; Wang et al., 2019; Kivaranovic et al., 2020). A high-quality PI should be as narrow as possible, whilst maintaining a pre-specified level of data coverage or marginal coverage (Pearce et al., 2018). Compared with PIs obtained based on coverage-only consideration, the “high-quality” criterion is beneficial in balancing between marginal coverage probability and interval width. However, the conditional coverage given a feature, which is critical for making reliable context-based decisions, is unassessed and missing in most existing works on high-quality PIs. In the presence of heteroskedasticity and model misspecification, the marginal coverage can be very different from the conditional coverage at a given point, which affects the downstream decision-making task that relies on the uncertainty information provided by the PI. Our main goal is to meaningfully incorporate and assess conditional coverages in high-quality PIs.
Conditional coverage estimation is challenging for two reasons. First is that the natural evaluation metric of conditional coverage error, an Lp distance between the estimated and ground-truth conditional coverages, is difficult to compute as it requires obtaining the conditional probability given feature x, which is arguably as challenging as the regression problem itself. Our first goal in this paper is to address this issue by developing a new metric called calibration-based conditional coverage error for conditional coverage estimation measurement. Our approach is inspired from the calibration notion in classification (Guo et al., 2017). The basic idea is to relax conditional coverage at any given point to being averaged over all points that bear the same estimated value. An estimator satisfying the relaxed property is regarded as well-calibrated. In regression, calibration-based conditional coverage error provides a middle ground between the enforcement of marginal coverage (lacking any conditional information) and conditional coverage (computationally intractable). Compared with conditional coverage, this middle-ground metric can be viewed as a “dimension reduction” of the conditioning variable from the original sample space to the space [0, 1], so that we can easily discretize to compute the empirical metric values.
The second challenge is the discontinuity in the above metrics that hinders efficient training of PIs that are both high-quality and possess reliable conditional coverage information. To address this, we design a new loss function based on a combination of the high-quality criterion and a coverage assessment loss. The latter can be flexibly added as a separate module to any neural network (NN) used to train PIs. It is based on an empirical version of a tight upper bound on the coverage error in terms of a Kullback–Leibler (KL) divergence, which can be readily employed for running gradient descent. We theoretically show how training with our proposed loss function attains this upperbounding value via a concentration bound. We also demonstrate the empirical performance of our approach in terms of PI quality and conditional coverage assessment compared with benchmark methods.
Summary of Contributions: (1) We identify the conditional coverage estimation problem as a new challenge for high-quality PIs and introduce a new evaluation metric for coverage estimation. (2) We propose an end-to-end algorithm that can simultaneously construct high-quality PIs and generate conditional coverage estimates. In addition, we provide theoretical justifications on the effectiveness of our algorithm by developing concentration bounds relating the coverage assessment loss and conditional coverage error. (3) By evaluating on benchmark datasets and synthetic examples, we empirically demonstrate that our approach not only achieves high performance on conditional coverage estimation, but also outperforms the state-of-the-art algorithms on high-quality PI generation.

heading: 2 EVALUATING CONDITIONAL COVERAGE FOR HIGH-QUALITY PIS
text: Let X ∈ X and Y ∈ Y ⊂ R be random variables denoting the input feature and label, where the pair (X,Y ) follows an (unknown) ground-truth joint distribution π(X,Y ). Let π(Y |X) be the conditional distribution of Y given X . We are given the training data D := {(xi, yi), i = 1, 2, · · · , n} where (xi, yi) are i.i.d. realizations of random variables (X,Y ). A PI refers to an interval [L(x), U(x)] where L, U are two functions mapping from X to Y trained on the data D. [L(x), U(x)] is called a PI at prediction level 1− α (0 ≤ α ≤ 1) if its marginal coverage is not less than 1 − α, i.e., P[Y ∈ [L(X), U(X)]|L,U ] ≥ 1 − α where P is with respect to a new test point (X,Y ) ∼ π. We say that [L(x), U(x)] is of high-quality if its marginal coverage attains a pre-specified target prediction level and has a short width on average. In particular, a best-quality PI at prediction level 1− α is an optimal solution to the following constrained optimization problem:
min L,U
E[U(X)− L(X)] subject to P[Y ∈ [L(X), U(X)]|L,U ] ≥ 1− α. (2.1)
The high-quality criterion has been widely adopted in previous work (see Section 6). However, this criterion alone may fail to carry important model uncertainty information at specific test points. Consider a simple example where x ∼ Uniform[0, 1], y = 0 for x ∈ [0, 0.95] and y|x ∼ Uniform[0, 1] for x ∈ (0.95, 1]. Then according to equation 2.1, a best-quality 95% PI is precisely L(x) = U(x) = 0 for all x ∈ [0, 1]. This PI has nonconstant coverage if we condition at different points (1 for x ∈ [0, 0.95] and 0 for x ∈ (0.95, 1]), and can deviate significantly from the overall coverage 95%. More examples to highlight the need of obtaining conditional coverage information can be found in our numerical experiments in Section 5.1.
To mitigate the drawback of the high-quality criterion, we define: Definition 2.1 (Conditional Coverage and Its Estimator). The conditional coverage associated with a PI [L(x), U(x)] is A(x) := P[Y ∈ [L(X), U(X)]|L,U,X = x] for a.e. x ∈ X , where P is taken with respect to π(Y |X). For a (conditional) coverage estimator P̂ , which is a measurable function from X to [0, 1], we define its Lp conditional coverage error (C̃Ep) as
C̃Ep := ∥∥∥P[Y ∈ [L(X), U(X)]|L,U,X]− P̂ (X)∥∥∥
Lp(X )
where the Lp-norm is taken with respect to the randomness of X (1 ≤ p ≤ +∞).
Note that evaluating C̃Ep relies on approximating the conditional coverage A(x), which can be as challenging as the original prediction problem. To address this, we leverage the similarity of estimatingA(x) to generating prediction probabilities in binary classification, which motivates us to
borrow the notion of calibration in classification. This idea is based on a relaxed error criterion by looking at the conditional coverage among all points that bear the same coverage estimator value, instead of conditioning at any given point. The resulting error metric then only relies on probabilities conditioned on variables in a much lower-dimensional space [0, 1] than X . To explain concretely, we introduce a “perfect-calibrated coverage estimator” as:
Definition 2.2 (Perfect Calibration). A coverage estimator P̂ is called a perfect-calibrated coverage estimator associated with [L(x), U(x)] if it satisfies
P̂ (x) = P[Y ∈ [L(X), U(X)]|L,U, P̂ (X) = P̂ (x)], a.e. P̂ (x) ∈ [0, 1]. (2.2)
where a.e. is with respect to the probability measure on [0, 1] induced by the random variable P̂ (X).
Equation 2.2 means that a point x with conditional coverage estimate P̂ (x) = p has an average coverage of precisely p, among all points in X that possess the same conditional coverage estimated value. That is, the average coverage of the PI restricted on the subset {x ∈ X : P̂ (x) = p} should be precisely p. Corresponding to Definition 2.2, we define: Definition 2.3 (Calibration-based Error). An Lp (1 ≤ p ≤ +∞) calibration-based conditional coverage error, or coverage error for short (CEp), of a coverage estimator P̂ is:
CEp := ∥∥∥P[Y ∈ [L(X), U(X)]|L,U, P̂ (X)]− P̂ (X)∥∥∥
Lp(X ) (2.3)
where Lp-norm is taken with respect to the randomness of P̂ (X).
In the above definition the conditional probability P[Y ∈ [L(X), U(X)]|L,U, P̂ (X)] is a measurable function of random variable P̂ (X), say γ(P̂ (X)). By a change of variable,
CEpp := ∥∥∥γ(P̂ (X))− P̂ (X)∥∥∥p
Lp(X ) = ∫ 1 0 |γ(t)− t|pdFP̂ (X)(t) (2.4)
where FP̂ (X) is a probability distribution of P̂ (X) on [0, 1]. Here, CEp only requires estimating γ(t) for t ∈ [0, 1], which can be done easily by discretizing [0, 1] for empirical calculation. We call the empirical Lp calibration-based conditional coverage error ECEp. More details of this empirical error can be found in Appendix A.4.
A calibration-based error CEp provides a middle ground between the enforcement of marginal coverage and conditional coverage. The ground-truth conditional coverage is perfectly calibrated, but not vice versa. However, if we enforce the perfect calibration criterion for a coverage estimator to hold when restricted to any measurable subset in X , then the choice of estimator will reduce uniquely to the conditional coverage (Definition A.1 and Lemma A.2). In this sense, CEp is an error metric that is a natural relaxation of C̃Ep, and although less precise, CEp is computationally much more tractable than C̃Ep.
Evaluation Metric for Coverage Estimator. We use ECE1 as the primary evaluation metric to measure the quality of a coverage estimator. A high ECE1 value of a coverage estimator indicates an unreliable coverage estimation while a small ECE1 value indicates that the coverage estimator is close to the perfect-calibrated property. Ideally, an effective algorithm should output a coverage estimator with a small ECE1 value.

heading: 3 METHODOLOGY
text: We propose a novel end-to-end algorithm, named coverage assessment network (CaNet), to simultaneously generate a coverage estimator along with the high-quality PI. As illustrated in Figure 1, our CaNet consists of two major modules: (1) predictor module and (2) coverage assessment module (Ca-Module). The predictor module provides the upper and lower bound of the estimated PIs. Meanwhile, the Ca-Module is added to the output layer to access the conditional coverage information of PIs from the predictor module. Our model is jointly optimized by three loss functions: coverage assessment loss LCA, intervals width loss LIW , and coverage probability loss LCP . Benefited from these powerful modules, the CaNet can generate and validate the coverage estimator at the same time without any requirement for further post-processing steps.
Coverage Assessment Module. Our Ca-Module consists of two neurons fully connected to the last hidden layer to estimate the conditional coverage. After passing through the softmax activation function, it outputs a two-point probability distribution P̂ (x), 1− P̂ (x) where P̂ (x) is the coverage estimator of the PI from the predictor module. Our Ca-Module can be easily integrated into the output layer of deep neural networks to estimate their conditional coverage.
Loss Function Design and Tuning Procedure. Our loss function is a sum of the predictor loss and the coverage loss. The predictor loss aims to narrow the prediction intervals as much as possible, while maintaining a specified marginal coverage of data. Inspired by (Khosravi et al., 2010; 2011; Pearce et al., 2018; Rosenfeld et al., 2018), our predictor loss is formed by the sum of intervals width (IW) loss LIW and coverage probability (CP) loss LCP :
LIW = 1
n n∑ i=1 (U(xi)− L(xi)), (3.1)
LCP = 1
n n∑ i=1 k̃i, CP = 1 n n∑ i=1 ki, (3.2)
where ki indicates whether each data point has been captured by the PIs: ki = 1 if L(xi) ≤ yi ≤ U(xi) and ki = 0 otherwise. k̃i is a soft version of ki, which is defined as: k̃i := σ(λ3(U(xi)−yi))· σ(λ3(yi−L(xi))), where λ3 ≥ 0 is a tunable parameter and σ(t) := 11+e−t is the sigmoid function. Therefore, LCP is a soft version of CP that can be used for gradient descent. Associated with the Ca-Module, we introduce a coverage assessment loss LCA to estimate the conditional coverage:
LCA = 1
n n∑ i=1 ( ki log(P̂ (xi)) + (1− ki) log(1− P̂ (xi)) ) (3.3)
We will show in Section 4 that the expectation of coverage assessment loss LCA provides an upper bound for both the conditional coverage error (Definition 2.1) and calibration-based conditional coverage error (Definition 2.3). Hence, minimizing LCA contributes to the recovery of the conditional coverage. In order to run gradient-based methods, we replace the discrete indicator (ki, 1 − ki) in LCA with its soft version (k̃i, 1− k̃i):
L̃CA = 1
n n∑ i=1 ( k̃i log(P̂ (xi)) + (1− k̃i) log(1− P̂ (xi)) ) (3.4)
Our total loss function for the CaNet is defined as:
Total Loss = LIW + λ1(1− LCP ) + λ2L̃CA (3.5)
where λ1 ≥ 0, λ2 ≥ 0 are tunable parameters. We propose an easy-to-implement yet effective tuning procedure to pick up these parameters. Please refer to Appendix D for more algorithm details.
Deep Ensembles. Following previous research (Lee et al., 2015; Lakshminarayanan et al., 2017; Pearce et al., 2018; Fort et al., 2019; Ovadia et al., 2019; Gustafsson et al., 2020; Pearce et al., 2020), we apply the deep ensemble technique to provide more robust and better results. During the training period, with the same hyperparameters λi, i = 1, 2, 3, m networks are trained with different
initializations. The prediction results from i-th network are denoted as: ([Li(x), Ui(x)], P̂i(x), 1 − P̂i(x)). Finally, the output from CaNet is:
Lower bound L̄ := ∑m i=1 1 mLi,
Upper bound Ū := ∑m i=1 1 mUi,
Coverage estimator ¯̂ P := ∑m i=1 1 m P̂i.
(3.6)

heading: 4 THEORETICAL ANALYSIS
text: In this section, we provide theoretical insights that minimizing LCA is equivalent to minimizing a tight upper bound of the conditional coverage error and thus can recover the true conditional coverage. To achieve this, we first prove that both CEp and C̃Ep are bounded above by the expectation of a Kullback–Leibler divergence-type random variable K1(X). Then we show that LCA is an empirical counterpart of K1(X) and establish a concentration bound between LCA and E[K1(X)]. Theorem 4.1. Let A(x) := P[Y ∈ [L(X), U(X)]|L,U,X = x] be the conditional coverage in Definition 2.1. Let K(x) = A(x) log ( A(x)
P̂ (x)
) + (1−A(x)) log ( 1−A(x) 1−P̂ (x) ) . Then
CEp ≤ C̃Ep ≤ ( 1
2 E[K(X)]
)αp/2 ,∀1 ≤ p ≤ +∞
where αp = 1,∀1 ≤ p ≤ 2 and αp = 2p ,∀2 ≤ p ≤ +∞. Moreover, the inequality is attainable if, e.g., P̂ (x) equals the conditional coverage A(x).
From Theorem 4.1, we see that minimizing E[K(x)] is equivalent to minimizing a tight upper bound for the coverage error. For every x, K(x) is the Kullback–Leibler divergence between the distributions represented by (P̂ (x), 1 − P̂ (x)) and (A(x), 1 − A(x)). K(x) = K0(x) + K1(x), where K0(x) = A(x) log(A(x)) + (1 − A(x)) log(1 − A(x)) and K1(x) = −A(x) log(P̂ (x)) − (1 − A(x)) log(1 − P̂ (x)). Minimizing E[K(X)] over P̂ is equivalent to minimizing E[K1(X)]. The type of results in Theorem 4.1 that bounds an Lp conditional coverage error via a Kullback-Leiblertype error is new as far as we know. Next, to show LCA approximates E[K1(X)], we need the following assumptions: Assumption 4.2. The four classes of functions ([L(x), U(x)], P̂ (x), 1− P̂ (x)) output by the neural network (NN) in Figure 1 have finite VC dimensions, say they are bounded above by V0.
Assumption 4.2 holds for a wide range of NNs (e.g., Theorem 8.14 in Anthony & Bartlett (2009), Theorem 7 in Bartlett et al. (2019)). In particular, it holds for the one we adopt in the experiments (where we use the ReLU-activated NN to construct ψi, i = 1, 2, 3, 4; see Section 5): Theorem 4.3. Suppose ψi, i = 1, 2, 3, 4 are the pre-activated output neurons of the NN in Figure 1 using the ReLU activation function. Then Assumption 4.2 holds. Moreover, suppose the NN has W parameters and U computation units (nodes). Then V0 = O(WU).
Assumption 4.4. | log(P̂ (x))| ≤M , | log(1− P̂ (x))| ≤M for all x and P̂ .
This is a natural assumption in practice because log(P̂ (x)) and log(1 − P̂ (x)) are replaced by log(P̂ (x) + ) and log(1− P̂ (x) + ) respectively to avoid explosion when implementing the algorithm. In particular, in our experiments in Section 5, = 0.16 and thus M = 14. Let
F = {f(x, y) = Iy∈[L(x),U(x)] : L,U are output by the NN}, G = {P̂ (x) : P̂ is output by the NN}.
Theorem 4.5. Suppose Assumptions 4.2 and 4.4 hold. The training data D = {(xi, yi), i = 1, 2, · · · , n} where (xi, yi) are i.i.d. samples∼ π. Recall that the (hard) coverage estimator assessment loss is LCA = − 1n ∑n i=1 ( f(xi, yi) log(P̂ (xi)) + (1− f(xi, yi)) log(1− P̂ (xi)) ) . Then for any t > 0, we have
P ( sup
f∈F,P̂∈G |LCA − E[K1(X)]| ≥ t
) ≤ C∗e− nt2 16M2 .
where C∗ only depends on V0 in Assumption 4.2.
Theorem 4.5 shows that the coverage assessment loss approximates E[K1(x)] well with an exponential tail bound. The difficulty in analyzing Theorem 4.5 lies in the fact that the hypothesis classes in Assumption 4.2 (which are constructed by the NN) are different from the hypothesis class used in LCA. To overcome this difficulty, we use the theory of VC-subgraph classes to connect the VC dimension among multiple hypothesis classes, including the class of ψi, the four classes of output functions, and F , log G. Then we establish the covering number bound for the class F and log G, and finally prove Theorem 4.5. To conclude, minimizing E[K1(X)] over P̂ is equivalent to minimizing E[K(X)], which in turn is minimizing a tight upper bound for the coverage assessment loss. Our coverage assessment loss empirically approximates E[K1(X)] well, so that its minimization can ultimately reduce the conditional coverage error.

heading: 5 EXPERIMENTS
text: Experimental Setup. We empirically verify the effectiveness of our proposed CaNet on both synthetic examples and benchmark regression datasets. These datasets have been widely used for the evaluation of methods in regression tasks (Hernández-Lobato & Adams, 2015; Gal & Ghahramani, 2016; Lakshminarayanan et al., 2017; Rosenfeld et al., 2018; Pearce et al., 2018; Zhu et al., 2019). In addition, we adopt the same experimental procedure in Pearce et al. (2018) for data normalization and dataset splitting. To avoid overfitting, we apply a simple network architecture with only 2 hidden layers and each hidden layer has 64 neurons. For each hidden layer, the ReLU activation function is applied to capture the non-linear features. We empirically set the ensemble number m to 5, as the smallest number leading to a stable prediction results. Please refer to Appendix E for implementation details, including those for baseline algorithms.
Evaluation Metrics. To evaluate the conditional coverage estimation of our CaNet, we examine the quality of our Ca-Module measured by the empirical coverage error ECE1 over a partition ∆. ∆ is constructed by equally dividing the width of [0, 1] into M sub-intervals. The value of M depends on the size of the dataset, which is determined by the following strategy: M = min{bthe number of data in validation/50c + 1, 20}. In addition, we also report the IW and CP , where IW = 1n ∑n i=1(U(xi)−L(xi)) (equation 3.1) and CP = 1 n ∑n i=1 ki, (equation 3.2) to show the effectiveness of our predictor module under the high-quality criteria.

heading: 5.1 CONDITIONAL COVERAGE ON SYNTHETIC EXAMPLES
text: In this section, we conduct a series of experiments on synthetic examples to directly compare our prediction results with the ground-truth conditional coverage. In these examples, the conditional coverage can be analytically calculated under the known data distribution. Figure 2 compares the conditional coverage with our predicted coverage under the following settings: x ∼ Uniform[−2, 2] and y|x is drawn from fi(x) = 13 sin(x) + εi(x), x ∈ [−2, 2] where
ε1(x) = 0.1×N(0, 1), ε2(x) = 0.1|x| ×N(0, 1), ε3(x) = 0.1|x| × t4.
N(0, 1) is the standard Gaussian variable and t4 is the standard t random variable with 4 degrees of freedom. Then, the conditional coverage in Definition 2.1 can be analytically calculated as:
P[Y ∈ [L(X), U(X)]|L,U,X = x] = Fi(U(x)− 1
3 sin(x))− Fi(L(x)−
1 3 sin(x))
where Fi is the cumulative distribution function of N(0, 0.12) for i = 1, N(0, (0.1x)2) for i = 2, and 0.1|x| × t4 for i = 3. As shown in Figure 2, the conditional coverages of high-quality PIs at different points diverge among each other and they deviate from the marginal coverage. Thus, having access to the marginal coverage for the whole dataset is not sufficient for decision making, which highlights the need of conditional coverage. In addition, the predicted coverage estimator of our model is highly consistent with the conditional coverage on all of the synthetic examples. These results confirm that our CaNet can accurately estimate the conditional coverage on noisy datasets.

heading: 5.2 PERFORMANCE OF PIS ON BENCHMARK DATASETS
text: In this section, we compare the predictor module of our CaNet on real-world benchmark datasets with following baseline algorithms: (1) nearest-neighbors kernel conditional density estimation (NNKCDE) (Dalmasso et al., 2020), (2) quantile regression forest (QRF) (Meinshausen, 2006), (3) split conformal learning (SCL) (Lei et al., 2018) and (4) the quality-driven PI method (QD-Ens) (Pearce et al., 2018). We quote the results from (Pearce et al., 2018) as a comparison since we share the same experiment setup. Table 1 reports the results of CP , IW and ECE1 for generating PIs at 95% prediction level on benchmark datasets. We employed the criteria in Pearce et al. (2018) to evaluate the performance of PIs: the best IW is achieved by the model with the smallest IW value among those with CP ≥ 95%. As can be seen, our model achieves the best performance on PI generation under the high-quality criteria. With special consideration on interval width quality, it obtains the smallest average interval width (IW ) while maintaining high coverage probability (CP ≥ 95%) on all datasets.

heading: 5.3 PERFORMANCE OF COVERAGE ESTIMATOR ON BENCHMARK DATASETS
text: Coverage for 95% PIs. We use ECE1 to evaluate the coverage estimation performance on realworld datasets as the conditional coverage is unknown. As shown in Table 1, the ECE1 on all experiments are generally around or less than 1%, with better performance on larger datasets. The coverage estimators produced by CaNet have small ECE1 values, which are very close to the perfect-calibrated coverage estimators (Definition 2.2). Compared withECE1 values obtained from the state-of-the-art algorithms in classification tasks (Guo et al., 2017; Kull et al., 2019), the ECE1 values from CaNet are similar and sometimes less than their post-calibrated ECE1 results (usually around 1% to 3%), even though the size of most regression datasets are smaller than the classification datasets. These results further demonstrate that our CaNet can accurately estimate the coverage information of 95% high-quality PIs on real-world regression tasks.
Coverage for PIs at Different Prediction Levels. We conduct multiple experiments on different PI prediction levels to show the robustness of our CaNet. By only modifying the parameter λ1 in Equation 3.5, our Ca-Module could get access to different levels of coverage probability. Table 2 reports the CP, IW and ECE1 values from the CaNet at different PI prediction levels on three benchmark datasets. Results for more datasets can be found in Appendix E. As can be seen, all ECE1 values in Table 2 are fairly small (∼ 1%), demonstrating the stability of our proposed model. Thus, our CaNet can provide accurate coverage estimation on PIs at different prediction levels. These results demonstrate the robustness of our CaNet on real-world datasets, further suggesting its broad applicability.

heading: 6 RELATED WORK
text: Prediction Interval Estimation. High-quality PIs, which can be viewed via a constrained optimization problem where the constraint concerns marginal coverage and the objective is the PI width, has been extensively studied in 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). Such intervals are in the same spirit as the highest density intervals in statistics (Box & Tiao, 2011). While powerful, these approaches could not directly provide the conditional coverage information investigated in this work. Coverageonly criteria, on the other hand, focus solely on coverage satisfaction as the guarantee. These approaches include conformal learning (CL) and its conditional variants (Vovk et al., 2005; 2009; Lei & Wasserman, 2014; Lei et al., 2015; 2018; Kuchibhotla & Ramdas, 2019; Romano et al., 2019; Barber et al., 2019a;b). CL is desirably distribution- or model-free, and in many cases enjoys finitesample guarantees. However, unlike high-quality PIs, they do not explicitly account for the interval width as a quality metric. We also mention conditional density estimation (Holmes et al., 2007; Dutordoir et al., 2018; Izbicki & Lee, 2016; Dalmasso et al., 2020; Freeman et al., 2017; Izbicki et al., 2017) and closely relatedly quantile regression (Koenker & Hallock, 2001; Meinshausen, 2006) as PI construction approaches by converting from the estimated conditional quantile function. These approaches focus on the quality of conditional distribution/quantile, instead of the high-quality criterion only. In this work, we use deep learning to construct high-quality PIs (Khosravi et al., 2010; Pearce et al., 2018; Kivaranovic et al., 2020).
Uncertainty Measurement in Deep Learning. The Bayesian framework offers principled approaches for model uncertainty measurement by computing the posterior distribution over the NN parameters (MacKay, 1992; Neal, 2012). These approaches can also be used to construct PIs. However, they provide a different perspective from the frequentist view taken in this paper, and they focus on the parameter uncertainty of the NN instead of the coverage over a test point. In addition, exact Bayesian inference is computationally intractable for deep neural networks, making it less practical to implement. Gal & Ghahramani (2016) applied a Monte Carlo dropout method to proxy the inference. Directly generated from the networks, softmax response is also commonly used for uncertainty measurement on deep learning models (Bridle, 1990; Lakshminarayanan et al., 2017; Geifman et al., 2018; Sensoy et al., 2018; Ozbulak et al., 2018). Moreover, Niculescu-Mizil & Caruana (2005) showed that NNs typically produce well-calibrated probabilities on binary classification tasks without the need of any post-hoc techniques. In this paper, we also follow this line of work and use the softmax output to access the model uncertainty information.

heading: 7 CONCLUSION
text: In this paper, we identify and investigate the conditional coverage estimation problem for highquality PIs, which is critical for risk-based decision making in regression fields. To address the challenge, we propose an end-to-end algorithm with two powerful modules: coverage assessment module and predictor module. Benefited from these powerful modules, our model can generate and validate the coverage estimator without any requirement for further post-processing steps. In addition, we conduct theoretical analysis to show the effectiveness of our proposed model. Experimental results on synthetic examples and benchmark datasets further demonstrate that our model can robustly provide accurate coverage estimation while simultaneously producing a high-quality PI. Moreover, our Ca-Module can be easily integrated into other deep-learning-based algorithms to get access to their coverage information, opening up more opportunities for broad applications. In the future, we will extend our work by conducting comparison studies with Bayesian methods.

heading: A FURTHER DETAILS ON COVERAGE ESTIMATOR AND COVERAGE ERROR

heading: A.1 COVERAGE PROBABILITY TYPES OF PIS
text: Zhang et al. (2019) introduce the following four coverage probability types of PIs. In general, most of the coverage in PIs considered in the literature falls into one of these types.
Type I: P[Y ∈ [L(X), U(X)]] (marginal coverage); Type II: P[Y ∈ [L(X), U(X)]|L,U ] (conditional coverage given the PI); Type III: P[Y ∈ [L(X), U(X)]|X = x] (conditional coverage given X = x); Type IV: P[Y ∈ [L(X), U(X)]|L,U,X = x] (conditional coverage given the PI and X = x). Note that in high-quality criteria, only Type II coverage is considered in the constraint but Type IV coverage is lacking. Since Type I and III are not considered in our paper, for simplicity Type II coverage is called the marginal coverage, and Type IV coverage is called the conditional coverage in Definition 2.1.
Throughout Appendix, we let A(x) := P[Y ∈ [L(X), U(X)]|L,U,X = x] denote the conditional coverage in Definition 2.1. Moreover, since we only concern Type II and IV coverage in this work, we make the following convention. Throughout the Appendix A-B, P and E should be understood as probability and expectation conditional on L,U . So for A(x), we could simply write A(x) := P[Y ∈ [L(X), U(X)]|X = x] and omit “conditional on L,U”.

heading: A.2 THE TERMINOLOGY “PERFECT-CALIBRATED”
text: The terminology “perfect-calibrated” is borrowed from the confidence calibration in classification tasks. We first review the name of “confidence” in classification.
Confidence calibration is the problem of predicting probability estimates representative of the true correctness likelihood (Guo et al., 2017). Intuitively, an reliable confidence should reflect the true correctness likelihood of the prediction (Kumar et al., 2019). For example, given 100 predictions, each with confidence of 0.8, we expect that 80 should be correctly classified (Guo et al., 2017).
Now let h be the prediction of any models, which is a map from X to Y , trained on the data D. According to their definition, the “best” confidence map P̂ should be the true probability of correctness: P̂ (x) = E[1h(X) is correct for Y |h,X = x], ∀x ∈ X . (A.1) which is a measureable function from X to [0, 1]. In particular, we consider h as a PI in regression tasks and “correctness” of a PI is naturally defined as the success of coverage on the outcome Y . Then the right hand side of Equation (A.1) becomes
E[1h(X) is correct for Y |h,X = x] = E[1Y ∈[L(X),U(X)]|L,U,X = x] = P[Y ∈ [L(X), U(X)]|L,U,X = x].
which is the conditional coverage in our Definition 2.1.
In addition, Guo et al. (2017) introduce the perfect-calibrated confidence as follows:
P(Ŷ = Y |Ŷ , P̂ = p) = p, ∀p ∈ [0, 1] where Ŷ is the class prediction and Ŷ = Y means that the predicted and true class label coincide. Obviously, if P̂ (X) = P(Ŷ (X) = Y |X), i.e., the “best” confidence, then the above equality holds. Transferring this idea into PIs, we can naturally define the perfect-calibrated coverage estimator as
p = E[1h(X) is correct for Y |L,U, P̂ = p] = P[Y ∈ [L(X), U(X)]|L,U, P̂ = p], ∀p ∈ [0, 1]. which is the conditional coverage in our Definition 2.2.

heading: A.3 DETAILS ON COVERAGE ESTIMATOR
text: A perfect-calibrated coverage estimator inherits some properties of a conditional coverage. For example, both of them have the following interpretation: If we have 1000 testing points for a PI, each with the same conditional/perfect-calibrated coverage 0.9, then approximately 900 of them are correctly covered by the PI. Note that the conditional coverage is uniquely defined, but a perfectcalibrated coverage estimator is not necessarily so. Moreover, we have the following facts: (a) The conditional coverage is always perfect-calibrated, but not vice versa. (b) A perfect-calibrated coverage estimator can be viewed as an averaged conditional coverage. (c) A perfect-calibrated coverage estimator is less “informative” than the conditional coverage.
A perfect-calibrated coverage estimator is an averaged conditional coverage. Let P̂ be a (general) coverage estimator. We have
P[Y ∈ [L(X), U(X)]|P̂ (X) = P̂ (x)] =P[Y ∈ [L(X), U(X)]|X ∈ P̂−1(P̂ (x))]
= P[Y ∈ [L(X), U(X)], X ∈ P̂−1(P̂ (x))]
P[X ∈ P̂−1(P̂ (x))]
=
∫ t∈P̂−1(P̂ (x)) E[1Y ∈[L(X),U(X)]|X = t]P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x))]
=
∫ t∈P̂−1(P̂ (x))A(t)P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x))]
Suppose P̂ is a perfect-calibrated coverage estimator. Then we have
P̂ (x) =
∫ t∈P̂−1(P̂ (x))A(t)P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x))]
which implies that P̂ (x) is a weighted average of A(t) over the set P̂−1(P̂ (x)) with weights based on the marginal distribution of X .
A conditional coverage is perfect-calibrated. If P̂ (x) = A(x), then A(t) = A(x) for any t ∈ A−1(A(x)) so∫
t∈A−1(A(x))A(t)P[X ∈ dt] P[X ∈ A−1(A(x))] = A(x)
∫ t∈A−1(A(x)) P[X ∈ dt]
P[X ∈ A−1(A(x))] = A(x).
This shows that A(x) must be a perfect-calibrated coverage estimator. Another way to see this is taking conditional expectation given A(X) = p in the Definition 2.1. Then we get
p = E[A(X)|A(X) = p] = E[P[Y ∈ [L(X), U(X)]|X]|A(X) = p] = P[Y ∈ [L(X), U(X)]|A(X) = p] by the tower property.
A perfect-calibrated coverage estimator may be less informative and may not be the conditional coverage. Suppose we have a PI [L(X), U(X)] at the exact prediction level of 1 − α, i.e., P[Y ∈ [L(X), U(X)]] = 1− α. Then the constant coverage estimator
P̂ (x) = 1− α, ∀x ∈ X
can be viewed as an average coverage estimator over the entire space X . It is a perfect-calibrated coverage estimator since by definition,
P[Y ∈ [L(X), U(X)]|P̂ (X) = 1− α] = P[Y ∈ [L(X), U(X)]] = 1− α = P̂ (X).
But it is not a conditional coverage in general (e.g., the second synthetic example in Section 5).
Now we give an extension of the definition of “perfect-calibrated” coverage estimator, allowing it to be defined on any measurable subsets.
Definition A.1. A coverage estimator P̂ is called a perfect-calibrated coverage estimator on a measurable subset S ⊂ X with P(S) > 0 associated with [L(x), U(x)] if it satisfies
P̂ (x) = P[Y ∈ [L(X), U(X)]|L,U, P̂ (X) = P̂ (x), X ∈ S], a.e. P̂ (x) ∈ [0, 1]. (A.2)
where a.e. is with respect to the probability measure on [0, 1] induced by the random variable P̂ (X)|S . Note that the conditional probability space is standard: (S,FS := {A ∩ S : A ∈ F},PS(A ∩ S) := P(A|S)).
(As our convention in Section A.1, we will omit “conditional on L,U” for simplicity.)
Lemma A.2. (a) A coverage estimator is the conditional coverage if and only if it is a perfectcalibrated coverage estimator on any positive-probability measurable subset S of X . (b) Suppose P̂ is a perfect-calibrated coverage estimator on two disjoint positive-probability measurable subsets S1, S2. Then P̂ is a perfect-calibrated coverage estimator on S1 ∪ S2.
Proof. (a) The proof can be found in Lemma A.4.
(b) We note that by law of total probability,
P[Y ∈ [L(X), U(X)]|P̂ (X) = p,X ∈ S1 ∪ S2] =P[Y ∈ [L(X), U(X)]|P̂ (X) = p,X ∈ S1]P[X ∈ S1|P̂ (X) = p,X ∈ S1 ∪ S2]
+ P[Y ∈ [L(X), U(X)]|P̂ (X) = p,X ∈ S2]P[X ∈ S2|P̂ (X) = p,X ∈ S1 ∪ S2] =p(P[X ∈ S1|P̂ (X) = p,X ∈ S1 ∪ S2] + P[X ∈ S2|P̂ (X) = p,X ∈ S1 ∪ S2]) =p.
Hence P̂ is a perfect-calibrated coverage estimator on S1 ∪ S2.
Lemma A.2(a) is motivated from a theoretical point of view. It provides a guidance that in order to well resemble the conditional coverage, an estimator should be perfect-calibrated on as many subsets on the feature space as possible.

heading: A.4 DETAILS ON COVERAGE ERROR
text: In Section 2, we have introduced CEp to quantify the discrepancy between a coverage estimator and a perfect-calibrated coverage estimator, and C̃Ep to quantify the discrepancy between a coverage estimator and the conditional coverage. We note that by Hölder’s inequality
CEp ≤ CEq, for 1 ≤ p ≤ q ≤ +∞.
A larger value of p corresponds to a larger CE value. Continuing Definition A.1, we can further introduce the calibration-based conditional coverage error on a measurable subset as follows:
Definition A.3. An Lp (1 ≤ p ≤ +∞) calibration-based conditional coverage error, or coverage error for short, of a coverage estimator P̂ on a measurable subset S ⊂ X with P(S) > 0 is defined as:
CEp(S) = ∥∥∥P[Y ∈ [L(X), U(X)]|L,U, P̂ (X), X ∈ S]− P̂ (X)∥∥∥
Lp(S) (A.3)
where Lp-norm is taken with respect to the randomness of P̂ (X) on the conditional probability space (S,FS := {A ∩ S : A ∈ F},PS(A ∩ S) := P(A|S)). In particular, we have CEp := CEp(X ).
Lemma A.4. A coverage estimator P̂ is the conditional coverage if and only if its coverage error CEp(S) = 0 for any measurable subset S ⊂ X with P(S) > 0. In particular, a coverage estimator is the conditional coverage if and only if it is a perfect-calibrated coverage estimator on any measurable subset S of X with P(S) > 0.
Proof. We first show that the conditional coverage is perfect-calibrated. Taking conditional expectation given {A(X) = p,X ∈ S} in the Definition 2.1. Then we get
p = E[A(X)|A(X) = p,X ∈ S] = E[P[Y ∈ [L(X), U(X)]|X]|A(X) = p,X ∈ S] = P[Y ∈ [L(X), U(X)]|A(X) = p,X ∈ S] by the tower property.
So A(x) is a perfect-calibrated coverage estimator on any measurable subset S with P(S) > 0. Hence CEp(S) = 0 for any measurable subsets S ⊂ X with P(S) > 0. On the other hand, similarly to Section A.3, we can express
P[Y ∈ [L(X), U(X)]|P̂ (X) = P̂ (x), X ∈ S] =
∫ t∈P̂−1(P̂ (x))∩S A(t)P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x)) ∩ S] .
Suppose P̂ (x) is not the conditional coverage, then P[P̂ (X) 6= A(X)] > 0. Without loss of generality, we assume P[P̂ (X) > A(X)] > 0. Let S0 := {x ∈ X : P̂ (x) > A(x)}. Note that S0 = ∪+∞n=1{x ∈ X : P̂ (x) > A(x) + 1n}. Since P(S0) > 0, there exists a n0 such that S := {x ∈ X : P̂ (x) > A(x) + 1n0 } and P(S) > 0.
Then for x ∈ S, we have∫ t∈P̂−1(P̂ (x))∩S A(t)P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x)) ∩ S] ≤
∫ t∈P̂−1(P̂ (x))∩S(P̂ (t)− 1 n0 )P[X ∈ dt]
P[X ∈ P̂−1(P̂ (x)) ∩ S] = P̂ (x)− 1 n0 .
Then we have
CEp(S) ≥ CE1(S) ≥ 1
n0 > 0
so CEp(S) > 0, which is a contradiction. Hence P̂ (x) is the conditional coverage.
Next, we describe in detail how to discretize the right-hand-side of equation 2.4 for empirical calculation. We construct a discrete version of (2.2) and then introduce an empirical counterpart of CEp (2.3), which we refer to as Lp empirical calibration-based conditional coverage error ECEp. The ideas behind these are natural extensions of the classification case (Guo et al., 2017; Kull et al., 2019; Kumar et al., 2019; Nixon et al., 2019) into PIs.
We consider the following partition ∆ of [0, 1]. Let [0, 1] be divided into M intervals Im = (am−1, am] (m = 1, · · · ,M ) where 0 = a0 ≤ a1 ≤ · · · ≤ aM = 1. Let Bm = {i = 1, · · · , n : P̂ (xi) ∈ Im}, i.e., the set (bin) of indices i of samples whose prediction coverage estimator P̂ (xi) falls into the interval Im. Note that coverage estimators that are close to each other will fall into the same interval. The coverage probability (i.e., the proportion of successful coverage) in Bm is defined as:
cp(Bm) = 1 |Bm| ∑ i∈Bm 1yi∈[L(xi),U(xi)]. (A.4)
The average estimated coverage in Bm is defined as:
cove(Bm) = 1 |Bm| ∑ i∈Bm P̂ (xi) (A.5)
where P̂ (xi) is the coverage estimator for sample i. cp(Bm) and cove(Bm) approximate the left and right hand sides of (2.2) respectively in the interval Im. A perfect-calibrated coverage estimator should satisfy cp(Bm) = cove(Bm) for all m ∈ {1, · · · ,M}. The diagram of cp(Bm) versus cove(Bm) for all m ∈ {1, · · · ,M} is called the reliability diagram in some literature (Guo et al., 2017).
Based on the partition ∆, we can introduce an empirical version of CEp (which we refer to as ECEp) as:
Definition A.5. The Lp empirical calibration-based conditional coverage error (ECEp) of a coverage estimator P̂ is defined as
ECEp = ∥∥(cp(Bm(i))− cove(Bm(i)))i=1,2,··· ,n∥∥lp . (A.6)
where the lp is the standard p-norm in Rn and Bm(i) is the bin containing sample i.
Equivalently,
ECEp = ( M∑ m=1 |Bm| n |cp(Bm)− cove(Bm)|p ) 1 p , 1 ≤ p < +∞.
ECE∞ = max m=1,2,··· ,M
|cp(Bm)− cove(Bm)| .
Note that ECEp is discontinuous and cannot be used easily for training with gradient-based methods. Therefore, we use the coverage assessment loss LCA introduced in Section 3 for conditional coverage estimation.
Finally, we give some explanations about why the conditional coverage error C̃Ep cannot be used easily for training. Unlike LCA which is unbiased and provides guaranteed estimation accuracy for E[K1(X)] (Theorem 4.5), it is in general not easy to establish an empirical calculation for C̃Ep. Take C̃E2 as an instance. A heuristic argument is to use
θ̂ := 1
n n∑ i=1 (ki − P̂ (xi))2
to approximate
C̃E 2
2 = E[|A(X)− P̂ (X)|2] where ki = 0 or 1 indicates whether each data point has been captured by the PIs; see Section 3. Unfortunately, θ̂ is in general not an unbiased estimator of C̃E 2
2 due to the following observations:
E[θ̂] = E[(k1 − P̂ (x1))2] = E [ E[(k1 − P̂ (x1))2|x1] ] = E [ E[k21|x1]− 2P̂ (x1)E[k1|x1] + P̂ (x1)2
] = E [ E[k1|x1]− 2P̂ (x1)E[k1|x1] + P̂ (x1)2
] = E [ A(x1)− 2P̂ (x1)A(x1) + P̂ (x1)2
] = E [ A(x1) 2 − 2P̂ (x1)A(x1) + P̂ (x1)2 ] + E[A(x1)−A(x1)2]
= E[|A(X)− P̂ (X)|2] + E[A(X)−A(X)2].
Hence E[θ̂] 6= C̃E 2
2 unless A(X) = 0 or 1 almost surely. Note that the gap E[A(X)−A(X)2] does not depend on n and thus not vanish as n grows. To obtain a reasonable estimator for C̃Ep, one need more information about A(x) besides the marginal coverage E[A(X)]. However, this information is usually hard to obtain locally because of the nature of (possibly high-dimensional) feature space. In fact, this has led us to the idea of “dimension reduction” of the sample space to the space [0, 1] and the definition of CEp.

heading: B MATHEMATICAL DEVELOPMENTS FOR THEOREM 4.1
text: This section proves Theorem 4.1: we show that both coverage error and conditional coverage error are tightly bounded above by the expectation of a Kullback–Leibler divergence-type random variable K1(x). This means that minimizing E[K1(X)] can recover the true conditional coverage and effectively reduce the coverage error. The proof consists of several inequalities regarding coverage errors and their relations. We first begin with the following connection between CEp and C̃Ep.
Theorem B.1. For any PI and its associated P̂ , the Lp coverage error is always less than or equal to the Lp conditional coverage error, i.e.,
CEp ≤ C̃Ep, ∀1 ≤ p ≤ +∞
Proof. We note that the function t 7→ |t|p is a convex function. We also note that σ(P̂ (X)) ⊂ σ(X) where σ(Y ) represents the σ-field generated by a random variable Y .
C̃E p p =E [∣∣∣A(X)− P̂ (X)∣∣∣p]
=E [ E|A(x)− P̂ (X)|p|P̂ (X)] ] ≥E
[∣∣∣E[A(x)− P̂ (X)|P̂ (X)]∣∣∣p] by Jensen’s inequality =E
[∣∣∣E[1Y ∈[L(X),U(X)]|P̂ (X)]− P̂ (X)∣∣∣p] by the tower property =CEpp
Therefore we have CEp ≤ C̃Ep, ∀1 ≤ p ≤ +∞.
Next we have the following bounds on Lp conditional coverage error:
Theorem B.2. The Lp conditional coverage error is bounded above by a power function of the L2 conditional coverage error. Formally,
C̃Ep ≤ C̃E αp 2 , ∀1 ≤ p ≤ +∞,
where αp = 1,∀1 ≤ p ≤ 2 and αp = 2p ,∀2 ≤ p ≤ +∞.
Proof. By Hölder’s inequality,
C̃Ep ≤ C̃E2, if 1 ≤ p ≤ 2.
Since 0 ≤ |A(x)− P̂ (x)| ≤ 1,
then, |A(x)− P̂ (x)|p ≤ |A(x)− P̂ (x)|2, ∀p ≥ 2
and thus
C̃E p p = E[|A(X)− P̂ (X)|p] ≤ E[|A(X)− P̂ (X)|2] ≤ C̃E 2 2, ∀p ≥ 2.
Next, recall that
K(x) = A(x) log
( A(x)
P̂ (x)
) + (1−A(x)) log ( 1−A(x) 1− P̂ (x) ) ,
K0(x) = A(x) log(A(x)) + (1−A(x)) log(1−A(x)),
K1(x) = −A(x) log(P̂ (x))− (1−A(x)) log(1− P̂ (x)). Theorem B.3. The L2 conditional coverage error is bounded above by the expectation of K(x). Formally,
C̃E αp 2 ≤ ( 1
2 E[K(X)]
)αp/2 ,
where αp is defined in Theorem B.2.
Proof. For any fixed x, consider two random variables with Bernoulli distributions:
W1 = { 1 w.p. A(x), 0 w.p. 1−A(x).
W2 = { 1 w.p. P̂ (x), 0 w.p. 1− P̂ (x).
Let Pi be the distribution ofWi. It follows from Pinsker’s inequality, e.g., Theorem 2.16 in (Massart, 2007), that
‖P1 − P2‖2TV ≤ 1
2 K(P1, P2).
where TV denotes the total variation distance and K denotes the KL divergence. Since Pi is the Bernoulli distribution, we can express it as
|A(x)− P (x)|2 ≤ 1 2 K(x)
Taking expectation, we obtain
E[|A(X)− P (X)|2] ≤ 1 2 E[K(X)].
Hence,
C̃E αp 2 ≤ ( 1
2 E[K(X)]
)αp/2 .
Combining Theorem B.1, B.2 and B.3, we immediately conclude that: Theorem B.4 (Restated Theorem 4.1).
CEp ≤ C̃Ep ≤ C̃E αp 2 ≤ ( 1
2 E[K(X)]
)αp/2 , ∀1 ≤ p ≤ +∞
where αp is defined in Theorem B.2. Moreover, all inequalities are attainable, e.g., if P̂ (x) is the conditional coverage A(x).

heading: C JUSTIFICATION OF ASSUMPTION 4.2 AND MATHEMATICAL DEVELOPMENTS FOR THEOREM 4.5
text: In this section, we analyze the rationality of Assumption 4.2 and build essential ingredients for proving Theorem 4.3 and 4.5. The difficulty in analyzing Theorem 4.5 lies in the fact that the hypothesis classes in Assumption 4.2 (which are constructed by the NN) are different from the hypothesis class used in LCA. To overcome this difficulty, we use the theory of VC-subgraph classes to analyze the connection between the VC dimension of the two hypothesis classes.

heading: C.1 REVIEW OF THE VC DIMENSION
text: For self-contained purpose, we first review the definition of the VC-subgraph class and VC dimension. Definition C.1. Consider an arbitrary collection {x1, · · · , xn} of points in a set X and a collection C of subsets of X . We say that C shatters {x1, · · · , xn} if all of 2n possible subsets of {x1, · · · , xn} can be written as A = C ∩ {x1, · · · , xn} for some C ∈ C. The VC dimension V (C) of the class C is the smallest n for which no set of size n {x1, · · · , xn} is shattered by C. If C shatters sets of arbitrarily large size, we set V (C) =∞. We say that C is a VC-class if V (C) <∞.
In some literature, the VC dimension V (C) of the class C is alternatively defined as the largest n for which there exists a set of size n {x1, · · · , xn} shattered by C, i.e., it is the value in definition C.1 minus 1. We can more formally define the VC dimension by the growth function as follows:
Definition C.2. Define the nth shatter coefficient (or growth function) of C as
ΠC(n) := max x1,··· ,xn
|{A : A = C ∩ {x1, · · · , xn} for some C ∈ C}|
Then V (C) := inf{n : ΠC(n) < 2n}.
Definition C.3. For a function f : X → R, the subset of X × R given by {(x, t) : t < f(x)} is the (open) subgraph of f . A collection F of measurable real functions on the sample space X is a VC-subgraph class or VC-class, if the collection of all subgraphs of functions in F forms a VC-class of sets (as sets in X × R). Let V (F) denote the VC dimension of the set of subgraphs of F . Lemma C.4. In Definition C.3, the open subgraph of f , {(x, t) : t < f(x)}, can be replaced by the close subgraph {(x, t) : t ≤ f(x)}, the close supergraph {(x, t) : t ≥ f(x)} or the open supergraph {(x, t) : t > f(x)}. All of them lead to the equivalent definition of the VC-class and the equal VC dimension.
Proof. This result follows from Lemma 9.33 and Lemma 9.9(iv) in Kosorok (2007).
For indicator functions of sets, we have the following equivalence. Lemma C.5. For any class C of sets in a set X , the class FC of indicator functions of sets in C is a VC-class if and only if C is a VC-class. Moreover, whenever at least one of C or FC is VC-class, the respective VC dimensions are equal.
Proof. This is Lemma 9.8 in Kosorok (2007). Note that the sets of C are in X while the subgraphs of functions of FC are in X × R.

heading: C.2 JUSTIFYING ASSUMPTION 4.2
text: We first restate the assumption:
Assumption C.6 (Restated Assumption 4.2). The four classes of functions ([L(x), U(x)], P̂ (x), 1− P̂ (x)) output by the neural network (NN) in Figure 1 have finite VC dimensions, say they are bounded above by V0.
In Figure 1, the output four neurons of the NN are denoted as (L(x), U(x), P̂ (x), 1 − P̂ (x)). We further let (ψ1(x), ψ2(x), ψ3(x), ψ4(x)) denote the pre-activated values of (L(x), U(x), P̂ (x), 1 − P̂ (x)). In other words,
L(x) = min(ψ1(x), ψ2(x)),
U(x) = max(ψ1(x), ψ2(x)),
P̂ (x) = softmax(ψ3(x)||ψ4(x)) = σ(ψ3(x)− ψ4(x)),
1− P̂ (x) = softmax(ψ4(x)||ψ3(x)) = σ(ψ4(x)− ψ3(x)), where σ is the sigmoid function. Let the function classes
H1 = {L(x) : L is output by the the NN},
H2 = {U(x) : U is output by the the NN},
G = {P̂ (x) : P̂ is output by the NN},
1− G = {1− P̂ (x) : P̂ is output by the NN}.
Assumption 4.2 holds for a wide range of NNs, in particular the one we adopt in the experiments (where we use the ReLU-activated NN to construct ψi, i = 1, 2, 3, 4 ; see Section 5). Our first result is to concretely show that the four NN outputs above,H1,H2, G and 1−G, under the ReLU setting, all have finite VC dimensions and thus satisfy Assumption 4.2. Theorem C.7 (Restated Theorem 4.3). Suppose ψi, i = 1, 2, 3, 4 are the pre-activated output neurons of the NN in Figure 1 using the ReLU activation function. Then Assumption 4.2 holds. Moreover, suppose the NN has W parameters and U computation units (nodes). Then V0 = O(WU).
Proof. First, we look at L(x) and U(x). Note that the class of ψi (i = 1, 2) is constructed by a NN with the ReLU activation function. Therefore by Theorem 8 in Bartlett et al. (2019) (see also C.8 below), V ({ψ1}) = O(WU) <∞, V ({ψ2}) = O(WU) <∞. By Lemma 9.9 (i) in Kosorok (2007), we have
V (H1) = V ({min(ψ1, ψ2)}) ≤ V ({ψ1}) + V ({ψ2})− 1 = O(WU) <∞.
By Lemma 9.9 (ii) in Kosorok (2007), we have
V (H2) = V ({max(ψ1, ψ2)}) ≤ V ({ψ1}) + V ({ψ2})− 1 = O(WU) <∞.
Next, we look at P̂ (x) and 1 − P̂ (x). We add an additional neuron after the layer where ψ3, ψ4 stand. This neuron is defined as ψ5 = ψ3 − ψ4 which is a linear combination of ψ3 and ψ4. Note that the class of ψ5 is constructed by a NN with the ReLU activation function and linear activation function (by adding one unit and two parameters in the originial NN). Therefore by Theorem 8 in Bartlett et al. (2019), V ({ψ5}) = O(WU) <∞ By Lemma 9.9 (viii) in Kosorok (2007), we have
V (G) = V ({σ(ψ5)}) ≤ V ({ψ5}) = O(WU) <∞.
since σ is a monotone function. Again, by Lemma 9.9 (viii) in Kosorok (2007), we have
V (1− G) ≤ V (G) = O(WU) <∞.
since t 7→ 1− t is a monotone function.
We also list some results for other activations here. From these results, and using the same argument as above, we see that Assumption 4.2 holds similarly for all these activations. Lemma C.8. Suppose the class of functions is constructed by a NN with W parameters and U units with activation functions that are piecewise polynomials with at most p pieces and of degree at most d. Then it has VC dimension O(WU log((d+ 1)p)).
Proof. This is Theorem 8 in Bartlett et al. (2019).
Note that the activation functions in Lemma C.8 include in particular the ReLU activation and linear activation. Lemma C.9. Suppose the class of functions is constructed by a NN with W parameters with binary as well as linear activation function. Then it has VC dimension O(W 2).
Proof. This is Theorem 5 in Sontag (1998).
Lemma C.10. Suppose the class of functions is constructed by a NN with W parameters and U units with activation function that is the standard sigmoid function (except that the output unit being a linear threshold unit). Then it has VC dimension O(W 2U2).
Proof. This is Theorem 8.13 in Anthony & Bartlett (2009).

heading: C.3 CONNECTIONS AMONG DIFFERENT HYPOTHESIS CLASSES
text: To prove Theorem 4.5, we need to study several building blocks on the relations between different hypothesis classes. Our first observation is: Theorem C.11. Suppose V (G) < +∞. Then all of the following classes have VC dimension ≤ V (G):
1− G := {1− P̂ (x) : P̂ is output by the NN}.
G′ := log(G) := {log(P̂ (x)) : P̂ is output by the NN}.
log(1− G) := {log(1− P̂ (x)) : P̂ is output by the NN}.
Proof. The result follows from Lemma 9.9 (viii) in Kosorok (2007) since all of the transformations are monotone functions.
Our second observation is about
F = {f(x, y) = Iy∈[L(x),U(x)] : L,U are output by the NN}.
Note that the domain of functions in F is different from the domain of functions inHi (i = 1, 2) as it includes the outcome space. Below we derive a result that connects the VC dimension ofHi with that of F . Theorem C.12. Suppose V (Hi) ≤ V0 (i = 1, 2). We have that
V (1−F) ≤ V (F) ≤ 10(V0 − 1) < +∞
where 1−F := {1− f(x, y) : f ∈ F}.
Proof. The first inequality follows from Lemma 9.9 (viii) in Kosorok (2007). We consider the following two classes:
F1 := {IL(x)≤t : L ∈ H1, t ∈ R},
F2 := {IU(x)≥t : U ∈ H2, t ∈ R}. Since the functions in F1 are all indicator functions, by Lemma C.5,
V (F1) = V ({{(x, t) : L(x) ≤ t} : L ∈ H1, t ∈ R}).
Note that the latter is the VC dimension of the close supergraphs of all functions in H1. Then by Definition C.3 and Lemma C.4, we have
V ({{(x, t) : L(x) ≤ t} : L ∈ H1, t ∈ R}) = V (H1)
Therefore we have V (F1) = V (H1) ≤ V0
Similarly, V (F2) = V (H2) ≤ V0
Note that we can write Iy∈[L(x),U(x)] = IL(x)≤yIU(x)≥y
By the definition of growth functions,
ΠF (m) := max (x1,y1),··· ,(xm,ym) ∣∣{(Iy1∈[L(x1),U(x1)], · · · , Iym∈[L(xm),U(xm)]) : L ∈ H1, U ∈ H2}∣∣ ≤ max
(x1,y1),··· ,(xm,ym) ∣∣{(IL(x1)≤y1 , · · · , IL(xm)≤ym) : L ∈ H1}∣∣× max (x1,y1),··· ,(xm,ym)
∣∣{(IU(x1)≥y1 , · · · , IU(xm)≥ym) : U ∈ H2}∣∣ =ΠF1(m)ΠF2(m)
≤ ( em
V0 − 1 )2(V0−1) for allm ≥ V0 where the last inequality is due to the Sauer–Shelah lemma. Takingm = 10(V0−1), we obtain (
em
V0 − 1
)2(V0−1) = (10e) 2(V0−1) ≤ 750V0−1 < 2m
Combining the above inequality, we have
ΠF (m) < 2 m.
This shows that V (F) ≤ m = 10(V0 − 1).

heading: C.4 PROOF OF THEOREM 4.5
text: This subsection proves Theorem 4.5. Recall that
F = {Iy∈[L(x),U(x)] : L,U are output by the NN},
G = {P̂ (x) : P̂ is output by the NN}. G′ := log(G) := {log(P̂ (x)) : P̂ is output by the NN}.
Let N( ,F , L2(Q)) denote the covering number, i.e., the minimal number of balls {g : ‖g − h‖L2(Q) < } of radius needed to cover the set F . We need the following bounds: Lemma C.13. Suppose F is a class of functions f : X × Y → [0, 1] with a finite VC dimension V (F). For every 0 < < 1,
sup Q
logN( ,F , L2(Q)) ≤ K2 ( 1 ) 1 e
where the constant K2 depends on V (F) only.
Proof. It follows from Theorem 2.6.7 in Van der Vaart & Wellner (1996) that there exists a universal constant K such that
sup Q N( ,F , L2(Q)) ≤ KV (F)(16e)V (F)
( 1 )V (F)−1 for any 0 < < 1. Since (
1
) 1 e ≥ max ( log ( 1 ) , 1 ) , ∀0 < < 1,
we have
sup Q
logN( ,F , L2(Q)) ≤ K3 + (V (F)− 1) log ( 1 ) ≤ K2 ( 1 ) 1 e
where K3 := log(KV (F)(16e)V (F)) and K2 = K3 + V (F)− 1 only depending on V (F).
We remark that a similar result can also be obtained for the class 1−F by Theorem C.12. Lemma C.14. Suppose G is a class of functions P̂ : X → [0, 1] with a finite VC dimension V (G) and | log(P̂ (x))| ≤M . Let G′ := {log(P̂ ) : P̂ ∈ G}. Then, for every 0 < < 1,
sup Q
logN( M,G′, L2(Q)) ≤ K2 ( 1 ) 1 e
where the constant K2 depends on V (G) only.
Proof. First note that φ(t) := log(t) is a monotone function. Hence G′ := {log(P̂ ) : P̂ ∈ G} is a VC-class with VC dimension ≤ V (G) by Lemma 9.9 (viii) in Kosorok (2007). The rest of the proof is similar to C.13.
We remark that a similar result can also be obtained for the class log(1− G) by Theorem C.11. Next we restate Assumption 4.4:
Assumption C.15 (Restated Assumption 4.4). | log(P̂ (x))| ≤ M , | log(1 − P̂ (x))| ≤ M for all x and P̂ .
As discussed in Section 4, this is a natural assumption in practice because log(P̂ (x)) and log(1 − P̂ (x)) are replaced by log(P̂ (x) + ) and log(1 − P̂ (x) + ) respectively to avoid explosion when implementing the algorithm. In particular, in our experiments in Section 5, = 0.16 and thus M = 14.
We are now ready to prove Theorem 4.5:
Theorem C.16 (Restated Theorem 4.5). Suppose Assumptions 4.2 and 4.4 hold. The training data D = {(xi, yi), i = 1, 2, · · · , n} where (xi, yi) are i.i.d. samples ∼ π. Recall that the (hard) coverage estimator assessment loss is
LCA = − 1
n n∑ i=1 ( f(xi, yi) log(P̂ (xi)) + (1− f(xi, yi)) log(1− P̂ (xi)) ) .
Then for any t > 0, we have
P ( sup
f∈F,P̂∈G |LCA − E[K1(X)]| ≥ t
) ≤ C∗e− nt2 16M2
where C∗ only depends on V0 in Assumption 4.2.
Proof. Note that E[f(xi, yi)|xi] = A(xi) for any fixed L and U . Taking expectation on LCA, we have
E[LCA] = E [E[LCA|x1, x2, · · · , xn]]
= E [ − 1 n n∑ i=1 ( A(xi) log(P̂ (xi)) + (1−A(xi)) log(1− P̂ (xi)) )] = E[K1(X)].
We consider the first part A(x) log(P̂ (x)). The second part can be done using the same argument. Note that by Theorem 9.15 in Kosorok (2007), we have
sup Q logN( M,F · G′, L2(Q)) ≤ sup Q logN( /2,F , L2(Q)) + sup Q logN( M/2,G′, L2(Q)).
Consider the class 12 + 1 2MF · G ′ := { 12 + 1 2M (f(x, y) log(P̂ (x))) : f ∈ F , log(P̂ ) ∈ G ′} which consists of functions taking values in [0, 1]. We have
sup Q
logN( , 1
2 +
1
2M F · G′, L2(Q))
= sup Q
logN(2 M,F · G′, L2(Q))
≤ sup Q logN( ,F , L2(Q)) + sup Q logN( M,G′, L2(Q))
≤K2 ( 1 )1/e where the last inequality follows from Lemma C.13 and Lemma C.14, and K2 only depends on V (F) and V (G). (Recall that we have shown V (G′) ≤ V (G) in Lemma C.14.) Moreover, by Theorems C.11 and C.12, we can claim that K2 only depends on V0. This inequality shows that 1 2 + 1 2MF · G
′ satisfies the conditions in Theorem 2.14.10 in Van der Vaart & Wellner (1996) and thus for every δ > 0 and t > 0,
P ( sup
φ∈ 12+ 1 2M F·G′ ∣∣∣∣∣ 1n n∑ i=1 φ(xi, yi)− E[φ(x, y)] ∣∣∣∣∣ ≥ t ) ≤ CeD( √ nt)U+δe−2nt 2
where U = 1 e (6− 1 e )
2+ 1e < 1 and the constants C and D depend on K2 and δ only. Let δ = 1−U . Note
that −2( √ nt)2 +D( √ nt) ≤ −( √ nt)2 + (D/2)2.
Hence we have
P ( sup
φ∈ 12+ 1 2M F·G′ ∣∣∣∣∣ 1n n∑ i=1 φ(xi, yi)− E[φ(x, y)] ∣∣∣∣∣ ≥ t ) ≤ C∗e−nt 2
where C∗ only depends on K2, or, only depends on V0.
This shows that
P ( sup
f∈F,P̂∈G ∣∣∣∣∣ 1n n∑ i=1 f(xi, yi) log(P̂ (xi))− E[A(x) log(P̂ (x))] ∣∣∣∣∣ ≥ t ) ≤ C∗e− nt2 4M2 .
A similar result can be established for the second part since the hypothesis classes there have been studied in Theorem C.11 and C.12:
P ( sup
f∈F,P̂∈G ∣∣∣∣∣ 1n n∑ i=1 (1− f(xi, yi)) log(1− P̂ (xi))− E[(1−A(x)) log(1− P̂ (x))] ∣∣∣∣∣ ≥ t )
≤C∗e− nt2 4M2 .
Combining the two parts and noting the following fact:
{sup |γ + β| ≥ t} ⊂{sup |γ|+ sup |β| ≥ t}
⊂{sup |γ| ≥ t 2 } ∪ {sup |β| ≥ t 2 },
we conclude that
P ( sup
f∈F,P̂∈G |LCA − E[K1(x)]| ≥ t
) ≤ C∗e− nt2 16M2
where C∗ only depends on V0.
Lastly, the following corollary explicitly connects our theoretical developments to the experimental setup:
Corollary C.17. Suppose the NN is designed as the one specified in the experiments (Section 5). The training data D = {(xi, yi), i = 1, 2, · · · , n} where (xi, yi) are i.i.d. samples ∼ π. Then for any t > 0, we have
P ( sup
f∈F,P̂∈G |LCA − E[K1(X)]| ≥ t
) ≤ C∗e− nt2 16M2 .
where C∗ only depends on V0 in Assumption 4.2.
Proof. We note that Assumptions 4.2 and 4.4 hold in this case by Theorem 4.3 and the observation after Assumption 4.4. So Theorem 4.5 implies Corollary C.17.

heading: D ALGORITHM DETAILS
text: We provide additional algorithm details for our framework in Section 3. Algorithm 1 is the description of our tuning procedure for hyper-parameters λ1, λ2, λ3. Let the marginal coverage probability CPD′ and the average coverage estimation ACD′ on the validation set D′ be defined as
CPD′ = 1 |D′| ∑ i∈D′ 1yi∈[L(xi),U(xi)], ACD′ = 1 |D′| ∑ i∈D′ P̂ (xi).
where ([L(x), U(x)], P̂ (x)) are prediction results from the deep ensemble. Then, λi, i = 1, 2, 3 are adjusted to ensure that CPD′ coincides roughly with ACD′ , and CPD′ attains the target prediction level.

heading: E EXPERIMENTAL DETAILS AND MORE RESULTS
text: This section illustrates experimental details and more experimental results from our proposed model.
Algorithm 1: Tuning algorithm Goal: Tune hyperparameters λ1, λ2, and λ3; Input: Prediction level 1− α, training dataset D, validation dataset D′; Procedure: (1) Initialize λi (i = 1, 2, 3) so that CPD′ is nontrivial, i.e., not (almost) 0 or 1. (2) While CPD′ is nontrivial: tune λ2 and λ3 so that |CPD′ −ACD′ | ≤ (e.g., = 1%.) (3) Otherwise tune λ1 such that CPD′ is nontrivial. Do step 2 again until we find λ2 and λ3. (4) Tune λ1 such that CPD′ > 1− α where λ2 and λ3 are fixed from (3). Output: λ1, λ2, and λ3.

heading: E.1 EXPERIMENTAL DETAILS
text: Table 3 gives a detailed description about the datasets we use. These open-access real-world benchmark regression datasets are widely used for the evaluation of methods in regression tasks (Hernández-Lobato & Adams, 2015; Gal & Ghahramani, 2016; Lakshminarayanan et al., 2017; Rosenfeld et al., 2018; Pearce et al., 2018; Zhu et al., 2019).
For synthetic datasets, 2000 i.i.d data are generated for each synthetic setting and randomly split into 1000 training data and 1000 testing data. For benchmark datasets, we first do the data normalization and then randomly split 80% data for training and 20% for testing. The choice of 80%/20% split, compared to the 90%/10% split in Pearce et al. (2018), is motivated from the need to increase the test size in order to get a meaningfulECE evaluation. The latter is due to that evaluating ECE requires binning, where using a larger number of bins approximates more closely CE, but also requires a larger test size to sustain enough statistical quality for the resulting ECE estimate. This delicate tradeoff motivates us to increase the share of the test set in our split. Following Pearce et al. (2018), our hyper-parameters are selected using the validation set from a random split. Then, they are fixed during the evaluation on other random splits.
As specified in Section 2 (Equation (2.4)) and Appendix A.4 (Equation (A.6)), ECE is evaluated based on dividing [0,1] into M sub-intervals. The larger M is, the more precise is in using ECE to approximate CE, the latter being the ideal conditional coverage error estimator. On the other hand, a larger test set size is needed to support the use of a larger M without deteriorating the statistical quality of the ECE. This delicate tradeoff motivates us to increase the share of the test set in our split.
For synthetic datasets, the hyperparameters and corresponding results in Figure 2 are: (a) λ1 = 1.7, λ2 = 10−5, λ3 = 1500, CP = 0.95, IW = 0.40, ECE1 = 0.62%. (b) λ1 = 1.9, λ2 = 10−5, λ3 = 1000, CP = 0.96, IW = 0.40, ECE1 = 0.12%. (c) λ1 = 3.4, λ2 = 10−5, λ3 = 1000, CP = 0.95, IW = 0.50, ECE1 = 0.65%.
For benchmark datasets, the implementation details for baseline algorithms in Table 1 are:
(1) Nearest-neighbors kernel conditional density estimation (NNKCDE). The algorithm is based on Section 2.1 in Dalmasso et al. (2020). We use the same Python code provided by Dalmasso et al. (2020) with the default Gaussian kernel. Two tuning parameters, i.e., the number of nearest neighbors k and the bandwidth h of the smoothing kernel, are chosen in a principled way by minimizing the CDE loss on validation data, the same way as in Dalmasso et al. (2020).
(2) Quantile regression forest (QRF). The algorithm is based on Meinshausen (2006). We use the RandomForestQuantileRegressor from the package scikit-garden in Python.
(3) Split conformal learning (SCL). The algorithm based on Algorithm 2 in Lei et al. (2018). The regression algorithm inside SCL that we use is a neural network with mean square loss. The neural network has the same structure of hidden layers as in Section 5.

heading: E.2 ADDITIONAL EXPERIMENTAL RESULTS
text: Table 4 gives additional experimental results.

heading: E.3 COMPARISONS WITH A TWO-STAGE APPROACH
text: We compare the performance of CaNet with a two-stage approach to further demonstrate the effectiveness of our Ca-Module. The two-stage approach is implemented with two separate steps: (1) given a regression dataset, we train a neural network to generate the prediction interval, (2) after getting the predictor, we train another network to estimate the conditional coverage of the PI from the previous stage using LCA loss.
Figure 3 compares the reliability diagrams (introduced in A.4) and the coverage histograms of CaNet and the two-stage approach on the dataset “Protein”. The coverage histograms demonstrate the percentage of samples in each bin Bm (equation A.5) for m ∈ {1, · · · ,M}. The average estimated coverage of our model closely matches its coverage probability, while the average estimated coverage of the two-stage algorithm is substantially lower than its coverage probability. In addition, the ECE1 from CaNet (0.77%) is much lower than the ECE1 from the two-stage approach (3.9%).

Review summaries:
In the submitted paper, the authors study high-quality prediction intervals (PIs). The paper proposes a novel design of loss functions to generate PIs and conditional coverage estimates. The theoretical justification for using the conditional coverage error (in Ca-module) is presented and the numerical experiments with promising results are provided on multiple benchmark datasets.

In the paper, the author addresses a calibration-based conditional coverage error in order to avoid the difficulty of conditional coverage, which provides a middle ground between marginal (no conditional information) and conditional coverage (high computational cost). The author generates the idea building on prior work and designs a new loss function combining the high-quality criterion and a coverage assessment loss. The theoretical framework about the loss function is laid out clearly, and the performances on benchmark datasets provide accurate results which outperform the other baseline algorithms on high-quality prediction intervals generation.