On Optimal Data Splitting for Split Conformal Prediction
arXiv:2606.31600v1 [math.ST] 30 Jun 2026
Sayan Das∗
Bahram Yaghooti†
Todd A. Kuffner∗
Soumendra N. Lahiri∗
Abstract Conformal prediction and its variants, including the split conformal prediction, provide a distribution-free framework for uncertainty quantification by constructing prediction intervals or sets with finite-sample coverage guarantees. The statistical efficiency of these intervals depends critically on how the data are split into training and calibration samples. Despite its practical importance, a principled characterization of the training-calibration split that minimizes prediction interval length while maintaining coverage has remained largely unresolved. In this paper, we develop a theoretical framework for optimal data splitting in split conformal prediction. We first analyze the problem in a general setting and derive analytical characterizations of the length-optimal split ratio under both symmetric and asymmetric regimes. We then show how the general results specialize to several commonly used regression settings, including linear regression, nonparametric regression, and neural networks, thereby demonstrating the scope of the framework. We also describe a data-based method for selecting the optimal proportion. Our analysis clarifies how model-related features govern the optimal allocation of samples between training and calibration and provides principled guidance for constructing shorter prediction intervals. Experiments on both synthetic and real-world datasets demonstrate the applicability of the proposed methodology across a variety of practical scenarios.
Keywords: Conformal prediction; optimal sample allocation; asymptotic optimality; distributionfree inference. MSC2020: Primary 62G15, 62G20; Secondary 62G08.
1
Introduction
Conformal prediction is a distribution-free framework for uncertainty quantification that constructs prediction intervals or sets with finite-sample marginal coverage guarantees, assuming only data exchangeability [42]. Its key advantage lies in providing rigorous, model-agnostic calibration through a post-hoc mechanism that converts point predictions into statistically valid confidence intervals or sets, without modifying the underlying learning algorithm [2]. Classical conformal prediction, also known as full conformal prediction, guarantees exact coverage but requires retraining the model for each possible output, limiting its practical applicability in large-scale settings due to substantial computational complexity. To address these computational challenges, several variants have emerged. Among these, inductive conformal prediction, also known as split conformal prediction, stands out for its practicality: ∗
Department of Statistics and Data Science, Washington University in St. Louis, St. Louis, MO 63130, USA Department of Electrical and Systems Engineering, Washington University in St. Louis, St. Louis, MO 63130, USA †
1
it partitions data into (proper) training- and calibration-subsets and achieves valid coverage with a single model fit [29]. Other extensions, such as cross-conformal prediction [41] and Jackknife+ [6], enhance the stability and sharpness of prediction intervals through resampling techniques. Additionally, online conformal prediction extends the framework to sequential or streaming settings, enabling real-time uncertainty quantification [4, 17]. Beyond computational efficiency, recent research has substantially expanded the applicability of conformal prediction. Adaptive techniques [46], blockwise calibration methods [45], and covariate-adjusted conformal inference [39] have been developed to address challenges such as nonexchangeable data and distributional shift. Further efforts focus on achieving stronger forms of validity, including class-conditional coverage guarantees [13]. In addition, [12] examines the relationship between calibration set size and predictive accuracy under computational constraints, underscoring the importance of principled data allocation strategies. Despite split conformal prediction’s theoretical guarantees and widespread adoption across a broad range of application domains, including medical diagnostics [36], natural language processing [8], computer vision [3, 28], finance [30], medical imaging [26], robotics and autonomous vehicles [27], reinforcement learning [16], and large language models [9], determining optimal split ratio and its impact on prediction interval length remains a major challenge. Prior work has analyzed the expected size of conformal prediction sets in both asymptotic and finite-sample regimes [12, 22, 24, 33, 43]. However, none of these studies provide procedures for deriving or achieving optimal-length prediction intervals, leaving open the question of how to systematically determine the split ratio that yields the smallest valid prediction interval. Another growing line of work has begun to address this gap by explicitly developing methods for prediction interval length optimization. Several approaches focus on designing improved conformity scores that yield sharper and more informative sets under various forms of conditional validity [11, 15, 25, 31, 32, 44]. Another direction leverages insights from level-set estimation, drawing on foundations from classical statistics and their adaptation to conformal prediction in [23, 33], to construct covariate-dependent thresholds that directly target shorter prediction sets while accommodating refined coverage notions. Complementing these developments, conformal training methods [5, 9, 38] optimize the conformity score itself to reduce prediction-set size, whereas alternative formulations keep the score fixed and instead optimize adaptive, covariate-informed thresholds. More recently, a principled minimax framework for constructing optimal length prediction sets under various conditional coverage requirements has been introduced, further emphasizing the growing focus on prediction-set efficiency [20]. Collectively, these advances mark an important shift toward methods that not only ensure proper coverage but also aim to make conformal prediction sets as compact and informative as possible. Nevertheless, split conformal prediction continues to face a critical limitation: its performance is highly sensitive to how the data are partitioned between training and calibration sets. In practice, most implementations rely on fixed or randomly chosen split ratios without principled justification, which can lead to suboptimal performance, especially in data-limited or imbalanced settings. This challenge introduces a fundamental trade-off: Allocating more data to training improves predictive accuracy, whereas increasing the size of the calibration set enhances the reliability of the uncertainty quantification. Yet, despite this inherent tension, no principled framework currently exists for selecting the optimal split ratio, leaving a significant gap in the practical deployment of split conformal prediction. In this paper, we address the length-optimal data splitting problem in split conformal prediction by introducing a theoretical framework for selecting the optimal training-calibration sample size ratio. Our contributions are as follows:
2
• We derive analytical expressions that characterize the optimal training-calibration split ratio required to jointly optimize the prediction interval length while preserving nominal coverage accuracy. These results are obtained under general distributional assumptions, including both symmetric and asymmetric settings, and offer insights into the dependence of the optimal split on key distributional parameters. • We instantiate these theoretical findings within three representative classes of regression models, namely linear regression, nonparametric regression, and neural network-based regression, thereby illustrating the broad applicability and model-agnostic nature of the proposed framework. We further analyze the theoretical and empirical implications of our results in each case. • We also provide a general data based method for determining the optimal training-calibration split ratio that is agnostic to the learning algorithm. Applied to the three specific classes of problems mentioned above, the proposed method provides practically useful and accurate recommendations of the split ratio. • We conduct extensive empirical evaluations on both synthetic and real-world datasets, highlighting the impact of principled data splitting on prediction interval length across diverse regression models. In the absence of prior work on optimal split selection, our experimental results underscore the importance of model-aware splitting strategies, demonstrating clear advantages over fixed split ratios while preserving finite-sample coverage guarantees. Collectively, these results provide the first principled solution to a fundamental challenge in split conformal prediction, enhancing both the efficiency and reliability of uncertainty quantification in practical applications. The remainder of the paper is organized as follows. Section 2 reviews split conformal prediction and its coverage guarantees. Section 3 presents our optimal splitting framework and theoretical results for various regression models. Section 4 reports empirical validation using synthetic and real-world datasets. We conclude in Section 5 with key findings and directions for future work.
2
Background on split conformal prediction
Conformal prediction provides a general, distribution-free framework for constructing predictive sets with finite-sample coverage guarantees [35, 42]. Let (Xi , Yi )n+1 i=1 be a random sample from an unknown distribution PXY , where Xi ∈ X and Yi ∈ Y. Based on the first n observations, the objective is to construct a measurable set C(Xn+1 ) ⊆ Y such that, for a prespecified miscoverage level α ∈ (0, 1), Pr{Yn+1 ∈ C(Xn+1 )} ≥ 1 − α, without imposing parametric assumptions on PXY . In this paper, we focus on the regression setting, where Y ∈ R and Y = f (X) + ϵ,
X ∼ PX ,
so that C(Xn+1 ) is an interval.
3
iid
ϵ ∼ Fϵ ,
2.1
Split conformal prediction
Split conformal prediction divides the observed sample (Xi , Yi )ni=1 into a proper training set, say D1 , of size n1 and a calibration set, say D2 , of size n2 = n − n1 . A regression estimator fˆ1 ≡ fˆn1 is fitted only on the proper training data D1 , while the calibration sample is used to compute conformity scores. The sample-splitting step ensures statistical independence between the calibration residuals, which is essential for achieving a finite-sample coverage guarantee [25, 29]. The choice of (n1 , n2 ) determines the trade-off between estimator accuracy (improved by a larger n1 ) and coverage accuracy (improved by a larger n2 ). In later sections, we analyze this trade-off through expansions of the expected length of the prediction interval and quantify its potential impact on the coverage error, which must be balanced for optimal performance of the split conformal prediction method.
2.2
Symmetric conformal prediction
In the standard regression setting, conformity scores are taken as absolute residuals [2]: ϵ̂i = |Yi − fˆ1 (Xi )|,
i ∈ D2 .
Let q̂ ≡ q̂n2 ,1−α denote the ⌈(n2 + 1)(1 − α)⌉-th order statistic of {ϵ̂i : i ∈ D2 } where, for x ∈ R, ⌈x⌉ (and ⌊x⌋) denote(s) the smallest integer not less than x (and the largest integer not exceeding x, respectively). The symmetric conformal prediction interval is then defined as ˆ n+1 ) = fˆ1 (Xn+1 ) − q̂, fˆ1 (Xn+1 ) + q̂ . C(X Exchangeability yields the finite-sample guarantee ˆ n+1 ) ≤ 1 − α + 1 − α ≤ Pr Yn+1 ∈ C(X
1 . n2 + 1
In symmetric conformal prediction, the interval length is L̂ = 2q̂ and its asymptotic behavior depends on different factors like the distribution of ϵ, and the choice of (n1 , n2 ), which we will investigate in the later sections.
2.3
Asymmetric conformal prediction
Symmetric intervals may be inefficient for skewed error distributions. Asymmetric split conformal prediction allows distinct upper and lower miscoverage levels [10]. Let αL , αU ≥ 0 satisfy αL + αU = α ∈ (0, 1). Define the signed conformity scores ϵ̃i = Yi − fˆ1 (Xi ),
i ∈ D2 .
Let q̃l ≡ q̃n2 ,αL and q̃u ≡ q̃n2 ,1−αU denote the ⌊(n2 + 1)αL ⌋-quantile of {ϵ̃i : i ∈ D2 } and the ⌈(n2 + 1)(1 − αU )⌉-quantile of {ϵ̃i : i ∈ D2 } respectively. Then the asymmetric interval is ˜ n+1 ) = fˆ1 (Xn+1 ) + q̃l , fˆ1 (Xn+1 ) + q̃u . C(X Exchangeability again ensures that ˜ n+1 ) ≤ 1 − α + 1 − α ≤ Pr Yn+1 ∈ C(X 4
2 , n2 + 1
and the interval length is L̃ = q̃u − q̃l . Later we will see that in many cases, the use of asymmetric intervals can be advantageous over the symmetric ones. iid Under the model Y = f (X) + ϵ, with X ∼ PX and ϵ ∼ Fϵ , both symmetric and asymmetric split conformal intervals provide exact finite-sample coverage bounds, which primarily depend on the calibration sample size n2 through the general upper bound n2k+1 (with k = 1 for the symmetric case and k = 2 for the asymmetric case). On the other hand, the length of the conformal interval is a random quantity based on both the training set size n1 (through the estimator fˆ1 ) as well as the calibration set size n2 (through the sample quantiles based on conformal scores). As a result, once the optimal split for the interval length is determined, the coverage accuracy/error admits a very definite nonasymptotic bound without further work. Therefore, in our analysis, we mainly focus on determining the optimal lengths of the prediction intervals by explicitly characterizing the asymptotic behaviors of their expected length and the variance. The relative contributions of these terms determine how the sample should be split into the training and calibration sets. By developing suitable expansions of E(L̂) (or E(L̃)) and Var(L̂) (or Var(L̃)) for general (n1 , n2 ), in the next section we derive the optimal split proportion that minimizes the interval length while preserving valid coverage.
3
Main results
We begin by introducing some notation that will be used throughout the paper, followed by a description of the assumptions required for our analysis. Section 3.3 presents results on the expectation and variance of the lengths of symmetric and asymmetric conformal prediction intervals. In Section 3.4, we specialize the results to three regression settings, namely linear regression, nonparametric regression, and neural network regression, and derive specific optimal split proportions in each case. Finally, in Section 3.5, we describe a data-driven procedure for selecting the optimal split based on a subsampling technique.
3.1
Notation
For a function g on the real line, we denote its first and second derivatives by g ′ and g ′′ , respectively. The inverse of the function g is denoted by g −1 . Throughout the paper, C1 , C2 , . . . denote constants that are independent of the sample size n. For two positive sequences an and bn , we write an = O(bn ) if there exists a constant C such that an ≤ Cbn for all n. Similarly, an = o(bn ) indicates that an /bn → 0, as n → ∞. Furthermore, by an ≍ bn , we denote that an = O(bn ) and bn = O(an ), simultaneously. Recall that the floor and ceiling functions, denoted by ⌊x⌋ and ⌈x⌉, represent the greatest integer less than or equal to x, and the smallest integer greater than or equal to x, respectively. For an estimator fˆn of f based on sample size n, define the estimation error at a point x by γ̂n (x) = fˆn (x) − f (x), which measures the discrepancy between the estimator fˆn (x) and the true regression function f (x). For notational conveniences, we simply write the powers (γ̂n (x))k as γ̂nk (x) for any k ≥ 2. We denote by E1 , Var1 , and Cov1 the conditional expectation, variance, and covariance, respectively, given the proper training set D1 . For example, the L2 error can be written as Z E |fˆn (x) − f (x)|2 PX (dx) = EE1 (γ̂n2 (X)).
5
Similarly, the conditional variance is written as EVar1 (γ̂n (X)) = EE1 (γ̂n2 (X)) − E{E1 (γ̂n (X))}2 . Furthermore, for notational brevity, we will write γ̂n1 = fˆn1 − f as γ̂1 .
3.2
Conditions
For the symmetric split conformal prediction interval, define mα = (n2 + 1)(1 − α) ,
µα =
mα , n2 + 1
σα2 =
µα (1 − µα ) . n2 + 2
For the asymmetric split conformal prediction interval, define m1 = (n2 + 1)αL , m2 = (n2 + 1)(1 − αU ) , and µi =
mi , n2 + 1
σi2 =
together with ρ12 =
µi (1 − µi ) , n2 + 2
i = 1, 2,
µ1 (1 − µ2 ) . n2 + 2
Let S denote the population conformity score, with distribution function FS . Throughout the paper, ( |ϵ|, for symmetric split conformal prediction, S= ϵ, for asymmetric split conformal prediction. Furthermore, let δ ∈ (0, 1] be a constant independent of n. We make the following assumptions. (A.1) For some integer K ≥ 1, suppose there exists a positive, strictly increasing sequence of constants {βk }K k=1 , independent of n, such that Z E |fˆn (x) − f (x)|k PX (dx) ≍ n−βk . (A.2)r Suppose there exists two sequences defined by a1n = δ −1 n−γ , for some γ ∈ (0, (βk −βr )/k), k > p r and a2n = δ −1 log(n)/n, such that for each relevant quantile level µ, −1 −1 ′ ′ δ< inf FS FS (µ + x) + y ≤ sup FS FS (µ + x) + y < δ −1 , |y|≤a1n ,|x|≤a2n
|y|≤a1n ,|x|≤a2n
sup |y|≤a1n ,|x|≤a2n
and sup |y|≤a1n ,|x|≤a2n
FS′′ FS−1 (µ + x) + y < δ −1 ,
FS′′ FS−1 (µ + x) + y − FS′′ FS−1 (µ) = o(1), as n → ∞.
For Theorem 3.1, corresponding to the symmetric interval, this condition is imposed with FS = F|ϵ| and µ = µα . For Theorem 3.5, corresponding to the asymmetric interval, it is imposed with FS = Fϵ and µ ∈ {µ1 , µ2 }. 6
(A.3) There exists a constant βϵ > 0, independent of n, such that for all x ≥ δ −1 , Fϵ (−x) + (1 − Fϵ (x)) ≤ δ −1 |x|−βϵ . Assumption (A.1) imposes polynomial moment bounds on the estimation error of the R regression estimator fˆn . Specifically, it requires that for some integer K ≥ 1, the Lk -risk E |fˆn (x) − f (x)|k PX (dx) decays at the rate n−βk for k = 1, . . . , K. Such conditions are standard in nonparametric regression and statistical learning theory, and are satisfied by a broad range of estimators, including kernel smoothers, local polynomial estimators, and certain neural network regressors; see, for example, [14, 18]. In the context of conformal inference, moment bounds of this form ensure an accurate approximation of the residual distribution, and have appeared in recent theoretical analyses of predictive inference [25]. Assumption (A.2)r imposes smoothness and local stability conditions on the noise distribution function FS (where S is either |ϵ| or ϵ). Specifically, it requires that the first and second derivatives of FS remain uniformly bounded away from zero and infinity within a shrinking neighborhood. These conditions ensure that the quantile function does not become excessively flat or steep. Such conditions are classical in asymptotic quantile theory and in deriving Bahadur-type expansions [34]. The final requirement in Assumption (A.2)r , which controls the difference between second derivatives evaluated at nearby points, is essentially a local Lipschitz condition on FS′′ . Such derivative-boundedness and local smoothness conditions are standard in the analysis of quantile estimators and residual-based inference; see, for example, [40], where smoothness of the noise density plays a key role in obtaining uniform expansions of sample quantiles. Assumption (A.3) imposes a polynomial tail condition on the noise distribution Fϵ . This assumption allows for moderately heavy-tailed noise and is substantially weaker than sub-Gaussian or sub-exponential assumptions. Polynomial tail bounds of this form ensure that the contribution of extreme residuals is asymptotically negligible, thereby facilitating higher-order asymptotic approximations of quantile-based statistics. Similar tail conditions are common in the asymptotic theory of order statistics; see, for example, [7] and [37]. In the present setting, Assumption (A.3) is crucial for guaranteeing that contributions from extreme residuals are asymptotically negligible in the expansion of the conformal quantile, so that the approximation error is dominated by the smooth interior behavior described in Assumption (A.2)r . Taken together, assumptions (A.1)–(A.3) define a flexible and realistic framework under which the regression estimator converges at a polynomial rate, the residual quantile function is locally smooth, and the noise distribution exhibits only moderate polynomial tail decay. These conditions are mild and broadly consistent with regularity assumptions commonly used in nonparametric regression and quantile inference. They ensure that the analytic approximation of the conformal prediction interval length holds uniformly and that higher-order remainder terms remain negligible.
3.3
Main results
In this subsection, we present the main theoretical results of the paper. Under Assumptions (A.1)– (A.3), we derive expressions for the expectation and variance of the length of the split conformal prediction interval and use these results to characterize the optimal data-splitting strategy. The following theorem provides the mean and variance of the length of the symmetric split conformal prediction interval. h Theorem 3.1. Let βϵ > 2/n2 and α ∈ (βϵ + 1)/(βϵ (n2 + 1)), 1 − 1/(βϵ (n2 + 1)) .
7
(a) Suppose Assumptions (A.1), (A.2)1 and (A.3) hold, then we have, 1 1 −1 −1 2 −1 ′′ . E L̂ = 2F|ϵ| (µα ) + 2H|ϵ| F|ϵ| (µα ) EE1 (γ̂1 (X)) + σα F|ϵ| (µα ) + o β1 + n2 n1 Furthermore, under Assumption (A.2)2 , n n o2 o2 ′ 1 1 −1 −1 2 VarE1 (γ̂1 (X)) + 4σα F|ϵ| (µα ) + o β2 + Var L̂ = 4 H|ϵ| F|ϵ| (µα ) , n2 n1 ′ (a). where H|ϵ| (a) = −{Fϵ′ (a) − Fϵ′ (−a)}/F|ϵ|
(b) If the distribution of ϵ is symmetric and Assumption (A.2)2 holds, then 1 1 −1 −1 2 2 −1 ′′ . E L̂ = 2F|ϵ| (µα ) + 2Gϵ F|ϵ| (µα ) EE1 (γ̂1 (X)) + σα F|ϵ| (µα ) + o β2 + n2 n1 Furthermore, under Assumption (A.2)4 , n o2 n o2 ′ 1 1 −1 −1 Var L̂ = 4 Gϵ F|ϵ| (µα ) VarE1 (γ̂12 (X)) + 4σα2 F|ϵ| (µα ) + o β4 + , n2 n1 where Gϵ (a) = −Fϵ′′ (a)/(2Fϵ′ (a)). −1 Remark 3.2. The leading term 2F|ϵ| (µα ) in the expansion of E(L̂) corresponds to the expected length obtained when the regression function is known exactly. Since
µα =
⌈(n2 + 1)(1 − α)⌉ ≥ 1 − α, n2 + 1
−1 and F|ϵ| is an increasing function, this quantity is always at least as large as the oracle interval length −1 (1 − α). Thus, even in the absence of estimation error, the split conformal prediction interval 2F|ϵ| exhibits a positive finite-sample bias arising from the calibration quantile. The remaining terms in Theorem 3.1 quantify the additional increase in interval length due to estimating the regression function. Consequently, minimizing these higher-order terms yields the shortest attainable split conformal prediction interval while preserving the finite-sample coverage guarantee. −1 Remark 3.3. In the case of a symmetric noise distribution, the leading term 2F|ϵ| (µα ) depends only on the calibration sample size through the empirical quantile level µα and is independent of the proper training sample size. Since its dependence on n2 contributes only to the constant of the calibration error of order n−1 2 , it does not affect the asymptotically optimal order of the training– calibration split. Consequently, the optimal allocation of samples is determined by balancing the regression estimation error, EE1 (γ̂12 (X)) and VarE1 (γ̂12 (X)),
which decrease as the proper training sample size n1 increases, against the overall calibration error of order n−1 2 , which decreases as the calibration sample size n2 increases. For example, suppose EE1 (γ̂12 (X)) ≍ n−β 1 ,
VarE1 (γ̂12 (X)) ≍ n−2β , 1
8
for some β > 0. Then a bias-optimal split satisfies −1 n−β 1 ≍ n2 ,
whereas an MSE-optimal split satisfies n−2β ≍ n−1 1 2 . The corresponding values of β for several commonly used regression models are derived in Section 3.4. Remark 3.4. The quantities EE1 (γ̂1 (X)) and VarE1 (γ̂1 (X)) appearing in Theorem 3.1(a) may 1 2 converge faster than the upper bounds n−β and n−β implied by Assumption (A.1). Thus, the 1 1 stated rates should be interpreted as worst-case controls rather than exact orders. Furthermore, the remainder term can be made arbitrarily small by imposing higher-order smoothness assumptions on the error distribution. In particular, if Assumption (A.2)r holds with r = n0 for some positive integer −βn n0 , then the remainder becomes o(n1 0 + n−1 2 ), showing that increasingly accurate asymptotic expansions can be obtained by assuming additional differentiability of the error distribution. The next theorem gives the mean and variance of the length of the asymmetric split conformal prediction interval. Theorem 3.5. Suppose βϵ ≥ 2/(n2 −1), αL ∈ ((βϵ +1)/(βϵ (n+1)), 1/2], and αU ∈ [(βϵ +1)/(βϵ (n+ 1)), 1/2), with αL + αU = α. Then, under Assumptions (A.1), (A.2)2 and (A.3), we have o n E L̃ = Fϵ−1 (µ2 ) − Fϵ−1 (µ1 ) n o + Gϵ Fϵ−1 (µ2 ) − Gϵ Fϵ−1 (µ1 ) EVar1 γ̂1 (X) n σ2 o σ12 −1 ′′ 1 1 −1 ′′ 2 + F (µ2 ) − F (µ1 ) + o β2 + . 2 ϵ 2 ϵ n2 n1 Furthermore, if Assumption (A.2)4 holds, then o2 n Var L̃ = Gϵ Fϵ−1 (µ2 ) − Gϵ Fϵ−1 (µ1 ) VarVar1 (γ̂1 (X)) n o2 n o2 ′ ′ + σ22 Fϵ−1 (µ2 ) + σ12 Fϵ−1 (µ1 ) 1 1 −1 ′ −1 ′ − 2ρ12 Fϵ (µ2 )Fϵ (µ1 ) + o β4 + , n2 n1 where Gϵ (a) = −Fϵ′′ (a)/(2Fϵ′ (a)). Remark 3.6. Theorem 3.5 shows that the approximation error for the expected length of the split conformal prediction interval can be substantially smaller for asymmetric intervals than for the classical symmetric interval. The key distinction is that the leading term in Theorem 3.5 depends on EVar1 (γ̂1 (X)), whereas the corresponding term in Theorem 3.1(b) depends on EE1 (γ̂12 (X)). Since Var1 (γ̂1 (X)) = E1 (γ̂12 (X)) − {E1 (γ̂1 (X))}2 , we always have
EVar1 (γ̂1 (X)) ≤ EE1 (γ̂12 (X)), 9
with strict inequality whenever the conditional bias of the regression estimator is nonzero. Consequently, the leading approximation error for asymmetric conformal prediction may be of strictly smaller order than that of the symmetric interval. Moreover, unlike the symmetric case, this improvement does not require the noise distribution to be symmetric, making asymmetric conformal prediction theoretically advantageous under considerably weaker assumptions. Similar conclusions apply to the variance expansion.
3.4
Examples
In this subsection, we derive the expressions that appeared in the main theoretical results, which enable practical determination of the optimal split ratio between the proper training and calibration samples. For illustration, we focus on a symmetric conformal prediction interval with a symmetric error distribution. We consider three regression settings: linear regression, nonparametric regression via the Nadaraya-Watson kernel estimator, and fully connected neural network regression. These examples illustrate how the general theory applies across different regression paradigms. Since 2 2 σα2 ≍ n−1 2 always holds, it suffices to characterize EE1 (γ̂1 (X)) and VarE1 (γ̂1 (X)), equivalently the expectation and variance of the L2 errors, as functions of the (proper) training size n1 for each model. 3.4.1
Linear regression
We begin by considering a random effects model with the linear regression function f (x) = x⊤ β, where x ∈ Rd and β ∈ Rd are the covariates and regression coefficients, respectively. For illustrative 1 purposes, let {(Xi , Yi )}ni=1 denote a collection of n1 independent and identically distributed samples. iid
We consider a Gaussian design with Xi ∼ N (0, Id ), where Id is the d-dimensional identity matrix, iid
and errors ϵi ∼ N (0, 1). Let X denote the design matrix and Y the response vector, then the ordinary least squares linear regression estimate of the regression function is given by fˆ(x) = x⊤ β̂, where β̂ = (X ⊤ X)−1 X ⊤ Y and we have EE1 (γ̂12 (X)) = and
d 1 ≍ , n1 − d − 1 n1
h i h i 1 VarE1 (γ̂12 (X)) = 2E tr(X ⊤ X)−2 + Var tr(X ⊤ X)−1 ≍ 2 , n1
provided n1 > d + 3. Then, in ordinary least squares linear regression, a bias-optimal split requires n1 ≍ n2 , that is, the training and calibration sample sizes are of the same order. For an MSE-optimal √ split, one requires n1 ≍ n2 . 3.4.2
Nonparametric regression
We consider a nonparametric regression framework in which the goal is to estimate the regression 1 function f (x) = E[Y | X = x] from n1 iid samples {(Xi , Yi )}ni=1 . We assume that the covariates d 1 Xi ∈ [−a, a] , for some a > 0 are uniformly distributed, and the noise terms {ϵi }ni=1 are iid with iid
ϵi ∼ N (0, 1). For estimation, we employ the Nadaraya-Watson kernel estimator, a local-constant regression method that smooths the data using a bandwidth parameter h > 0 and a kernel function
10
K : Rd → R+ that satisfies
R
K(u) du = 1. The estimator at a point x ∈ Rd is defined as Pn1 x−Xi Yi i=1 K h , fˆh,n1 (x) = P n1 x−Xi i=1 K h
which can be viewed as a locally weighted average of the observed responses Yi , with larger weights assigned to observations where the covariates Xi are closer to x. For illustrative purposes, we assume that the kernel K is bounded, symmetric, and continuously differentiable with compact support. −1/(4+d) Finally, if we choose the bandwidth h ≍ n1 , standard calculation shows EE1 (γ̂12 (X)) ≍
1 4/(4+d) n1
and VarE1 (γ̂12 (X)) ≍
,
1 8/(4+d) n1
,
provided the regression function has bounded first and second derivatives. Then, in the NadarayaWatson kernel regression, the theoretically justified choice for the bias-optimal data splitting is 1+d/4 n1 ≍ n 2 , implying that the training sample should be considerably larger than the calibration 1/2+d/8 sample and n1 ≍ n2 for the MSE-optimal choice. 3.4.3
Neural network regression
For the final example, we consider a fully connected neural network with the Rectified Linear Unit (ReLU) activation function. In this setting, exact expressions for the expectation and variance of the L2 -error are not yet available in the literature. However, bounds on the expected L2 -error have been derived in [21], under certain structural assumptions on the regression function. To state these results, we rely on the following two definitions. Definition 3.7 ((p, C)-smoothness, [21]). Let p = q + s for some q ∈ N0 and 0 < s ≤ 1. A function m : Rd → R is called (p, C)-smooth if, for every multi-index α = (α1 , . . . , αd ) ∈ Nd0 with Pd j=1 αj = q, the partial derivative ∂qm ∂xα1 1 · · · ∂xαd d exists and satisfies
∂qm
∂qm
(x) − α1 (z) ∂xα1 1 · · · ∂xαd d ∂x1 · · · ∂xαd d
≤ C∥x − z∥s ,
for all x, z ∈ Rd , where ∥ · ∥ denotes the Euclidean norm. With the definition of (p, C)-smoothness, we define the class of hierarchical composition models as follows. Definition 3.8 (Hierarchical composition models, [21]). For l = 1 and smoothness constraint P ⊆ (0, ∞) × N, the space of hierarchical composition models is defined as H(1, P) := h : Rd → R : h(a) = m aπ(1) , . . . , aπ(K) , where m : RK → R is (p, C)-smooth for some (p, K) ∈ P and π : {1, . . . , K} → {1, . . . , d} . 11
For l > 1, we recursively define := H(l, P) h : Rd → R : h(x) = m f1 (a), . . . , fK (a) , where m : RK → R is (p, C)-smooth for some (p, K) ∈ P and fi ∈ H(l − 1, P) . Let the corresponding regression function f belong to the class H(l, p) for some l ∈ N and P ⊆ (0, ∞) × N. Then according to Theorem 1 in [21], and under some regularity conditions, we have − 2p EE1 (γ̂12 (X)) = O (log(n))6 max n 2p+K . (p,K)∈P
Although the exact order is not available, for a bias-optimal split, we can use the approximate rela1+K̄/2p̄ tion n1 ≍ n2 , where (p̄, K̄) ∈ P is chosen such that (p̄, K̄) = arg min(p,K)∈P p/K. This implies that the relative size of the proper training set to the calibration set depends on the smoothness and order constraint P.
3.5
Data based optimal choice of training ratio
When the sample size n is very large, it may be difficult to determine the optimal split, as doing so requires repeatedly training the model under different splits, which may be computationally infeasible in practice. However, using our results, one can employ a data-driven procedure to select the optimal split via subsampling (cf. [19]). Based on the preceding results, we obtain that the optimal split satisfies n1 ≍ nb2 for some known constant b > 0. Let ĉℓ denote the empirical estimate of the ratio of the training sample size to the total sample size, computed from a subsample of size ℓ(< n). The corresponding estimate ĉn for the full sample size n is then determined as the solution to the equation x=
nb−1 ĉℓ (1 − x)b . ℓb−1 (1 − ĉℓ )b
(3.1)
The following algorithm provides a way to approximate the ratio ĉℓ where ℓ < n. Algorithm Description. Given α ∈ [1/(n2 + 1), 1), and ℓ(< n), we choose ĉℓ as follows: 1. Divide the dataset randomly into K disjoint subsets {Qk }K k=1 , each containing ℓ samples, such that Kℓ ≤ n. 2. Choose a grid Gg = {p1 , . . . , pg } ∈ (0, 1)g of size g (e.g., G6 = {0.4, 0.5, . . . , 0.9}). 3. For each k = 1, 2, . . . , K: (a) For each candidate proportion pi ∈ Gg , i = 1, . . . , g, further split the subsets Qk into a proper training set Tk,i of size ⌊ℓpi ⌋ and a calibration set Ck,i of size ℓ − ⌊ℓpi ⌋. (b) Using Tk,i and Ck,i , compute prediction interval length L̂(Tk,i , Ck,i , pi ) and select the target proportion p̂k = arg minpi ∈Gg L̂(Tk,i , Ck,i , pi ). P 4. Set ĉℓ = K1 K k=1 p̂k 5. (Optional refinement) If greater precision is required, repeat the Steps 2–5 using a finer grid centered around the current estimate ĉℓ , and update ĉℓ accordingly. 12
Remark 3.9. In practice, one can choose ℓ such that K = 5 or 10 or more based on the sample sizes and model computation time. In this way, we need to train the model only on a dataset of size ℓp, for p ∈ G, for which we can have ℓp << n1 . However, it is important to note that, for the algorithm to function correctly, at least one value in Gg must be less than p∗ such that α > 1/(⌈ℓ(1 − p∗ )⌉ + 1). Remark 3.10. To the best of our knowledge, existing empirical practice has largely treated the training-calibration split in split conformal prediction as fixed rather than tunable. A common convention in the literature is to use equal-sized proper training and calibration sets, as adopted in empirical evaluations of split conformal methods with ridge regression, random forests, and neural networks in [31]. This convention effectively promotes a 1:1 split as a default choice, despite the lack of theoretical or empirical justification for its optimality in terms of interval length. Related recent work has explicitly targeted the problem of prediction set length, but from a different perspective. In particular, recent work formulates length optimization as a constrained problem that minimizes expected prediction set size subject to coverage constraints, and develops a minimax-based procedure that learns covariate-dependent thresholds to construct shorter sets [20]. Their framework is implemented on top of diverse predictive models, including linear models and deep neural network architectures such as ResNet and large language model backbones. However, across these experiments, the allocation between training and calibration data is not treated as a parameter to be optimized and is fixed within each setup, irrespective of the underlying model. In contrast, our results indicate that a 1:1 split is not generally optimal. The allocation of samples between training and calibration interacts with the underlying learning algorithm, as different models induce different trade-offs between estimation accuracy and calibration precision. Consequently, the split ratio should be regarded as an algorithm-dependent design parameter, and selecting it adaptively can yield shorter prediction intervals while preserving nominal coverage.
4
Experimental results
In this section, we evaluate our theoretical results using both synthetic and real-world datasets across three regression settings: linear regression, nonparametric regression, and neural networks. In synthetic datasets, each experiment involves n + 1 samples, consisting of n samples used for training and calibration, and one sample for testing. For each underlying model, we vary the training sample size n and randomly partition the data into training and calibration sets. To ensure reliability, all experiments are repeated 1,000 times for each training-calibration split ratio. We employ split conformal prediction to construct prediction intervals with a nominal coverage level of 1 − α = 0.9 and report the average prediction interval length.
4.1
Linear regression
We examine the following linear regression model: yi = x⊤ i β + ϵi ,
i = 1, . . . , n,
(4.2)
where yi ∈ R is the response variable, xi ∈ [0, 1]5 denotes the regressors, β ∈ R5 represents the regression coefficients, and ϵi ∈ R is the error term. Using equation (4.2), we generate synthetic datasets with varying sample sizes {100, 300, 500, 800, 2000}. The regressors xi are sampled independently from a uniform distribution on [0, 1]5 . We consider three types of error distributions: • Normal distribution: ϵi ∼ N (0, 1) 13
• Student’s t-distribution with 5 degrees of freedom (variance-normalized): t(5) ϵi ∼ p 5/3 • Lognormal distribution (standardized to zero mean and unit variance): ϵi ∼
Lognormal(0, 1) − exp(0.5) p (exp(1) − 1) exp(1)
We repeat the experiments 1,000 times for each error setting and sample size and report the average prediction interval lengths in the first row of Figure 1. As shown in Section 3.4.1, for linear regression models a bias-optimal split requires the training and calibration set sizes to satisfy n1 ≍ n2 . That is, the training set size should be approximately equal to the calibration set size. Equivalently, this corresponds to setting b = 1 in Equation 3.1. Solving the equation for b = 1 yields x = ĉℓ , which is consistent with the numerical results in the first row of Figure 1. As the total sample size increases, the optimal prediction interval length continues to occur when the training and calibration sets are approximately equal in size.
4.2
Nonparametric regression
In this section, we investigate two settings of a nonparametric regression model of the form yi,j = fj (xi ) + λj ϵi ,
i = 1, . . . , n,
j = 1, 2, (1)
(5)
where yi,j ∈ R is the response variable, fj is the regression function, xi = (xi , . . . , xi )⊤ ∈ R5 is a five-dimensional predictor which is uniformly distributed over the region [0, 1]5 independently, ϵi is the independent error term, which is also independent of xi as defined in Subsection 4.1, and λ1 , λ2 ≥ 0 are fixed scaling parameters. For reasons of comparability, the scaling parameters λj are chosen to match 20% of the typical range of variation of fj (x) under the distribution of x. This range is chosen via the interquartile range (IQR) of fj (x) computed from many simulated draws of x, so that the noise level is scaled relative to the signal magnitude. These scaling parameters are set to λ1 = 0.340 and λ2 = 0.262. The regression functions under consideration are: • Function f1 : f1 (x) = 3x(1) + tan(x(2) ) + (x(3) )3 + log(x(4) + 0.1) +
p x(5) + 0.1
• Function f2 : f2 (x) = exp(∥x∥), where x = (x(1) , . . . , x(5) )⊤ . Following the experimental setup in the linear regression setting, we generate synthetic datasets with sample sizes {100, 200, 500, 800, 2000}. For each sample size, we consider three types of error distributions: normal, lognormal, and Student’s t-distribution. Each experiment is repeated 1,000 times to ensure statistical reliability. The average prediction interval lengths computed via split conformal prediction for the two functions are shown in the second and third rows of Figure 1. From our theoretical analysis of Nadaraya-Watson kernel regression, the bias-optimal data split 14
Prediction Interval Length
Comparison of Regression Methods n = 100
6
n = 300 5
5
3.5
3.5
3.0
3.0
3
3
2.5
2.5
2.5
2
2
2.0
2.0
2.0
60
n1/n (%)
80
100
40
60
n1/n (%)
80
1.5 1.0 40
60
n1/n (%)
80
100
40
60
n1/n (%)
80
2.5
1.5 1.0 0.5 40
60
n1/n (%)
80
7 6 5 4 3 2 1
100
40
60
n1/n (%)
80
n1/n (%)
80
80
100
40
60
80
1.0
1.0
n1/n (%)
80
40
60
80
4
1.00
2
0.75 60
n1/n (%)
80
0.50 100
40
60
n1/n (%)
80
100
n1/n (%) Neural Network Regression for function f2 1.8 1.6 1.4 1.2 1.0 0.8 0.6
1.25
100
80
0.6 100
40
80
1.1 1.0 0.9 0.8 0.7 0.6 0.5 100
40
60
n1/n (%)
40
80
100
60
80
100
60
80
100
0.8 40
60
n1/n (%)
80
0.6 100
1.6
1.6
1.4
40
n1/n (%)
1.2 1.0 0.8
0.8
0.6
0.6
n1/n (%)
60
n1/n (%)
1.0
1.8
100
100
1.2
1.0
80
80
1.4
1.2
60
60
n1/n (%)
1.6
1.4
40
100
1.8
2.0
1.0 60
60
n1/n (%)
2.5
1.5
1.50
40
100
1.5
40
40
0.6 80
80
0.8
0.8
60
60
n1/n (%)
1.0
1.0
40
40
1.2
1.2
2.0
1.75
6
0.6 100
100
1.4
0.8
1.5
100
80
1.0
n1/n (%) Neural Network Regression for function f1
0.5 100
60
n1/n (%)
1.2
n1/n (%) Nonparametric Regression for function f2
100
40
1.4
2.00 8
100
1.6
2.5
2.0
60
60
n1/n (%) Nonparametric Regression for function f1
3.0
2.5
40
40
1.8 1.6 1.4 1.2 1.0 0.8 0.6
1.8 1.6 1.4 1.2 1.0 0.8 0.6
2.0
100
2.25 2.00 1.75 1.50 1.25 1.00 0.75
2.25 2.00 1.75 1.50 1.25 1.00 0.75
2.0
0.5
Prediction Interval Length
3.5
Student's t-distribution
n = 2000
4.0
3.0
2.5
Prediction Interval Length
n = 800
4.0
4
3.0
Prediction Interval Length
Linear Regression n = 500
Lognormal
4
40
Prediction Interval Length
4.0
Normal
40
60
n1/n (%)
80
100
40
n1/n (%)
Figure 1: Prediction interval lengths for synthetic datasets as a function of the data split ratio. Solid lines represent the experimental results. Dashed lines indicate the locations of the minimum prediction interval length. Dash-dot lines correspond to the true prediction interval lengths (i.e., the oracle prediction lengths). 15
1+d/4
satisfies n1 ≍ n2 , where d denotes the covariate dimension. In our experiments, d = 5, which 9/4 yields n1 ≍ n2 . This scaling indicates that nonparametric regression requires a larger training set than calibration set, since accurate nonparametric estimation relies primarily on having more training data. Under this relationship, the corresponding ratio parameter in Equation 3.1 is b = 9/4. Substituting this value into Equation 3.1 and considering the setting n = 2000 and ℓ = 300, together with the empirically observed minimum prediction interval length at ĉℓ ≈ 0.75, yields x ≈ 0.91. This theoretical value aligns closely with the experimental findings reported in the last columns of the second and third rows of Figure 1. Overall, these results provide empirical support for the proposed algorithm for selecting the split ratio and validate the practical relevance of Equation 3.1 in the nonparametric regression setting.
4.3
Neural network regression
In this section, we evaluate our theoretical results on neural network regression using synthetic and real-world datasets. 4.3.1
Synthetic dataset
For the synthetic dataset, we again consider the setting in Subsection 4.2: yi,j = fj (xi ) + λj ϵi ,
i = 1, . . . , n,
j = 1, 2.
We generate synthetic datasets for sample sizes {100, 300, 500, 800, 1000} and consider three types of error distributions: normal, lognormal, and Student’s t-distribution. Each configuration is evaluated over 1,000 repetitions. We report the average prediction interval lengths for each setup, with the results displayed in the fourth and fifth rows of Figure 1. In our simulations, the neural network consists of one hidden layer for sample sizes of {100, 300} and two hidden layers for sample sizes of {500, 800, 2000}. The network is trained for a maximum of 100 epochs in all configurations. From our theoretical analysis of neural networks, a bias-optimal split approximately satisfies n1 ≍ 1+K̄/2p̄ 7/2 n2 . In our setting, p̄ = 1 and K̄ = 5 for the functions f1 and f2 , which gives n1 ≍ n2 . As introduced in Section 3.4.3, (p̄, K̄) denote the smoothness and order parameters in P selected as the minimizer of p/K. Similar to the nonparametric regression case, this scaling indicates that neural networks also require a larger training set than a calibration set. Under this relationship, the corresponding ratio parameter in Equation 3.1 is b = 7/2. Substituting this value into Equation 3.1 with n = 2000 and ℓ = 100, and using the empirically observed minimum prediction interval length at ĉℓ ≈ 0.6, yields x ≈ 0.93 (similar results are obtained if we consider higher values of ℓ). This theoretical value closely matches the experimental results reported in the last columns of the fourth and fifth rows of Figure 1. These findings further support the proposed algorithm for selecting the split ratio. 4.3.2
Concrete compressive strength dataset
To evaluate the empirical performance of the proposed algorithm on a benchmark dataset, we consider the Concrete Compressive Strength dataset [1]. The neural network used in our experiments has two hidden layers and is trained for up to 100 epochs. The dataset is partitioned into training and test sets using an 80:20 split, yielding 824 training and 206 test observations. The training portion is further divided into training and calibration sets for conformal prediction. To apply Equation 3.1 and the proposed procedure for selecting a split ratio that minimizes prediction interval length, we
16
set p̄ = 1 and K̄ = 8 (see Section 3.4.3). Although these parameters can be chosen in a datadependent way, we only choose them heuristically. Under these choices, n1 ≍ n52 , corresponding to b = 5. To implement Equation 3.1, the training data are partitioned into four subsets, giving ℓ = 206 in Equation 3.1. Experiments are conducted over a range of split ratios within each subset, and split conformal prediction is applied to the test data. Results are averaged across the four subsets, as shown in the left panel of Figure 2. The estimate ĉℓ = 0.66 is obtained. Substituting these values into Equation 3.1 yields x ≈ 0.88. This theoretical value is consistent with the empirical results shown in the right panel of Figure 2, which are based on the full training dataset. This agreement provides further support for the proposed split-ratio selection procedure.
Average over 4 subsets of the training dataset (206 samples per subset)
Full training dataset (824 samples) 42
Prediction Interval Length
Prediction Interval Length
70 65 60 55 50 45
50
60
70
n1/n (%)
80
40 38 36 34
90
50
60
70
n1/n (%)
80
90
Figure 2: Prediction interval length for different training samples on the concrete compressive strength dataset.
5
Discussion and future directions
While our experiments are conducted on datasets of moderate size, where one could, in principle, evaluate multiple train-calibration splits to search for shorter prediction intervals, such an approach is often impractical in real-world deployments. In many safety-critical and high-stakes domains, such as healthcare, autonomous systems, and financial risk management, datasets often contain a large number of samples, making model training computationally demanding, time-consuming, and subject to operational or regulatory constraints. Repeated retraining across numerous candidate split ratios is therefore inefficient and, in some cases, infeasible. Additionally, in certain applications, data acquisition itself is costly or limited, further reducing the viability of brute-force split selection. Our method instead provides a principled approach to selecting the data split with the objective of minimizing prediction interval length without requiring repeated retraining on a large dataset. To the best of our knowledge, existing work that aims to reduce conformal prediction interval length typically relies on random or heuristic data splits rather than explicitly optimizing the traincalibration partition. This highlights the practical relevance and novelty of our approach.
17
6
Proof of the main results
First, we present the proof of Theorem 3.5 and then move on to prove Theorem 3.1. h Proof of Theorem 3.5. Fix α2 ∈ (βϵ + 1)/βϵ (n2 + 1)), 1 − 1/βϵ (n2 + 1)) . Let F̃ be the conditional CDF of ϵ̃i = ϵi − γ̂(Xi ), i ∈ D2 , given D1 . Then, E(q̃u ) = E(E(q̃u | D1 )), Var(q̃u ) = Var(E(q̃u | D1 )) + E(Var(q̃u | D1 )), where q̃u is the m2 = ⌈(n2 + 1)(1 − α2 )⌉-th quantile of {ϵ̃i }i∈D2 (same argument follows for m1 = ⌊(n2 +1)α1 ⌋-th order statistic of {ϵ̃i }i∈D2 ). Let U1 , . . . , Un2 be a set of uniformly distributed random variables on (0, 1). Define U(m2 ) as the mth order statistic of U1 , . . . , Un2 , with mean µ2 and variance p 2 σ22 . Recall that, we have a2n2 = C log(n2 )/n2 , for some constant C > 0. Then, using Taylor’s theorem, we have E(q̃u | D1 ) = E1 F̃ −1 (U(m2 ) ) = E1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | ≤ a2n2 ) + E1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | > a2n2 ) ′ = E1 F̃ −1 (µ2 ) + (U(m2 ) − µ2 )F̃ −1 (µ2 ) 1 2 −1′′ (µ2 ) 1(|U(m2 ) − µ2 | ≤ a2n2 ) + (U(m2 ) − µ2 ) F̃ 2 1 ′′ ′′ + E1 (U(m2 ) − µ2 )2 F̃ −1 (µ∗ ) − F̃ −1 (µ2 ) 1(|U(m2 ) − µ2 | ≤ a2n2 ) 2 + E1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | > a2n2 )
(6.3)
σ22 −1′′ F̃ (µ2 ) + R1n , 2 where µ∗ is some point between U(m2 ) and µ2 and R1n is the remainder term defined as ′ R1n = −E1 F̃ −1 (µ2 ) + (U(m2 ) − µ2 )F̃ −1 (µ2 ) 1 2 −1′′ + (U(m2 ) − µ2 ) F̃ (µ2 ) 1(|U(m2 ) − µ2 | > a2n2 ) 2 1 ′′ ′′ + E1 (U(m2 ) − µ2 )2 F̃ −1 (µ∗ ) − F̃ −1 (µ2 ) 1(|U(m2 ) − µ2 | ≤ a2n2 ) 2 + E1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | > a2n2 ). = F̃ −1 (µ2 ) + 0 +
Then, we obtain |R1n | ≤ |F̃ −1 (µ2 )|P(|U(m2 ) − µ2 | > a2n2 ) + |F̃ −1 (µ2 )| E|U(m2 ) − µ2 |1(|U(m2 ) − µ2 | > a2n2 ) ′
+ |F̃ −1 (µ2 )| E|U(m2 ) − µ2 |2 1(|U(m2 ) − µ2 | > a2n2 ) ′′
+ σ22
sup |F̃
−1′′
(µ2 + x) − F̃
−1′′
(µ2 )|
|x|≤a2n2
+ E1 |F̃ −1 (U(m2 ) )|1(|U(m2 ) − µ2 | > a2n2 ). 18
(6.4)
A standard Hoeffding bound for Beta random variables gives, P(|U(m2 ) − µ2 | > a2n2 ) ≤ 2 exp(−2n2 a22n2 ) = o and later we show that,
1 , n2
F̃ −1 (µ2 ) = Fϵ−1 (µ2 ) + R2n ,
(6.5)
with E|R2n | = o(1). Therefore, by assumptions (A.1), (A.2)2 and (A.3), we get that E|R1n | = o(1/n2 ). Next, we show that the quantiles (conditional on D1 ) of ϵ̃i are close to the quantiles of ϵi . First, note that for some fixed a and h, by Taylor’s theorem F̃ (a + h) = E1 Fϵ (a + γ̂1 (X) + h) n o = Fϵ (a) + Fϵ′ (a) E1 (γ̂1 (X)) + h n o 1 + Fϵ′′ (a) E1 (γ̂12 (X)) + 2hE1 (γ̂1 (X)) + h2 + R3n (a, h), 2
(6.6)
where R3n (a, h) is a remainder term given as n o R3n (a, h) = −E1 Fϵ (a) + γ̂1 (X) + h Fϵ′ (a) o2 1n ′′ γ̂1 (X) + h Fϵ (a) 1 |γ̂1 (X) + h| > a1n1 + 2 1 ′′ ′′ ∗ Fϵ a + θ γ̂1 (X) + h − Fϵ (a) + 2 n o2 E1 γ̂1 (X) + h 1 |γ̂1 (X) + h| ≤ a1n1 + E1 Fϵ a + γ̂1 (X) + h 1 |γ̂1 (X) + h| > a1n1 . for some θ∗ ∈ (0, 1) and a1n1 is as defined in (A.2)2 . Therefore, using Hölder inequality |R3n (a, h)| ≤ 2P |γ̂1 (X) + h| > a1n1 | D1 n o1/2 n 2 o1/2 + P |γ̂1 (X) + h| > a1n1 | D1 Fϵ′ (a) E1 γ̂1 (X) + h n 4 o1/2 ′′ + Fϵ (a) E1 γ̂1 (X) + h + sup |y|≤a1n1
(6.7)
n o2 Fϵ′′ (a + y) − Fϵ′′ (a) E1 γ̂1 (X) + h .
Note that, by Markov’s inequality and using assumption (A.1), we can get CEE1 |γ̂1 (X)|k 1 = o β2 . EP(|γ̂1 (X) − E1 (γ̂1 (X))| > a1n1 | D1 ) ≤ ak1n1 n1
(6.8)
Then, by choosing h = −E1 (γ̂1 (X)) + Gϵ (a)Var1 (γ̂1 (X)), where Gϵ (a) = −Fϵ′′ (a)/(2Fϵ′ (a)) and a = Fϵ−1 (µ2 ) and we obtain, F̃ −1 (µ2 ) = Fϵ−1 (µ2 ) − E1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ2 ) Var1 (γ̂1 (X)) + R4n , (6.9) 19
where using (6.7), (6.8) and assumptions (A.1), we get E|R4n | = o(1/nβ1 2 ). This shows that (6.5) follows from (6.9). Therefore, we have E(q̃u ) = EE(q̃u | D1 ) = Fϵ−1 (µ2 ) − EE1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ2 ) EVar1 (γ̂1 (X)) σ22 −1 ′′ 1 1 . + Fϵ (µ2 ) + o β2 + 2 n2 n1
(6.10)
4 Also, by using (A.1) and (A.2)4 , from (6.7) and (6.8), we get E|R4n | = o(n−β 1 ) and consequently we get
n o2 VarE(q̃u | D1 ) = VarE1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ2 ) VarVar1 (γ̂1 (X)) 1 1 −1 − 2Gϵ Fϵ (µ2 ) Cov E1 (γ̂1 (X)), Var1 (γ̂1 (X)) + o β4 + . n2 n1
(6.11)
Next, we calculate EVar(q̃u | D1 ). Using Taylor’s theorem, we have Var(q̃u | D1 ) = Var1 F̃ −1 (U(m2 ) ) −1 −1′ = Var1 F̃ (µ2 ) + (U(m2 ) − µ2 )F̃ (µ2 ) −1 −1′ − F̃ (µ2 ) + (U(m2 ) − µ2 )F̃ (µ2 ) 1(|U(m2 ) − µ2 | > a2n2 ) −1′ −1′ + (U(m2 ) − µ2 )(F̃ (µ∗ ) − F̃ (µ2 ))1(|U(m2 ) − µ2 | ≤ a2n2 )
(6.12)
+ Var1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | > a2n2 ) −1 − E1 F̃ (U(m2 ) )1(|U(m2 ) − µ2 | ≤ a2n2 ) E1 F̃ −1 (U(m2 ) )1(|U(m2 ) − µ2 | > a2n2 ) n o2 ′ = σ22 F̃ −1 (µ2 ) + R5n , where R5n is the remainder term such that, following similar calculations as in (6.4), under assumptions (A.2)4 and (A.3), E|R5n | = o(1/n2 ). Therefore, EVar(q̃u | D1 ) = EVar1 (F̃ −1 (µ2 )) n o2 1 2 −1′ = σ2 E F̃ (µ2 ) + o n2 n o2 1 1 ′ = σ22 Fϵ−1 (µ2 ) + o β4 + . n2 n1
20
(6.13)
Finally, we get n o2 VarVar1 (γ̂1 (X)) Var(q̃u ) = VarE1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ2 ) − 2Gϵ Fϵ−1 (µ2 ) Cov E1 (γ̂1 (X)), Var1 (γ̂1 (X)) n o2 ′ 1 1 + σ22 Fϵ−1 (µ2 ) + o β4 + . n2 n1
(6.14)
For the covariance, note that if α1 < (1 − α2 ), then we have Cov(q̃l , q̃u ) = Cov E(q̃l | D1 ), E(q̃u | D1 ) + ECov(q̃l , q̃u | D1 ), where q̃l and q̃u are the m1 = ⌊(n2 + 1)α1 ⌋ and m2 = ⌈(n2 + 1)(1 − α2 )⌉-th quantiles of {ϵ̃i }i∈D2 , as iid
th previously defined and U(m1 ) , U(m2 ) are the mth 1 and m2 order statistic of U1 , . . . , Un2 ∼ U nif (0, 1), with means µ1 = m1 /(n2 + 1) and µ2 = m2 /(n2 + 1), and covariance ρ12 = m1 (n2 + 1 − m2 )/((n2 + 1)2 (n2 + 2)). Then, following (6.3) and (6.9), we have Cov E(q̃l | D1 ), E(q̃u | D1 ) = VarE1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ1 ) Gϵ Fϵ−1 (µ2 ) VarVar1 (γ̂1 (X)) −1 −1 − Gϵ Fϵ (µ1 ) + Gϵ Fϵ (µ2 ) Cov E1 (γ̂1 (X)), Var1 (γ̂1 (X)) 1 1 + o β4 + . n2 n1
Furthermore, following similar calculations as in (6.12) and (6.13), it is straightforward to show that ECov(q̃l , q̃u | D1 ) = ECov1 (F̃ −1 (µ1 ), F̃ −1 (µ2 )) 1 −1′ −1′ = E F̃ (µ1 )F̃ (µ2 ) Cov1 (U(m1 ) , U(m2 ) ) + o n2 1 1 ′ ′ = ρ12 Fϵ−1 (µ1 )Fϵ−1 (µ2 ) + o β4 + . n 2 n1 Then we get, Cov(q̃l , q̃u ) = VarE1 (γ̂1 (X)) + Gϵ Fϵ−1 (µ1 ) Gϵ Fϵ−1 (µ2 ) VarVar1 (γ̂1 (X)) −1 −1 − Gϵ Fϵ (µ1 ) + Gϵ Fϵ (µ2 ) Cov E1 (γ̂1 (X)), Var1 (γ̂1 (X)) 1 1 −1 ′ −1 ′ . + ρ12 Fϵ (µ1 )Fϵ (µ2 ) + o β4 + n2 n1
(6.15)
Note that, L̃ = q̃u − q̃l , where q̃u ≡ q̃n2 ,1−αU and q̃l ≡ q̃n2 ,αL such that αL < (1 − αU ). Therefore, E(L̃) = E(q̃u ) − E(q̃l ), Var(L̃) = Var(q̃u ) + Var(q̃l ) − 2Cov(q̃u , q̃l ). Then the result follows from (6.10), (6.14) and (6.15) by considering α1 = αL and α2 = αU . 21
Proof of Theorem 3.1. Proof of this theorem closely follows the proof of Theorem 3.5. For α ∈ h (βϵ + 1)/(βϵ (n2 + 1)), 1 − 1/(βϵ (n2 + 1)) , we have E(q̂) = EE(q̂ | D1 ), Var(q̂) = VarE(q̂ | D1 ) + EVar(q̂ | D1 ). Let F̂ be the CDF of |ϵ̂i | = |ϵi − γ̂(Xi )|, i = n1 + 1, . . . , n, given D1 . Assume U1 , . . . , Un2 be a set of uniformly distributed random variables on (0, 1). Define, U(m) as the m = ⌈(n2 + 1)(1 − α)⌉-th iid
2 order statistic of U1 , . . . , Un2 ∼ U nif p (0, 1), with mean µα and variance σα . Recall that, we have a2n2 = C log(n2 )/n2 , for some constant C > 0. Then, following similar calculations as in (6.3), and (6.4), we have
E(q̂ | D1 ) = E1 F̂ −1 (U(m) ) = F̂ −1 (µα ) + 0 +
σα2 −1′′ F̂ (µα ) + R6n , 2
(6.16)
where
′
R6n = −E1 F̂ −1 (µα ) + (U(m) − µα )F̂ −1 (µα ) 1 ′′ + (U(m) − µα )2 F̂ −1 (µα ) 1(|U(m) − µα | > a2n2 ) 2 1 ′′ ′′ + E1 (U(m) − µα )2 F̂ −1 (µ∗ ) − F̂ −1 (µα ) 1(|U(m) − µα | ≤ a2n2 ) 2 + E1 F̂ −1 (U(m) )1(|U(m) − µα | > a2n2 ), where µ∗ is some point between U(m) and µα and we have by assumptions (A.2)1 and (A.3), E|R6n | = o(1/n2 ). Next, we show that the quantiles (conditional on D1 ) of |ϵ̂i | are close to the quantiles of |ϵi |. Note that for some fixed a ad h, by Taylor’s theorem F̂ (a + h) = F̃ (a + h) − F̃ (−a − h) ′ = F|ϵ| (a) + hF|ϵ| (a) + E1 (γ̂1 (X)){Fϵ′ (a) − Fϵ′ (−a)} + R7n (a, h),
where F̃ is the CDF of ϵ̂i conditional on D1 and R7n (a, h) is a remainder term which can be bounded by following similar derivations in (6.7), (6.8) and using assumption (A.2)1 . Therefore, by choosing ′ (a) and a = F −1 (µ ), we obtain, h = H|ϵ| (a)E1 (γ̂1 (X)), where H|ϵ| (a) = −{Fϵ′ (a) − Fϵ′ (−a)}/F|ϵ| α |ϵ| −1 −1 F̂ −1 (µα ) = F|ϵ| (µα ) + H|ϵ| F|ϵ| (µα ) E1 (γ̂1 (X)) + R9n , 1 where E|R9n | = o(n−β 1 ). Therefore, we have
E(q̂) = EE(q̂ | D1 ) −1 = F|ϵ| (µα ) + H|ϵ|
−1 F|ϵ| (µα )
σα2 −1 ′′ 1 1 EE1 (γ̂1 (X)) + F (µα ) + o β1 + . 2 |ϵ| n2 n1
and by similar argument as in (6.11), (6.12) and (6.13), we get n o2 n o2 ′ 1 1 −1 −1 2 Var(q̂) = H|ϵ| F|ϵ| (µα ) VarE1 (γ̂1 (X)) + σα F|ϵ| (µα ) + o β2 + . n2 n1 22
Note that, if the distribution of ϵ is symmetric, then we would have H|ϵ| (a) = 0, for all a and the new correction factor h would be h = Gϵ (a)E1 (γ̂12 (X)) where Gϵ (a) = −Fϵ′′ (a)/(2Fϵ′ (a)), and we would get σ 2 −1 ′′ EE1 (γ̂12 (X)) + α F|ϵ| (µα ) + o
1 . + 2 n2 nβ1 2 n o2 n o2 ′ 1 1 −1 −1 2 2 VarE1 (γ̂1 (X)) + σα F|ϵ| (µα ) + o β4 + Var(q̂) = Gϵ F|ϵ| (µα ) . n2 n1 −1 E(q̂) = F|ϵ| (µα ) + Gϵ
−1 F|ϵ| (µα )
1
Finally, it follows that E(L̂) = 2E(q̂) and Var(L̂) = 4Var(q̂).
References [1] Concrete compressive strength data set. http://archive.ics.uci.edu/ml/datasets/ concrete+compressive+strength. Accessed: January, 2019. [2] Anastasios N Angelopoulos and Stephen Bates. Conformal prediction: A gentle introduction. Foundations and Trends in Machine Learning, 16(4):494–591, 2023. [3] Anastasios Nikolas Angelopoulos, Stephen Bates, Michael Jordan, and Jitendra Malik. Uncertainty sets for image classifiers using conformal prediction. In International Conference on Learning Representations, 2021. [4] Anastasios Nikolas Angelopoulos, Rina Barber, and Stephen Bates. Online conformal prediction with decaying step sizes. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 1616–1630, 21–27 Jul 2024. [5] Yu Bai, Song Mei, Huan Wang, Yingbo Zhou, and Caiming Xiong. Efficient and differentiable conformal prediction with general function classes. In International Conference on Learning Representations, 2022. [6] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. Predictive inference with the jackknife+. The Annals of Statistics, 49(1):486–507, 2021. [7] Peter J Bickel. Some contributions to the theory of order statistics. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics, volume 5, pages 575–592. University of California Press, 1967. [8] Margarida Campos, António Farinhas, Chrysoula Zerva, Mário AT Figueiredo, and André FT Martins. Conformal prediction for natural language processing: A survey. Transactions of the Association for Computational Linguistics, 12:1497–1516, 2024. [9] John J Cherian, Isaac Gibbs, and Emmanuel J Candès. Large language model validity via enhanced conformal prediction methods. Advances in Neural Information Processing Systems, 37:114812–114842, 2024. [10] Thibault Cordier, Vincent Blot, Louis Lacombe, Thomas Morzadec, Arnaud Capitaine, and Nicolas Brunel. Flexible and systematic uncertainty estimation with conformal prediction via the mapie library. In Conformal and Probabilistic Prediction with Applications, pages 549–581. PMLR, 2023. 23
[11] Nicolas Deutschmann, Mattia Rigotti, and Maria Rodriguez Martinez. Adaptive conformal regression with split-jackknife+ scores. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. [12] Guneet S Dhillon, George Deligiannidis, and Tom Rainforth. On the expected size of conformal prediction sets. In International Conference on Artificial Intelligence and Statistics, pages 1549–1557. PMLR, 2024. [13] Tiffany Ding, Anastasios Angelopoulos, Stephen Bates, Michael Jordan, and Ryan J Tibshirani. Class-conditional conformal prediction with many classes. Advances in neural information processing systems, 36:64555–64576, 2023. [14] Jianqing Fan. Local polynomial modelling and its applications: monographs on statistics and applied probability 66. Routledge, 2018. [15] Shai Feldman, Stephen Bates, and Yaniv Romano. Improving conditional coverage via orthogonal quantile regression. Advances in neural information processing systems, 34:2060–2071, 2021. [16] Feichen Gan, Lu Youcun, Yingying Zhang, and Yukun Liu. Conformal prediction beyond the horizon: Distribution-free inference for policy evaluation. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. [17] Isaac Gibbs and Emmanuel J Candès. Conformal inference for online prediction with arbitrary distribution shifts. Journal of Machine Learning Research, 25(162):1–36, 2024. [18] László Györfi, Michael Kohler, Adam Krzyżak, and Harro Walk. A distribution-free theory of nonparametric regression. Springer, 2002. [19] Peter Hall, Joel L Horowitz, and Bing-Yi Jing. On blocking rules for the bootstrap with dependent data. Biometrika, 82(3):561–574, 1995. [20] Shayan Kiyani, George J Pappas, and Hamed Hassani. Length optimization in conformal prediction. Advances in Neural Information Processing Systems, 37:99519–99563, 2024. [21] Michael Kohler and Sophie Langer. On the rate of convergence of fully connected deep neural network regression estimates. The Annals of Statistics, 49(4):pp. 2231–2249, 2021. [22] Jing Lei and Larry Wasserman. Distribution-free prediction bands for non-parametric regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 76(1):71–96, 2014. [23] Jing Lei, James Robins, and Larry Wasserman. Efficient nonparametric conformal prediction regions. arXiv preprint arXiv:1111.1418, 2011. [24] Jing Lei, James Robins, and Larry Wasserman. Distribution-free prediction sets. Journal of the American Statistical Association, 108(501):278–287, 2013. [25] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman. Distribution-free predictive inference for regression. Journal of the American Statistical Association, 113(523):1094–1111, 2018.
24
[26] Charles Lu, Andréanne Lemay, Ken Chang, Katharina Höbel, and Jayashree Kalpathy-Cramer. Fair conformal predictors for applications in medical imaging. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 12008–12016, 2022. [27] Rachel Luo, Shengjia Zhao, Jonathan Kuck, Boris Ivanovic, Silvio Savarese, Edward Schmerling, and Marco Pavone. Sample-efficient safety assurances using conformal prediction. The International Journal of Robotics Research, 43(9):1409–1424, 2024. [28] Paul Melki, Lionel Bombrun, Boubacar Diallo, Jérôme Dias, and Jean-Pierre Da Costa. Groupconditional conformal prediction via quantile regression calibration for crop and weed classification. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 614–623, 2023. [29] Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In European conference on machine learning, pages 345–356. Springer, 2002. [30] Olivier C Pasche, Henry Lam, and Sebastian Engelke. Extreme conformal prediction: Reliable intervals for high-impact events. arXiv preprint arXiv:2505.08578, 2025. [31] Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019. [32] Yaniv Romano, Matteo Sesia, and Emmanuel Candes. Classification with valid and adaptive coverage. Advances in neural information processing systems, 33:3581–3591, 2020. [33] Mauricio Sadinle, Jing Lei, and Larry Wasserman. Least ambiguous set-valued classifiers with bounded error levels. Journal of the American Statistical Association, 114(525):223–234, 2019. [34] Robert J Serfling. Approximation theorems of mathematical statistics. John Wiley & Sons, 2009. [35] Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008. [36] Akshai Parakkal Sreenivasan, Aina Vaivade, Yassine Noui, Payam Emami Khoonsari, Joachim Burman, Ola Spjuth, and Kim Kultima. Conformal prediction enables disease course prediction and allows individualized diagnostic uncertainty in multiple sclerosis. npj Digital Medicine, 8 (1):224, 2025. [37] Stephen Mack Stigler. Linear Functions of Order Statistics. The Annals of Mathematical Statistics, 40(3):770 – 788, 1969. [38] David Stutz, Krishnamurthy Dj Dvijotham, Ali Taylan Cemgil, and Arnaud Doucet. Learning optimal conformal classifiers. In International Conference on Learning Representations, 2022. [39] Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas. Conformal prediction under covariate shift. Advances in neural information processing systems, 32, 2019. [40] Aad W Van der Vaart. Asymptotic statistics, volume 3. Cambridge university press, 2000. [41] Vladimir Vovk. Cross-conformal predictors. Annals of Mathematics and Artificial Intelligence, 74:9–28, 2015. 25
[42] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. Springer, 2005. [43] Vladimir Vovk, Valentina Fedorova, Ilia Nouretdinov, and Alexander Gammerman. Criteria of efficiency for conformal prediction. In Symposium on conformal and probabilistic prediction with applications, pages 23–39. Springer, 2016. [44] Ran Xie, Rina Barber, and Emmanuel Candes. Boosted conformal prediction intervals. Advances in Neural Information Processing Systems, 37:71868–71899, 2024. [45] Chen Xu and Yao Xie. Conformal prediction for time series. IEEE transactions on pattern analysis and machine intelligence, 45(10):11575–11587, 2023. [46] Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Dieuleveut. Adaptive conformal predictions for time series. In International Conference on Machine Learning, pages 25834–25866. PMLR, 2022.
26