Keywords: Conformal prediction; optimal sample allocation; asymptotic optimality; distribution-free 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 [1]. 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: it partitions data into (proper) training- and calibration-subsets and achieves valid coverage with a single model fit [3]. Other extensions, such as cross-conformal prediction [4] and Jackknife+ [5], 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 [6], [7].

Beyond computational efficiency, recent research has substantially expanded the applicability of conformal prediction. Adaptive techniques [8], blockwise calibration methods [9], and covariate-adjusted conformal inference [10] have been developed to address challenges such as non-exchangeable data and distributional shift. Further efforts focus on achieving stronger forms of validity, including class-conditional coverage guarantees [11]. 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 [13], natural language processing [14], computer vision [15], [16], finance [17], medical imaging [18], robotics and autonomous vehicles [19], reinforcement learning [20], and large language models [21], 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]–[25]. 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 [26]–[31]. Another direction leverages insights from level-set estimation, drawing on foundations from classical statistics and their adaptation to conformal prediction in [23], [32], to construct covariate-dependent thresholds that directly target shorter prediction sets while accommodating refined coverage notions. Complementing these developments, conformal training methods [21], [33], [34] 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 [35].

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:

• 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 [1], [36]. Let \((X_i,Y_i)_{i=1}^{n+1}\) be a random sample from an unknown distribution \(\mathbf{P}_{XY},\) where \(X_i\in\mathcal{X}\) and \(Y_i\in\mathcal{Y}\). Based on the first \(n\) observations, the objective is to construct a measurable set \(\mathcal{C}(X_{n+1})\subseteq\mathcal{Y}\) such that, for a prespecified miscoverage level \(\alpha \in (0,1),\) \[\Pr\{Y_{n+1}\in \mathcal{C}(X_{n+1})\} \ge 1-\alpha,\] without imposing parametric assumptions on \(\mathbf{P}_{XY}\). In this paper, we focus on the regression setting, where \(Y\in\mathbb{R}\) and \[Y = f(X) + \epsilon,\qquad X\sim \mathbf{P}_X,\qquad \epsilon \overset{\mathrm{iid}}{\sim} F_\epsilon,\] so that \(\mathcal{C}(X_{n+1})\) is an interval.

2.1 Split conformal prediction↩︎

Split conformal prediction divides the observed sample \((X_i,Y_i)_{i=1}^n\) into a proper training set, say \(\mathcal{D}_1\), of size \(n_1\) and a calibration set, say \(\mathcal{D}_2\), of size \(n_2=n-n_1\). A regression estimator \(\hat{f}_1 \equiv \hat{f}_{n_1}\) is fitted only on the proper training data \(\mathcal{D}_1\), 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 [3], [26]. The choice of \((n_1,n_2)\) determines the trade-off between estimator accuracy (improved by a larger \(n_1\)) and coverage accuracy (improved by a larger \(n_2\)). 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]: \[\hat{\epsilon}_i = |Y_i - \hat{f}_1(X_i)|, \qquad i \in \mathcal{D}_2.\] Let \(\hat{q} \equiv \hat{q}_{n_2,1-\alpha}\) denote the \(\lceil (n_2+1)(1-\alpha)\rceil\)-th order statistic of \(\{\hat{\epsilon}_i : i\in \mathcal{D}_2\}\) where, for \(x\in \bbr\), \(\lceil x\rceil\) (and \(\lfloor x \rfloor\)) 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 \[\hat{\mathcal{C}}(X_{n+1}) = \big[\hat{f}_1(X_{n+1}) - \hat{q},\; \hat{f}_1(X_{n+1}) + \hat{q} \big].\] Exchangeability yields the finite-sample guarantee \[1-\alpha \leq \Pr\Big(Y_{n+1}\in \hat{\mathcal{C}}(X_{n+1})\Big) \leq 1-\alpha + \frac{1}{n_2+1}.\] In symmetric conformal prediction, the interval length is \[\mathcal{\hat{L}}= 2 \hat{q}\] and its asymptotic behavior depends on different factors like the distribution of \(\epsilon,\) and the choice of \((n_1,n_2)\), 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 [37]. Let \(\alpha_L,\alpha_U\ge 0\) satisfy \(\alpha_L+\alpha_U=\alpha \in (0,1)\). Define the signed conformity scores \[\tilde{\epsilon}_i = Y_i - \hat{f}_1(X_i) , \qquad i \in \mathcal{D}_2.\] Let \(\tilde{q}_l \equiv \tilde{q}_{n_2,\alpha_L}\) and \(\tilde{q}_u \equiv \tilde{q}_{n_2,1-\alpha_U}\) denote the \(\lfloor (n_2+1)\alpha_L\rfloor\)-quantile of \(\{\tilde{\epsilon}_i : i\in \mathcal{D}_2\}\) and the \(\lceil (n_2+1)(1-\alpha_U)\rceil\)-quantile of \(\{\tilde{\epsilon}_i : i\in \mathcal{D}_2\}\) respectively. Then the asymmetric interval is \[\tilde{\mathcal{C}}(X_{n+1}) = \big[\hat{f}_1(X_{n+1}) + \tilde{q}_l,\; \hat{f}_1(X_{n+1}) + \tilde{q}_u \big].\] Exchangeability again ensures that \[1-\alpha \leq \Pr\Big(Y_{n+1}\in \tilde{\mathcal{C}}(X_{n+1})\Big) \leq 1-\alpha + \frac{2}{n_2+1} ,\] and the interval length is \[\mathcal{\tilde{L}} = \tilde{q}_u - \tilde{q}_l.\] Later we will see that in many cases, the use of asymmetric intervals can be advantageous over the symmetric ones.

Under the model \(Y = f(X) + \epsilon,\) with \(X \sim \mathbf{P}_X\) and \(\epsilon\overset{\mathrm{iid}}{\sim} F_\epsilon,\) both symmetric and asymmetric split conformal intervals provide exact finite-sample coverage bounds, which primarily depend on the calibration sample size \(n_2\) through the general upper bound \(\frac{k}{n_2+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 \(n_1\) (through the estimator \(\hat{f}_{1}\)) as well as the calibration set size \(n_2\) (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(\mathcal{\hat{L}})\) (or \(\E(\mathcal{\tilde{L}})\)) and \(\Var(\mathcal{\hat{L}})\) (or \(\Var(\mathcal{\tilde{L}})\)) for general \((n_1,n_2)\), 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, \(C_1, C_2,\dots\) denote constants that are independent of the sample size \(n\). For two positive sequences \(a_n\) and \(b_n,\) we write \(a_n=O(b_n)\) if there exists a constant \(C\) such that \(a_n \leq Cb_n\) for all \(n\). Similarly, \(a_n = o(b_n)\) indicates that \(a_n/b_n \to 0,\) as \(n\to\infty\). Furthermore, by \(a_n \asymp b_n\), we denote that \(a_n=O(b_n)\) and \(b_n=O(a_n),\) simultaneously. Recall that the floor and ceiling functions, denoted by \(\lfloor x \rfloor\) and \(\lceil x \rceil,\) represent the greatest integer less than or equal to \(x,\) and the smallest integer greater than or equal to \(x,\) respectively.

For an estimator \(\hat{f}_n\) of \(f\) based on sample size \(n\), define the estimation error at a point \(x\) by \(\hat{\gamma}_n(x) = \hat{f}_n(x) - f(x),\) which measures the discrepancy between the estimator \(\hat{f}_n(x)\) and the true regression function \(f(x)\). For notational conveniences, we simply write the powers \((\hat{\gamma}_n(x))^k\) as \(\hat{\gamma}_n^k(x)\) for any \(k\geq2\). We denote by \(\E_1,\) \(\Var_1,\) and \(\Cov_1\) the conditional expectation, variance, and covariance, respectively, given the proper training set \(\mathcal{D}_1\). For example, the \(L_2\) error can be written as \[\E\int |\hat{f}_n(x) - f(x)|^2 \mathbf{P}_X(dx)=\E\E_1(\hat{\ga}_n^2(X)).\] Similarly, the conditional variance is written as \[\E\Var_1(\hat{\ga}_n(X)) = \E\E_1(\hat{\ga}_n^2(X)) - \E\{\E_1(\hat{\ga}_n(X))\}^2.\] Furthermore, for notational brevity, we will write \(\hat{\ga}_{n_1} = \hat{f}_{n_1} - f\) as \(\hat{\gamma}_1\).

3.2 Conditions↩︎

For the symmetric split conformal prediction interval, define \[m_\alpha=\left\lceil (n_2+1)(1-\alpha)\right\rceil,\qquad \mu_\alpha=\frac{m_\alpha}{n_2+1},\qquad \sigma_\alpha^2=\frac{\mu_\alpha(1-\mu_\alpha)}{n_2+2}.\]

For the asymmetric split conformal prediction interval, define \[m_1=\left\lfloor (n_2+1)\alpha_L\right\rfloor,\qquad m_2=\left\lceil (n_2+1)(1-\alpha_U)\right\rceil,\] and \[\mu_i=\frac{m_i}{n_2+1},\qquad \sigma_i^2=\frac{\mu_i(1-\mu_i)}{n_2+2}, \qquad i=1,2,\] together with \[\rho_{12}=\frac{\mu_1(1-\mu_2)}{n_2+2}.\]

Let \(S\) denote the population conformity score, with distribution function \(F_S\). Throughout the paper, \[S= \begin{cases} |\epsilon|, & \text{for symmetric split conformal prediction},\\ \epsilon, & \text{for asymmetric split conformal prediction}. \end{cases}\] Furthermore, let \(\delta\in(0,1]\) be a constant independent of \(n\). We make the following assumptions.

1. For some integer \(K\geq 1,\) suppose there exists a positive, strictly increasing sequence of constants \(\{\beta_k\}_{k=1}^K,\) independent of \(n,\) such that \[\E\int|\hat{f}_{n}(x) - f(x)|^k \;\mathbf{P}_X(dx) \asymp n^{-\beta_k}.\]

2. \(_r\) Suppose there exists two sequences defined by \(a_{1n} = \delta^{-1}n^{-\gamma},\) for some \(\gamma\in(0,(\beta_k-\beta_r)/k),k>r\) and \(a_{2n} = \delta^{-1}\sqrt{\log(n)/n},\) such that for each relevant quantile level \(\mu\), \[\delta<\inf_{|y|\leq a_{1n},|x|\leq a_{2n}} F_S'\Big(F_S^{-1} (\mu+x)+y\Big) \leq \sup_{|y|\leq a_{1n},|x|\leq a_{2n}} F_S'\Big(F_S^{-1} (\mu+x)+y \Big) <\delta^{-1},\] \[\sup_{|y|\leq a_{1n},|x|\leq a_{2n}} \Big|F_S''\Big(F_S^{-1}(\mu+x)+y\Big)\Big|<\delta^{-1},\] and \[\sup_{|y|\leq a_{1n},|x|\leq a_{2n}} \Big|F_S''\Big(F_S^{-1}(\mu+x)+y\Big)-F_S''\Big(F_S^{-1}(\mu)\Big) \Big| = o(1), \text{ as } n\to\infty.\] For Theorem 1, corresponding to the symmetric interval, this condition is imposed with \(F_S=F_{|\epsilon|}\) and \(\mu=\mu_\alpha\). For Theorem 5, corresponding to the asymmetric interval, it is imposed with \(F_S=F_\epsilon\) and \(\mu\in\{\mu_1,\mu_2\}\).

3. There exists a constant \(\beta_\epsilon>0,\) independent of \(n,\) such that for all \(x\geq \delta^{-1},\) \[F_\epsilon(-x) + (1-F_\epsilon(x)) \leq \delta^{-1} |x|^{-\beta_\epsilon}.\]

Assumption [A3] imposes polynomial moment bounds on the estimation error of the regression estimator \(\hat{f}_n\). Specifically, it requires that for some integer \(K \ge 1,\) the \(L_k\)-risk \(\E\!\int |\hat{f}_n(x)-f(x)|^k\,\mathbf{P}_X(dx)\) decays at the rate \(n^{-\beta_k}\) for \(k=1,\dots,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, [38], [39]. 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 [26].

Assumption [A1]\(_r\) imposes smoothness and local stability conditions on the noise distribution function \(F_S\) (where \(S\) is either \(|\epsilon|\) or \(\epsilon\)). Specifically, it requires that the first and second derivatives of \(F_S\) 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 [40]. The final requirement in Assumption [A1]\(_r,\) which controls the difference between second derivatives evaluated at nearby points, is essentially a local Lipschitz condition on \(F_S''\). Such derivative-boundedness and local smoothness conditions are standard in the analysis of quantile estimators and residual-based inference; see, for example, [41], where smoothness of the noise density plays a key role in obtaining uniform expansions of sample quantiles.

Assumption [A2] imposes a polynomial tail condition on the noise distribution \(F_\epsilon\). 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, [42] and [43]. In the present setting, Assumption [A2] 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 [A1]\(_r\).

Taken together, assumptions [A3]–[A2] 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 [A3]–[A2], 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.

The next theorem gives the mean and variance of the length of the asymmetric split conformal prediction interval.

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 \(\sigma_\alpha^2 \asymp n_2^{-1}\) always holds, it suffices to characterize \(\E\E_1(\hat{\gamma}_1^2(X))\) and \(\Var\E_1(\hat{\gamma}_1^2(X)),\) equivalently the expectation and variance of the \(L_2\) errors, as functions of the (proper) training size \(n_1\) for each model.

3.4.1 Linear regression↩︎

We begin by considering a random effects model with the linear regression function \(f(x) = x^\top\beta,\) where \(x\in\mathbb{R}^d\) and \(\beta\in\mathbb{R}^d\) are the covariates and regression coefficients, respectively. For illustrative purposes, let \(\{(X_i,Y_i)\}_{i=1}^{n_1}\) denote a collection of \(n_1\) independent and identically distributed samples. We consider a Gaussian design with \(X_i\overset{iid}{\sim} \mathcal{N}(0, I_d),\) where \(I_d\) is the \(d\)-dimensional identity matrix, and errors \(\epsilon_i \overset{iid}{\sim} \mathcal{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 \(\hat{f}(x) = x^\top\hat{\beta},\) where \(\hat{\beta} = (X^\top X)^{-1}X^\top Y\) and we have \[\E\E_1(\hat{\gamma}_1^2(X)) = \frac{d}{n_1-d-1} \asymp \frac{1}{n_1},\] and \[\Var\E_1(\hat{\gamma}_1^2(X)) = 2\E\Big[\text{tr}(X^\top X)^{-2}\Big] + \Var\Big[\text{tr}(X^\top X)^{-1}\Big] \asymp \frac{1}{n_1^2},\] provided \(n_1 > d+3\). Then, in ordinary least squares linear regression, a bias-optimal split requires \(n_1 \asymp n_2\), that is, the training and calibration sample sizes are of the same order. For an MSE-optimal split, one requires \(n_1 \asymp \sqrt{n_2}\).

3.4.2 Nonparametric regression↩︎

We consider a nonparametric regression framework in which the goal is to estimate the regression function \(f(x)=\E[Y \mid X=x]\) from \(n_1\) iid samples \(\{(X_i, Y_i)\}_{i=1}^{n_1}\). We assume that the covariates \(X_i\in[-a, a]^d,\) for some \(a>0\) are uniformly distributed, and the noise terms \(\{\epsilon_i\}_{i=1}^{n_1}\) are iid with \(\epsilon_i\overset{iid}{\sim}\mathcal{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 \(\mathcal{K}:\mathbb{R}^d \to \mathbb{R}_+\) that satisfies \(\int\mathcal{K}(u)\,du=1\). The estimator at a point \(x\in\mathbb{R}^d\) is defined as \[\hat{f}_{h,n_1}(x) = \frac{\sum_{i=1}^{n_1} \mathcal{K}\Big(\frac{x - X_i}{h}\Big) Y_i}{\sum_{i=1}^{n_1} \mathcal{K}\Big(\frac{x - X_i}{h}\Big)},\] which can be viewed as a locally weighted average of the observed responses \(Y_i,\) with larger weights assigned to observations where the covariates \(X_i\) are closer to \(x\). For illustrative purposes, we assume that the kernel \(\mathcal{K}\) is bounded, symmetric, and continuously differentiable with compact support. Finally, if we choose the bandwidth \(h\asymp n_1^{-1/(4+d)},\) standard calculation shows \[\E\E_1(\hat{\gamma}_1^2(X)) \asymp \frac{1}{n_1^{4/(4+d)}},\] and \[\Var\E_1(\hat{\gamma}_1^2(X)) \asymp \frac{1}{n_1^{8/(4+d)}},\] provided the regression function has bounded first and second derivatives. Then, in the Nadaraya-Watson kernel regression, the theoretically justified choice for the bias-optimal data splitting is \(n_1\asymp n_2^{1+d/4},\) implying that the training sample should be considerably larger than the calibration sample and \(n_1\asymp n_2^{1/2+d/8}\) 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 \(L_2\)-error are not yet available in the literature. However, bounds on the expected \(L_2\)-error have been derived in [44], under certain structural assumptions on the regression function. To state these results, we rely on the following two definitions.

With the definition of \((p,C)\)-smoothness, we define the class of hierarchical composition models as follows.

Let the corresponding regression function \(f\) belong to the class \(\mathcal{H}(l,p)\) for some \(l\in\mathbb{N}\) and \(\mathcal{P}\subseteq (0,\infty)\times\mathbb{N}\). Then according to Theorem 1 in [44], and under some regularity conditions, we have \[\E\E_1(\hat{\gamma}_1^2(X)) = O\Big( (\log(n))^6 \max_{(p,K)\in \mathcal{P}} n^{-\frac{2p}{2p+K}} \Big).\] Although the exact order is not available, for a bias-optimal split, we can use the approximate relation \(n_1 \asymp n_2^{1+\bar{K}/2\bar{p}},\) where \((\bar{p},\bar{K})\in \mathcal{P}\) is chosen such that \((\bar{p},\bar{K})=\arg\min_{(p,K)\in \mathcal{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 \(\mathcal{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. [45]).

Based on the preceding results, we obtain that the optimal split satisfies \(n_1 \asymp n_2^b\) for some known constant \(b>0\). Let \(\hat{c}_\ell\) denote the empirical estimate of the ratio of the training sample size to the total sample size, computed from a subsample of size \(\ell (<n)\). The corresponding estimate \(\hat{c}_n\) for the full sample size \(n\) is then determined as the solution to the equation \[\label{eq:ratio} x = \frac{n^{b-1}\hat{c}_\ell}{\ell^{b-1}(1-\hat{c}_\ell)^b} (1-x)^b.\tag{1}\]

The following algorithm provides a way to approximate the ratio \(\hat{c}_{\ell}\) where \(\ell<n\).

Algorithm Description. Given \(\alpha \in [1/(n_2+1),1),\) and \(\ell(<n),\) we choose \(\hat{c}_\ell\) as follows:

1. Divide the dataset randomly into \(K\) disjoint subsets \(\{\mathcal{Q}_{k}\}_{k=1}^K,\) each containing \(\ell\) samples, such that \(K\ell\leq n\).

2. Choose a grid \(\mathcal{G}_g = \{ p_1, \ldots, p_g \} \in (0,1)^g\) of size \(g\) (e.g., \(\mathcal{G}_6 = \{0.4, 0.5, \ldots, 0.9\}\)).

3. For each \(k = 1, 2, \ldots, K\):

   1. For each candidate proportion \(p_i \in \mathcal{G}_g,\) \(i=1,\dots,g,\) further split the subsets \(\mathcal{Q}_{k}\) into a proper training set \(\mathcal{T}_{k,i}\) of size \(\lfloor \ell p_i \rfloor\) and a calibration set \(\mathcal{C}_{k,i}\) of size \(\ell - \lfloor \ell p_i \rfloor\).

   2. Using \(\mathcal{T}_{k,i}\) and \(\mathcal{C}_{k,i},\) compute prediction interval length \(\hat{\mathcal{L}}(\mathcal{T}_{k,i},\mathcal{C}_{k,i},p_i)\) and select the target proportion \(\hat{p}_{k}= \arg\min _{p_i \in \mathcal{G}_g}\hat{\mathcal{L}}(\mathcal{T}_{k,i},\mathcal{C}_{k,i},p_i).\)

4. Set \(\hat{c}_\ell = \frac{1}{K}\sum_{k=1}^{K}\hat{p}_k\)

5. (Optional refinement) If greater precision is required, repeat the Steps 2–5 using a finer grid centered around the current estimate \(\hat{c}_\ell,\) and update \(\hat{c}_\ell\) accordingly.

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-\alpha=0.9\) and report the average prediction interval length.

4.1 Linear regression↩︎

We examine the following linear regression model: \[\begin{align} \label{eq:linear95regression} y_i = x_i^\top\beta + \epsilon_i, \qquad i=1,\ldots,n, \end{align}\tag{2}\] where \(y_i\in\mathbb{R}\) is the response variable, \(x_i\in[0,1]^5\) denotes the regressors, \(\beta\in\mathbb{R}^5\) represents the regression coefficients, and \(\epsilon_i\in\mathbb{R}\) is the error term. Using equation 2 , we generate synthetic datasets with varying sample sizes \(\{100,300,500,800,2000\}\). The regressors \(x_i\) are sampled independently from a uniform distribution on \([0,1]^5\). We consider three types of error distributions:

• Normal distribution: \[\begin{align} \epsilon_i \sim \mathcal{N}(0, 1) \end{align}\]

• Student’s t-distribution with 5 degrees of freedom (variance-normalized): \[\begin{align} \epsilon_i \sim \frac{t(5)}{\sqrt{5/3}} \end{align}\]

• Lognormal distribution (standardized to zero mean and unit variance): \[\begin{align} \epsilon_i \sim \frac{\text{Lognormal}(0, 1) - \exp(0.5)}{\sqrt{(\exp(1) - 1)\exp(1)}} \end{align}\]

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 \(n_1 \asymp n_2\). That is, the training set size should be approximately equal to the calibration set size. Equivalently, this corresponds to setting \(b=1\) in Equation 1 . Solving the equation for \(b=1\) yields \(x=\hat{c}_{\ell}\), 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 \[\begin{align} y_{i,j} = f_j(x_i) + \lambda_j\epsilon_i, \quad i = 1, \ldots, n, \quad j=1,2, \end{align}\] where \(y_{i,j}\in\mathbb{R}\) is the response variable, \(f_j\) is the regression function, \(x_i=(x_i^{(1)},\ldots,x_i^{(5)})^\top\in\mathbb{R}^5\) is a five-dimensional predictor which is uniformly distributed over the region \([0,1]^5\) independently, \(\epsilon_i\) is the independent error term, which is also independent of \(x_i\) as defined in Subsection 4.1, and \(\lambda_1,\lambda_2\ge 0\) are fixed scaling parameters. For reasons of comparability, the scaling parameters \(\lambda_j\) are chosen to match \(20\%\) of the typical range of variation of \(f_j(x)\) under the distribution of \(x\). This range is chosen via the interquartile range (IQR) of \(f_j(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 \(\lambda_1=0.340\) and \(\lambda_2=0.262\). The regression functions under consideration are:

• Function \(f_1\): \[\begin{align} f_1(x) = 3x^{(1)}+\tan (x^{(2)})+(x^{(3)})^3+\log(x^{(4)}+0.1)+\sqrt{x^{(5)}+0.1} \end{align}\]

• Function \(f_2\): \[\begin{align} f_2(x) = \exp(\|x\|), \end{align}\] where \(x=(x^{(1)},\ldots,x^{(5)})^\top\).

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 satisfies \(n_1 \asymp n_2^{1+d/4}\), where \(d\) denotes the covariate dimension. In our experiments, \(d=5\), which yields \(n_1 \asymp n_2^{9/4}\). 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 1 is \(b=9/4\). Substituting this value into Equation 1 and considering the setting \(n=2000\) and \(\ell=300\), together with the empirically observed minimum prediction interval length at \(\hat{c}_{\ell}\approx 0.75\), yields \(x\approx 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 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: \[\begin{align} y_{i,j} = f_j(x_i) + \lambda_j\epsilon_i, \qquad i = 1, \ldots, n, \quad j=1,2. \end{align}\] 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 \(n_1 \asymp n_2^{1+\bar{K}/2\bar{p}}\). In our setting, \(\bar{p}=1\) and \(\bar{K}=5\) for the functions \(f_1\) and \(f_2\), which gives \(n_1 \asymp n_2^{7/2}\). As introduced in Section 3.4.3, \((\bar{p},\bar{K})\) denote the smoothness and order parameters in \(\mathcal{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 1 is \(b=7/2\). Substituting this value into Equation 1 with \(n=2000\) and \(\ell=100\), and using the empirically observed minimum prediction interval length at \(\hat{c}_{\ell}\approx 0.6\), yields \(x\approx 0.93\) (similar results are obtained if we consider higher values of \(\ell\)). 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 [46]. 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 1 and the proposed procedure for selecting a split ratio that minimizes prediction interval length, we set \(\bar{p}=1\) and \(\bar{K}=8\) (see Section 3.4.3). Although these parameters can be chosen in a data-dependent way, we only choose them heuristically. Under these choices, \(n_1 \asymp n_2^5\), corresponding to \(b=5\). To implement Equation 1 , the training data are partitioned into four subsets, giving \(\ell=206\) in Equation 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 \(\hat{c}_{\ell}=0.66\) is obtained. Substituting these values into Equation 1 yields \(x \approx 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.

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 train-calibration partition. This highlights the practical relevance and novelty of our approach.

6 Proof of the main results↩︎

First, we present the proof of Theorem 5 and then move on to prove Theorem 1.

References↩︎

---

1. Department of Statistics and Data Science, Washington University in St. Louis, St. Louis, MO 63130, USA↩︎

2. Department of Electrical and Systems Engineering, Washington University in St. Louis, St. Louis, MO 63130, USA↩︎