Paper deep dive
PICPIs: Prediction-Interval-Conditional Prediction Intervals
Xuelin Yang, Baihe Huang, Yilong Hou, Guido Imbens, Michael I. Jordan
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 9/23/2026, 4:26:52 AM
Summary
The paper introduces Prediction-Interval-Conditional Prediction Intervals (PICPIs), a framework for uncertainty quantification that conditions on prediction values rather than covariates. PICPIs satisfy a self-consistency condition where the mean outcome within a prediction interval lies within that same interval. The authors provide algorithms for constructing these intervals, theoretical guarantees on their width convergence rate of O(n^{-1/3}), and inference procedures for probabilistic prediction and multi-class classification, demonstrating superior local calibration compared to standard marginal conformal prediction methods.
Entities (11)
Relation Signals (10)
Michael I. Jordan → affiliatedwith → University of California, Berkeley
confidence 95% · Michael I. Jordan1... 1University of California, Berkeley
Guido Imbens → affiliatedwith → Stanford University
confidence 95% · Guido Imbens2... 2Stanford University
Xuelin Yang → affiliatedwith → University of California, Berkeley
confidence 95% · Xuelin Yang1... 1University of California, Berkeley
Baihe Huang → affiliatedwith → University of California, Berkeley
confidence 95% · Baihe Huang1... 1University of California, Berkeley
Yilong Hou → affiliatedwith → University of California, Berkeley
confidence 95% · Yilong Hou1... 1University of California, Berkeley
PICPIs → satisfies → Self-consistency condition
confidence 95% · Formally, a PICPI is an interval I satisfying a self-consistency condition: E[Y | p(X) in I] in I
Michael I. Jordan → affiliatedwith → Inria Paris
confidence 90% · Michael I. Jordan1,3... 3Inria Paris
Guido Imbens → fundedby → Office of Naval Research
confidence 90% · Imbens’ work is supported by the Office of Naval Research under Grant N00014-17-1-2131
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:A classical question in statistics is which observable quantities to condition on when drawing inferences about unobservable targets. For conformal prediction in nonparametric uncertainty quantification, standard marginal validity offers limited resolution at the prediction values on which decisions are based, and fully conditional guarantees with respect to the covariates are provably unattainable. We address this gap by introducing a prediction-based conditioning framework that we refer to as Prediction-Interval-Conditional Prediction Intervals (PICPIs). Formally, a PICPI is an interval $I$ satisfying a self-consistency condition: $$\mathbb{E} [Y \mid p(X) \in I] \in I,$$ for predictive model $p$, contextual covariate $X$, and outcome $Y$. Thus, an interval simultaneously defines a stratum of prediction values and certifies that the mean outcome in that stratum lies in the same interval. This self-consistency condition yields data-adaptive strata without altering the original prediction. Such intervals can be constructed using practical algorithms. Under regularity of the prediction distribution, the constructed intervals cover all but an arbitrarily small fraction of prediction values and have widths that decrease at rate $n^{-1/3}$, up to logarithmic factors and the prediction error. Moreover, identifying these locally calibrated intervals can, in turn, inform downstream decision-making. We derive inference procedures for PICPIs in probabilistic prediction and multi-class classification, accompanied by theoretical guarantees. Empirical results are provided that compare PICPIs with existing interval-based baselines.
Tags
Links
- Source: https://arxiv.org/abs/2609.25388v1
- Canonical: https://arxiv.org/abs/2609.25388v1
Trouble viewing inline? Open PDF directly →
Full Text
141,023 characters extracted from source content.
Expand or collapse full text
PICPIs: Prediction-Interval-Conditional Prediction Intervals Xuelin Yang1, Baihe Huang1,11footnotemark: 1 Yilong Hou1,11footnotemark: 1 Guido Imbens2 Michael I. Jordan1,3 1University of California, Berkeley 2Stanford University 3Inria Paris †thanks: Equal contribution. Imbens’ work is supported by the Office of Naval Research under Grant N00014-17-1-2131; Jordan’s work is supported by European Union under Grant ERC-2022-SYG-OCEAN-101071601. Email correspondence: xuelin@berkeley.edu. Code: https://github.com/xyang23/picpi. Abstract A classical question in statistics is which observable quantities to condition on when drawing inferences about unobservable targets. For conformal prediction in nonparametric uncertainty quantification, standard marginal validity offers limited resolution at the prediction values on which decisions are based, and fully conditional guarantees with respect to the covariates are provably unattainable. We address this gap by introducing a prediction-based conditioning framework that we refer to as Prediction-Interval-Conditional Prediction Intervals (PICPIs). Formally, a PICPI is an interval I satisfying a self-consistency condition: [Y∣p(X)∈I]∈I,E[Y p(X)∈ I]∈ I, for predictive model p, contextual covariate X, and outcome Y. Thus, an interval simultaneously defines a stratum of prediction values and certifies that the mean outcome in that stratum lies in the same interval. This self-consistency condition yields data-adaptive strata without altering the original prediction. Such intervals can be constructed using practical algorithms. Under regularity of the prediction distribution, the constructed intervals cover all but an arbitrarily small fraction of prediction values and have widths that decrease at rate n−1/3n^-1/3, up to logarithmic factors and the prediction error. Moreover, identifying these locally calibrated intervals can, in turn, inform downstream decision-making. We derive inference procedures for PICPIs in probabilistic prediction and multi-class classification, accompanied by theoretical guarantees. Empirical results are provided that compare PICPIs with existing interval-based baselines. Keywords: Prediction-conditioned inference; conformal prediction; uncertainty quantification; calibration; distribution-free inference 1. Introduction A classical problem in statistical inference is that of conditioning—on what observable quantities should one condition in order to make inferences about unobservable quantities (Fisher, 1934; Cox, 1958; Cox, 1988; Reid, 1995)? This problem has received renewed attention in recent work on nonparametric uncertainty quantification, specifically the area of conformal prediction, which aims to produce confidence intervals for nonparametric prediction (Vovk et al., 2005; Shafer and Vovk, 2008). Consider for example the following formal model for binary prediction. Let X∈X denote a vector of covariates, let Y∈0,1Y∈\0,1\ denote an outcome, and let p:→[0,1]p:X→[0,1] be a fitted model whose output p(X)p(X) is intended to approximate the conditional probability p∗(X):=ℙ(Y=1∣X)p^*(X):=P(Y=1 X). Assuming exchangeability for the observable data, (Xi,Yi)i=1n\(X_i,Y_i)\_i=1^n, but making no other structural or probabilistic assumptions, conformal prediction yields distribution-free prediction sets that come with a marginal guarantee of the form: ℙ(Y∈C(X))≥1−α, (Y∈ C(X) )≥ 1-α, (1) where C(X)C(X) is a prediction set constructed from p(X)p(X), and α∈(0,1)α∈(0,1) is a user-specified failure probability. As has been emphasized in recent years, this marginal guarantee is unconditional with respect to the covariate X, and thus does not provide direct assessment of uncertainty about a specific prediction; moreover, fully conditional coverage with respect to X is generally unattainable without strong distributional assumptions. More formally, Lei and Wasserman (2014) and Barber et al. (2021a) demonstrate that any prediction interval C(X)C(X) satisfying the strict conditional validity guarantee ℙ(Y∈C(x)∣X=x)≥1−αfor almost all x, (Y∈ C(x) X=x )≥ 1-α almost all x, (2) must be vacuous. A significant effort has been made to temper the message in this result by enforcing conditional or group-wise validity via coarsening of the covariate (Ding et al., 2023; Martinez Gil et al., 2024; Gibbs et al., 2025; Duchi, 2025). We present an alternative perspective on the problem of obtaining conditional estimates of uncertainty. Rather than conditioning on the covariate X, we suggest conditioning on the prediction p(X)p(X). The idea of prediction-based conditioning is not new—indeed, it is the foundation of the literature on probabilistic calibration (Dawid, 1982; Foster and Vohra, 1998). Recall that a calibrated model is defined as one that satisfies the following conditional criterion: ℙ(Y=1∣p(X)=v)=v. (Y=1 p(X)=v )=v. (3) That is, among all instances where the predicted probability is approximately v, the outcome should occur with frequency approximately v. Our approach to conditional uncertainty quantification is built on a similar “self-consistency” property. Specifically, we ask that the probability of a positive outcome, conditioning on the model’s prediction falling within a specific interval I, also lies within I. Formally, we have: ℙ(Y=1∣p(X)∈I)∈I. (Y=1 p(X)∈ I )∈ I. (4) Effectively, we are using the model’s predictions to stratify risk. We refer to this approach as Prediction-Interval-Conditional Prediction Intervals (PICPIs). As illustrated in Figure 1, each PICPI serves a dual role: it simultaneously defines the stratum of predictions being evaluated and provides a calibrated uncertainty statement for that specific stratum. In this way, uncertainty quantification is aligned with prediction, producing a summary of model reliability across the prediction distribution. In contrast to the intervals defined in Eq. (2),we show (in Theorem 3.2) that the width of the PICPI intervals converges to zero under suitable assumptions, providing increasingly precise characterizations of risk. Figure 1: Illustration of the PICPI methodology. Each blue interval IjI_j is both a stratum of predicted probabilities and a certificate that its empirical outcome frequency, shown in red, lies in IjI_j. The dashed diagonal line corresponds to perfect calibration. This example is generated under model misspecification (quadratic ground truth with a linear fitted model). We make the following two distinctions in comparison with standard methodology in calibration. First, calibration is usually assessed via binning. Binning schemes are typically based on fixed, external partitions of the prediction space, with overall calibration error computed as an aggregation of the errors within each bin. While useful for general model assessment, fixed bins need not align with the decision boundaries that arise when optimizing a downstream risk function. The proposed PICPIs provide a local, data-adaptive diagnostics for calibration. Rather than forcing the model to conform to pre-defined bins through recalibration, we characterize regions of the prediction space where the model is “self-consistent” and thus reliable. Second, calibration error computed from binning schemes is often used to guide post-hoc recalibration procedures such as Platt scaling, isotonic regression, or self-calibration (Platt et al., 1999; van der Laan and Alaa, 2024; Berta et al., 2024). Recalibration, however, can be undesirable (Guo et al., 2017; Pleiss et al., 2017; Carriero et al., 2025). With PICPIs, a practitioner can instead identify regions of the prediction space that are directly actionable. We further demonstrate how the knowledge of these locally calibrated intervals can inform downstream decision-making, including probabilistic prediction and classification. An extended discussion of these issues is provided in Section 6. Our specific contributions are as follows: 1. We propose PICPIs as a novel approach to conditional uncertainty quantification, we design algorithms to construct PICPIs, and we provide theoretical characterizations of these intervals. In particular, we show that the PICPIs constructed via our algorithm cover all but a user-specified prediction mass have widths controlled according to a ~(n−1/3) O(n^-1/3) term plus the prediction error (Theorem 3.2). Consequently, these widths vanish when the base predictor is consistent. 2. We provide inference procedures for probabilistic prediction and multi-class classification, and present both theoretical guarantees and empirical results. Specifically, we show that probabilistic prediction using PICPIs has marginal miscoverage that is bounded by the within-stratum variation of p∗(X)p^*(X) together with the prediction mass not covered by the PICPIs (Theorem 4.1). Multi-class classification using PICPIs has worst-case class-conditional coverage certificates (Theorem 4.3), a set-size optimality result given a prescribed coverage threshold (Theorem 4.4), which allows data re-use when constructing intervals (Theorem 4.5). The remainder of the paper is organized as follows: Section 2 introduces the PICPI definition and presents its algorithmic instantiation. Section 3 derives theoretical properties of PICPIs. Section 4 presents an inference framework using PICPIs for probabilistic prediction and classification. Section 5 discusses the empirical performance of the PICPI approach in comparison with several standard interval-based baselines. Section 6 presents related work and we present our conclusions in Section 7. 2. Definition, Algorithmic Construction, and Validity Let Y∈0,1Y∈\0,1\, let X∈X , and let p:→[0,1]p:X→[0,1] be a fitted model for p∗(x):=ℙ(Y=1∣X=x)p^*(x):=P(Y=1 X=x). The model may be trained on a separate sample or supplied as an off-the-shelf predictor. A PICPI is an interval in the prediction space that both groups similar predictions and contains the group’s true underlying outcome probability. Definition 2.1 (Prediction-Interval-Conditional Prediction Interval (PICPI)). A collection ℐp=Ijj=1JI_p=\I_j\_j=1^J is called a collection of Prediction-Interval-Conditional Prediction Intervals (PICPIs) with respect to p if, for every j∈[J]j∈[J] with ℙ(p(X)∈Ij)>0P(p(X)∈ I_j)>0, ℙ(Y=1∣p(X)∈Ij)∈Ij. (Y=1 p(X)∈ I_j)∈ I_j. (5) For simplicity of exposition we focus on binary outcomes in this section. We note that the PICPI notion extends directly to regression with a bounded real-valued outcome, as we will formalize in Remark 3.6. Algorithm 1 PICPI Construction (population mode with theoretical validity) 1: procedure Calibration(data points (xi,yi)i=1n,\(x_i,y_i)\_i=1^n, predictor p, failure probability δ) ⊳ Obtain a set of PICPIs 2: Let 0,1/K,2/K,⋯,1\0,1/K,2/K,·s,1\ be a discretization of [0,1][0,1]. 3: Initialize PICPIs ℐp←∅I_p← 4: for 0≤i<j≤K0≤ i<j≤ K do 5: Let [a,b]←[i/K,j/K][a,b]←[i/K,j/K] ⊳ Evaluate the current interval 6: if [a,b][a,b] satisfies ∑i=1n(p(xi)∈[a,b])>0 _i=1^n 1(p(x_i)∈[a,b])>0 and ∑i=1nyi⋅(p(xi)∈[a,b])∑i=1n(p(xi)∈[a,b])∈[a+2log(K2/δ)∑i=1n(p(xi)∈[a,b]),b−2log(K2/δ)∑i=1n(p(xi)∈[a,b])] _i=1^ny_i· 1(p(x_i)∈[a,b]) _i=1^n 1(p(x_i)∈[a,b])∈ [a+2 (K^2/δ) _i=1^n 1(p(x_i)∈[a,b]),b-2 (K^2/δ) _i=1^n 1(p(x_i)∈[a,b]) ] then 7: ℐp←ℐp∪[a,b]I_p _p∪\[a,b]\ 8: end if 9: end for 10: Return PICPIs ℐpI_p 11: end procedure 12: 13: procedure Inference(covariate input x,x, predictor p, failure probability δ) ⊳ Find the shortest PICPI that contains p(x)p(x) 14: ℐp←Calibration((xi,yi)i=1n,p,δ)I_p← Calibration (\(x_i,y_i)\_i=1^n,p,δ ) 15: Let [a,b][a,b] be the shortest interval in ℐpI_p that contains p(x)p(x) 16: Return [a,b][a,b] 17: end procedure We now present an algorithm to construct PICPIs with an accompanying validity statement. The population-level algorithm, presented in pseudocode as Algorithm 1, comprises two procedures. Calibration returns a collection of intervals that we will show to satisfy the PICPI condition (Theorem 2.2). Once this collection has been computed, Inference reports the shortest returned interval containing a new prediction p(x)p(x). If a disjoint summary is preferred, we select a disjoint subset with maximum covered length through the DisjointCalibration procedure provided in Appendix C. The theoretical PICPI validity of the produced intervals—that is, whether they satisfy Definition 2.1 defined on the population—follows from a concentration argument. In particular, for a candidate interval I=[a,b]I=[a,b], the PICPI condition asks whether the unknown conditional mean ℙ(Y=1∣p(X)∈I)P(Y=1 p(X)∈ I) lies between a and b. We estimate this quantity by the average label among calibration observations whose predictions fall in I. Because this average could be noisy, it is not enough for it merely to lie in I: we require it to be separated from both endpoints by a concentration margin (Line 6). The margin is larger for intervals containing fewer calibration observations, precisely where the empirical average is less reliable. This ensures that every interval admitted by the algorithm satisfies the population-level PICPI property. We formalize this statement of validity in the following theorem: Theorem 2.2. Suppose (xi,yi)i=1n\(x_i,y_i)\_i=1^n are drawn i.i.d. from the distribution of (X,Y)(X,Y) and p is fixed independently of these observations. Then, with probability at least 1−δ1-δ over the calibration sample, every interval returned by Calibration in Algorithm 1 is a PICPI with respect to p. Proof sketch. We use Hoeffding’s inequality on ℙ(Y=1∣p(X)∈I)P(Y=1 p(X)∈ I) and a union bound over K2K^2 discretized candidate intervals. For the full proof see Appendix A.1. ∎ Theorem 2.2 guarantees the PICPI property under the true underlying distribution. An empirical version of the algorithm that we will discuss in a later section (Algorithm 3, Section 5.4) produces PICPIs with respect to the empirical distribution of the calibration sample. Conceptually, their relationship is analogous to the distinction between a population CDF and an empirical CDF: while Theorem 2.2 enforces the PICPI property with respect to the underlying population measure, Algorithm 3 satisfies it under the plug-in sample measure. 3. Properties of PICPIs In this section we show that, when a predictor has error ϵε and its predicted probabilities satisfy a λ-regularity condition, most predictions will be enclosed by PICPIs. As the calibration sample size grows, these intervals shrink while their union continues to cover most of the prediction mass; the rate depends on ϵε, λ, and the calibration sample size. Our result requires a regularity condition ensuring that the predicted probabilities are sufficiently spread across [0,1][0,1]. Definition 3.1 (λ-regularity). We say that a random variable Z over [0,1][0,1] is λ-regular if ℙ(Z∈A)≥λ⋅μB(A)P(Z∈ A)≥λ· _B(A) for any measurable set A, where μB _B is the Lebesgue measure on [0,1][0,1]. Theorem 3.2 (Finite-sample width guarantee for PICPIs). Consider (xi,yi)i=1m+n\(x_i,y_i)\_i=1^m+n i.i.d. data points such that n≥2−1/2λ−1log(3n/δ)n≥ 2^-1/2λ^-1 (3n/δ). Use m data points to construct the predictor p:→[0,1]p:X→[0,1] and n data points to produce PICPIs ℐpI_p via Algorithm 1. Suppose X∼μX μ, p(X)p(X) is λ-regular, and μ[|p(X)−ℙ(Y=1∣X)|]≤ϵ, _μ [ |p(X)-P(Y=1 X) | ]≤ε, for some ϵ>0ε>0. For any δ,δ1>0δ, _1>0, with probability at least 1−δ1-δ over the randomness of (xi,yi)i=1n\(x_i,y_i)\_i=1^n, ℙx∼μ(∃[a,b]∈ℐps.t.|b−a|≤C(n,ϵ,δ,δ1)andp(x)∈[a,b])≥1−δ1, _x μ (∃[a,b] _p \ s.t. \ |b-a|≤ C(n,ε,δ, _1) \ and \ p(x)∈[a,b] )≥ 1- _1, (6) where C(n,ϵ,δ,δ1)=c0log(1/δ1)⋅((log(n/δ)λn)1/3+ϵδ1) C(n,ε,δ, _1)=c_0 (1/ _1)· ( ( (n/δ)λ n )^1/3+ ε _1 ) (7) for a universal constant c0c_0. Proof sketch. First, we construct an initial collection of intervals of controlled length and desired coverage by applying Lemma 3.4 to the distribution of p(X)p(X). That is, there exist intervals [aj,bj]\[a_j,b_j]\ where the conditional mean prediction [p(X)∣p(X)∈[aj,bj]]E[p(X) p(X)∈[a_j,b_j]] is bounded away from its own endpoints. Then we transfer this property from the conditional mean prediction [p(X)∣p(X)∈[aj,bj]]E[p(X) p(X)∈[a_j,b_j]] to the conditional label probability [Y∣p(X)∈[aj,bj]]E[Y p(X)∈[a_j,b_j]] (by setting U=p∗(X)U=p^*(X) and V=p(X)V=p(X) in Lemma 3.5). Upon discarding a small amount of mass, the retained intervals satisfy the desired PICPI property. We then show that these intervals can be produced by Algorithm 1. Define A1A_1 to be the event that every interval contains sufficiently many calibration samples, and A2A_2 the event that the empirical label average close to the true conditional label probability. Under A1∩A2A_1∩ A_2, the empirical estimate is sufficiently far from both endpoints for the algorithm to accept. By Lemma 3.3 and Hoeffding’s inequality, A1∩A2A_1∩ A_2 occurs with probability at least 1−δ1-δ. We note that Eq. (7) serves as an upper bound on the widths, while obtaining adaptive bounds that exploit additional model properties or parametric structure is an open direction for future work. For the full proof see Appendix A.2. ∎ Figure 2: Demonstration of width shrinkage. For implementation details see Appendix B.5. 3.1. Supporting lemmas As shown in the proof sketch, Theorem 3.2 relies on three supporting lemmas. We present these lemmas here to illustrate the key arguments needed to establish the theorem. Lemma 3.3 (Sample size in each interval). Let ν be a λ-regular distribution over [0,1][0,1] and let [aj,bj],j∈[J][a_j,b_j],j∈[J] are be intervals such that |bj−aj|≥c|b_j-a_j|≥ c for all j∈[J]j∈[J]. For any n∈Z+n∈ Z_+, let z1,…,znz_1,…,z_n be i.i.d. samples from ν, then with probability at least 1−Jexp(−λcn/8)1-J (-λ cn/8 ) , we have that ∑i=1n(zi∈[aj,bj])≥λcn/2 _i=1^n 1(z_i∈[a_j,b_j])≥λ cn/2 holds for all j∈[J]j∈[J]. Proof sketch. Write the sample count in each interval as a sum of Bernoulli variables and lower-bound their expectations using the λ-regularity condition. The result follows from a lower-tail Chernoff bound. The full proof can be found in Appendix A.3. ∎ Lemma 3.4 (Partition of bounded conditional mean). For any ϵ,δ>0ε,δ>0, and any distribution ν over [0,1][0,1], there exist disjoint intervals [ai,bi],i∈[L][a_i,b_i],i∈[L] such that: 1. (Well-centered conditional means) The conditional mean of Z on [ai,bi][a_i,b_i] lies away from the interval boundaries by at least ϵε (except for the endpoints 00 or 11). That is, for any i∈[L]i∈[L], ai+ϵ⋅(ai≠0)≤Z∼ν[Z∣Z∈[ai,bi]]≤bi−ϵ⋅(bi≠1). a_i+ε· 1(a_i≠ 0) _Z ν[Z Z∈[a_i,b_i]]≤ b_i-ε· 1(b_i≠ 1). (8) 2. (Coverage guarantee) The union of these PICPIs covers all but a δ-fraction of the probability mass. That is, ℙZ∼ν(Z∉∪i∈[L][ai,bi])≤δ. _Z ν (Z∉ _i∈[L][a_i,b_i] )≤δ. (9) 3. (Bounded interval lengths) Each interval has a bounded length. That is, for any i∈[L]i∈[L], 4ϵ≤|bi−ai|≤56ϵ⋅log(3/δ). 4ε≤ |b_i-a_i |≤ 56ε· (3/δ). (10) Proof Sketch. We provide a constructive proof. Discretize [0,1][0,1] into K=⌊1/(4ϵ)⌋K= 1/(4ε) subintervals of length 4ϵ4ε and form a directed graph G=(V,E)G=(V,E) on V=[K]V=[K]. We draw a directed edge i→ji→ j (|i−j|=1|i-j|=1) if either the conditional mean of Z in bin i [Z∣Z∈Bi]E[Z Z∈ B_i] is closer to the neighboring endpoint on the side of bin j, or the probability mass of Z in bin i is not substantially greater than that of bin j (i.e., ℙ(Z∈Bi)≤3ℙ(Z∈Bj)P(Z∈ B_i)≤ 3P(Z∈ B_j)). By construction, every vertex has at least one outgoing edge, and there exists bidirectional edges. Group consecutive bidirectional edges into maximal blocks. Partition each block into groups of two or three vertices. Notice that the edges between blocks split into a leftward chain followed by a rightward chain. Attach them to the first and last groups of the neighboring blocks, truncated at a radius of ⌈log(3/δ)⌉ (3/δ) . This would give the desired intervals. We then show these intervals satisfy the stated properties. The boundedness of widths follows from truncation of the chain. The conditional mean property comes from the edge construction and total expectation. Finally, outside these intervals, the absence of reverse edges implies exponential mass decay, bounding total uncovered probability mass by δ. For the full proof see Appendix A.4. ∎ Lemma 3.5 (Bounds under conditional distribution shift). Let U,VU,V be random variables jointly over [0,1][0,1], such that [|U−V|]≤ϵ>0E[|U-V|]≤ε>0. Then for any δ>0δ>0, integer L and intervals [ai,bi],i=1,…,L[a_i,b_i],i=1,…,L with disjoint interiors such that ℙ(V∈∪j∈[L][aj,bj])>1−δP (V∈ _j∈[L][a_j,b_j] )>1-δ, there exists a subset i1,…,iJ\i_1,…,i_J\ of L such that 1. (Conditional distribution on selected intervals) Conditioning on a selected interval, the discrepancy of the expected value between U and V is controlled by 2ϵ/δ2ε/δ. That is, for any j∈[J]j∈[J], |[U∣V∈[aij,bij]]−[V∣V∈[aij,bij]]|≤2ϵ/δ. |E [U V∈[a_i_j,b_i_j] ]-E [V V∈[a_i_j,b_i_j] ] |≤ 2ε/δ. (11) 2. (Coverage guarantee) The intervals not selected carry at most a δ-fraction of the probability mass. That is, ℙ(V∈∪j∈[L]∖i1,…,iJ[aj,bj])≤δ. (V∈ _j∈[L] \i_1,…,i_J\[a_j,b_j] )≤δ. (12) Consequently, the union of the selected intervals has probability mass at least 1−2δ1-2δ. Proof sketch. By applying the triangle inequality to each interval w violating Eq. (11), the conditional mean absolute deviation satisfies [|U−V|∣V∈[aw,bw]]>2ϵ/δE[|U-V| V∈[a_w,b_w]]>2ε/δ. Accounting for endpoint overlaps (∑w(V∈[aw,bw])≤2 _w 1(V∈[a_w,b_w])≤ 2), we have 2ϵ≥∑w[|U−V|(V∈[aw,bw])]>ℙ(V∈⋃w[aw,bw])⋅(2ϵ/δ)2ε≥ _wE[|U-V| 1(V∈[a_w,b_w])]>P(V∈ _w[a_w,b_w])·(2ε/δ), yielding total violating mass strictly less than δ. The full proof is provided in Appendix A.5. ∎ Remark 3.6 (Extension to regression). If p is a regression predictor intended to approximate the conditional mean [Y∣X=x],E[Y X=x], Definition 2.1 becomes [Y∣p(X)∈Ij]∈Ij.E\! [Y p(X)∈ I_j ]∈ I_j. A binary outcome Y can be viewed as a special case where [Y∣p(X)∈Ij]=ℙ(Y=1∣p(X)∈Ij).E\! [Y p(X)∈ I_j ]=P\! (Y=1 p(X)∈ I_j ). For bounded outcomes where Y∈[a,b]Y∈[a,b] almost surely, the constructions in Section 2 still apply by replacing the empirical positive-outcome frequency with the empirical average of Y. The finite-sample argument in Theorem 3.2 holds via an affine rescaling of Y and p(X)p(X) from [0,1][0,1] to [a,b][a,b]. For unbounded outcomes, the definition remains valid whenever the relevant conditional expectations exist, but the distribution-free concentration guarantees developed here do not apply without additional tail or moment assumptions. 4. Inference via PICPIs We now show how to use PICPIs as building blocks for prediction sets and inference in two tasks: interval prediction and label set prediction. 4.1. Interval prediction for probabilities In this section we show how PICPIs can be used to construct intervals for probabilistic prediction with coverage guarantees. Assume that we are given a prediction model p and PICPIs ℐpI_p. Define functions l,u:[0,1]→[0,1]l,u:[0,1]→[0,1] such that the interval [l(z),u(z)][l(z),u(z)] represents a PICPI that contains z. In other words, [l(z),u(z)]∈ℐp[l(z),u(z)] _p and z∈[l(z),u(z)]z∈[l(z),u(z)]. If there is no such interval, then let [l(z),u(z)]=[0,1][l(z),u(z)]=[0,1]. Define the expansion operator EϵE^ε as follows: for any set A⊂[0,1]A⊂[0,1], Eϵ(A)=b∈[0,1]:d(b,A)≤ϵE^ε(A)=\b∈[0,1]:d(b,A)≤ε\ where d is the Euclidean distance. For x∈x , we obtain a prediction sets based on PICPIs ℐpI_p by: Ctp(x,ϵ)=Eϵ([l(p(x)),u(p(x))]). C_tp(x;ε)=E^ε ( [l(p(x)),u(p(x)) ] ). (13) This can be viewed as an ϵε-expanded version of a PICPI evaluated at p(x)p(x). The following theorem establishes on how well Ctp(x,ϵ)C_tp(x;ε) captures the true conditional probability ℙ(Y=1∣X=x)P(Y=1 X=x). Specifically, the coverage error is upper-bounded by the fluctuation of true probabilities within each interval and the uncovered prediction mass. By Lemma 3.4 the second term is at most δ. Theorem 4.1 (Marginal coverage of the true probability). Let p∗(x)=ℙ(Y=1∣X=x)p^*(x)=P(Y=1 X=x) denote the true probability that Y=1Y=1 given x. Assume we have disjoint PICPIs Ijj=1J\I_j\_j=1^J, and write mj=ℙ(p(X)∈Ij)m_j=P(p(X)∈ I_j). For any ϵ>0ε>0, it holds that ℙ(p∗(X)∉Ctp(X;ϵ))≤ϵ−2∑j:mj>0mjVar(p∗(X)∣p(X)∈Ij)⏟“in-bin” variation+(1−∑j=1Jmj)⏟predictions not covered by ℐp. (p^*(X)∉ C_tp(X;ε) )≤ε^-2 _j:m_j>0m_j (p^*(X) p(X)∈ I_j )_``in-bin′ variation+ (1- _j=1^Jm_j )_predictions not covered by I_p. The proof can be found in Appendix A.6. When the predictor and interval collection are constructed from data that is independent of the test observation, the theorem applies conditionally on those data whenever the resulting intervals satisfy its hypotheses. The bound in Theorem 4.1 can be sharpened when additional structure is available. The following remark shows that if the predictor p approximates the true probability p∗p^* uniformly or satisfies a Lipschitz-type condition, the in-bin variation term admits explicit upper bounds in terms of the interval lengths. Remark 4.2 (Corollaries under smoothness conditions). Under the assumptions of Theorem 4.1, suppose first that ‖p∗−p‖∞≤η\|p^*-p\|_∞≤η for some η≥0η≥ 0. For every positive-mass interval Ij=[aj,bj]I_j=[a_j,b_j], conditioning on p(X)∈Ijp(X)∈ I_j gives p∗(X)∈[max0,aj−η,min1,bj+η]almost surely. p^*(X)∈[ \0,a_j-η\, \1,b_j+η\] surely. The length of this range is at most |Ij|+2η|I_j|+2η. Since a random variable supported on an interval of length d has variance at most d2/4d^2/4, Theorem 4.1 yields ℙ(p∗(X)∉Ctp(X,ϵ)) (p^*(X)∉ C_tp(X;ε) ) ≤ϵ−2∑j:mj>0mjVar(p∗(X)∣p(X)∈Ij)+(1−∑j=1Jmj) \ ≤ε^-2 _j:m_j>0m_j (p^*(X) p(X)∈ I_j )+ (1- _j=1^Jm_j ) ≤14ϵ2∑j:mj>0mj(|Ij|+2η)2+(1−∑j=1Jmj). \ ≤ 14ε^2 _j:m_j>0m_j(|I_j|+2η)^2+ (1- _j=1^Jm_j ). Here ϵ>0ε>0 is the expansion radius in Eq. (13). Alternatively, suppose that |p∗(x)−p∗(x′)|≤L⋅|p(x)−p(x′)||p^*(x)-p^*(x )|≤ L·|p(x)-p(x )| for all x,x′∈x,x . When p(x),p(x′)∈Ijp(x),p(x )∈ I_j, their difference is at most |Ij||I_j|, so the conditional range of p∗(X)p^*(X) has length at most L|Ij|L|I_j|. Applying the same variance bound gives ℙ(p∗(X)∉Ctp(X,ϵ)) (p^*(X)∉ C_tp(X;ε) ) ≤ϵ−2∑j:mj>0mjVar(p∗(X)∣p(X)∈Ij)+(1−∑j=1Jmj) \ ≤ε^-2 _j:m_j>0m_j (p^*(X) p(X)∈ I_j )+ (1- _j=1^Jm_j ) ≤L24ϵ2∑j:mj>0mj|Ij|2+(1−∑j=1Jmj). \ ≤ L^24ε^2 _j:m_j>0m_j|I_j|^2+ (1- _j=1^Jm_j ). 4.2. Label set prediction for multi-class classification We now use PICPIs for label set prediction in classification with G classes. In particular, we establish three guarantees: (i) an exact finite-sample certificate for class-conditional coverage, providing the tightest worst-case miscoverage bound available from the interval information (Theorem 4.3); (i) population size-optimality, showing that our interval selection minimizes expected label-set size among deterministic unions of the given classwise intervals satisfying the prescribed class-conditional coverage guarantees (Theorem 4.4); and (i) a finite-sample same-sample guarantee, showing that the empirically feasible procedure retains valid class-conditional coverage when the same calibration data are used to construct the intervals and estimate their masses (Theorem 4.5). Let p:→[0,1]Gp:X→[0,1]^G be a prediction model that outputs class probabilities, where its g-th entry pg(x)p_g(x) predicts ℙ(Y(X)=g∣X=x)P(Y(X)=g X=x). For each class g∈[G]g∈[G], we construct disjoint PICPIs ℐpg=[ag,j,bg,j]j=1JgI_p_g= \[a_g,j,b_g,j] \_j=1^J_g for the one-vs-rest problem 1Y(X)=g1\Y(X)=g\ using the interval-construction procedure in Section 2. Suppose we select a subset of interval indices Sg⊆[Jg]S_g [J_g] for class g (by some selection rules of our choice), and denote their collection as S=(S1,…,SG)S=(S_1,…,S_G). For any given S, we define the classwise acceptance region as the union of intervals in SgS_g: Ag(Sg):=⋃j∈Sg[ag,j,bg,j].A_g(S_g):= _j∈ S_g[a_g,j,b_g,j]. Given input x, we predict a set of labels which collects the classes whose predicted probabilities fall within their respective acceptance regions: Cmc(x,S)=g∈[G]:pg(x)∈Ag(Sg). C_mc(x;S)= \g∈[G]:p_g(x)∈ A_g(S_g) \. (14) For example, SgS_g can be constructed by a simple upper-endpoint threshold rule with threshold τ: Sg(τ):=j∈[Jg]:bg,j≥τ,S(τ):=(S1(τ),…,SG(τ)), S_g(τ):= \j∈[J_g]:b_g,j≥τ \, S(τ):= (S_1(τ),…,S_G(τ) ), (15) so that Cmc(x;S(τ))=g∈[G]:pg(x)∈⋃j:bg,j≥τ[ag,j,bg,j]. C_mc (x;S(τ) )= \g∈[G]:p_g(x)∈ _j:\,b_g,j≥τ[a_g,j,b_g,j] \. (16) For notational convenience, define mg,j:=ℙ(pg(X)∈[ag,j,bg,j]),πg:=ℙ(Y(X)=g),m_g,j:=P (p_g(X)∈[a_g,j,b_g,j] ), _g:=P (Y(X)=g ), tg:=1−∑j=1Jgmg,j,Δpg:=max1≤j≤Jg(bg,j−ag,j).t_g:=1- _j=1^J_gm_g,j, _p_g:= _1≤ j≤ J_g (b_g,j-a_g,j ). Throughout, we assume πg>0 _g>0 for every class g∈[G]g∈[G]. The following theorem gives an exact class-conditional coverage certificate for any rule of the form Cmc(x,S)C_mc(x;S). Theorem 4.3 (Class-conditional coverage of the true label). Fix a class g∈[G]g∈[G] and a subset Sg⊆[Jg]S_g [J_g]. It holds that ℙ(Y(X)∉Cmc(X,S)∣Y(X)=g)≤Γg(Sg), (Y(X)∉ C_mc(X;S) Y(X)=g )≤ _g(S_g), (17) where Γg(Sg):=minΓgup(Sg),Γglow(Sg), _g(S_g):= \ _g^up(S_g), _g^low(S_g) \, and Γgup(Sg):=tg+∑j∉Sgmg,jbg,jπg,Γglow(Sg):=1−∑j∈Sgmg,jag,jπg, _g^up(S_g):= t_g+ _j∉ S_gm_g,jb_g,j _g, _g^low(S_g):=1- _j∈ S_gm_g,ja_g,j _g, Moreover, this bound is exact in the following sense: for fixed πg,tg,mg,j,ag,j,bg,jj=1Jg. _g, t_g, \m_g,j,a_g,j,b_g,j \_j=1^J_g. the quantity Γg(Sg) _g(S_g) is the worst-case class-conditional miscoverage over all one-vs-rest laws satisfying the PICPIs definition ℙ(Y(X)=g|pg(X)∈[ag,j,bg,j])∈[ag,j,bg,j]for all j.P (Y(X)=g\, |\,p_g(X)∈[a_g,j,b_g,j] )∈[a_g,j,b_g,j] all j. For the proof see Appendix A.7. Population-optimal label set. Ideally, we want our predicted label set to be as small as possible while still guaranteeing that the true label is caught within our error budget. To formalize this trade-off, let α=(α1,…,αG)∈[0,1]Gα=( _1,…, _G)∈[0,1]^G denote a vector of target class-conditional miscoverage levels specified by the user. For each class g∈[G]g∈[G], we find the optimal subset of prediction intervals to keep by minimizing the expected size of the prediction set, subject to the error guarantee. Formally, we define the population-optimal retained interval set as: Sg⋆(αg)∈argmin∑j∈S⊆[Jg]mg,jsubject toΓg(S)≤αg. S_g ( _g)∈ _S [J_g] _j∈ Sm_g,j to _g(S)≤ _g. (18) The resulting population-optimal label set simply collects all classes whose predicted probabilities fall into these optimally selected intervals: Cmc⋆(x,α):=Cmc(x,S⋆(α)),S⋆(α):=(S1⋆(α1),…,SG⋆(αG)).C_mc (x;α):=C_mc (x;S (α) ), S (α):= (S_1 ( _1),…,S_G ( _G) ). Because Γg(S)=minΓgup(S),Γglow(S), _g(S)= \ _g^up(S), _g^low(S) \, the feasible set in (18) is the union of two simpler feasible sets. Equivalently, one may solve Sgkeep(αg)∈argminS⊆[Jg] S_g^keep( _g)∈ _S [J_g] ∑j∈Smg,j _j∈ Sm_g,j (19) subject to ∑j∈Sag,jmg,j≥(1−αg)πg, _j∈ Sa_g,jm_g,j≥(1- _g) _g, and Sgdrop(αg)∈argminS⊆[Jg] S_g^drop( _g)∈ _S [J_g] ∑j∈Smg,j _j∈ Sm_g,j (20) subject to tg+∑j∉Smg,jbg,j≤αgπg, t_g+ _j∉ Sm_g,jb_g,j≤ _g _g, and then keep whichever feasible solution has the smaller objective value. Algorithm 2 Constructing population-optimal label set 1: Disjoint PICPIs ℐpg=[ag,j,bg,j]j=1JgI_p_g=\[a_g,j,b_g,j]\_j=1^J_g, population interval masses mg,j\m_g,j\, class priors πg\ _g\, and targets α=(α1,…,αG)α=( _1,…, _G) 2: Population-optimal label set Cmc⋆(⋅,α)C_mc (·;α) 3: for g=1g=1 to G do 4: Compute tg=1−∑j=1Jgmg,jt_g=1- _j=1^J_gm_g,j. 5: Solve the keep-side problem (19); if infeasible, assign objective value +∞+∞. 6: Solve the drop-side problem (20); if infeasible, assign objective value +∞+∞. 7: Let Sg⋆(αg)S_g ( _g) be the feasible solution with smaller objective value. 8: end for 9: return Cmc⋆(x,α)=g∈[G]:pg(x)∈⋃j∈Sg⋆(αg)[ag,j,bg,j]C_mc (x;α)= \g∈[G]:p_g(x)∈ _j∈ S_g ( _g)[a_g,j,b_g,j] \. Theorem 4.4 (Size-optimality at fixed guaranteed coverage). Let Cmc⋆(x,α)C_mc (x;α) be defined by (18). Then, for every class g∈[G]g∈[G], ℙ(Y(X)∉Cmc⋆(X,α)∣Y(X)=g)≤αg.P (Y(X)∉ C_mc (X;α) Y(X)=g )≤ _g. Moreover, [|Cmc⋆(X,α)|]=min[|Cmc(X,S)|]:Γg(Sg)≤αg for all g∈[G].E [ |C_mc (X;α) | ]= \E [|C_mc(X;S)| ]: _g(S_g)≤ _g for all g∈[G] \. For the proof see Appendix A.8. Empirically feasible version of the algorithm. The population-optimal algorithm presented above is not directly implementable because the interval masses mg,jm_g,j and the class priors πg _g are unknown. To obtain a practical algorithm we thus replace them by one-sided confidence bounds computed from the same calibration set that is already used to construct the PICPIs. Let Dcal=(Xk,Yk)k=1ncalD_cal= \(X_k,Y_k) \_k=1^n_cal be the calibration set, and let ℐ^pg=I^g,j=[a^g,j,b^g,j]j=1J^g I_p_g= \ I_g,j=[ a_g,j, b_g,j] \_j=1 J_g be the resulting disjoint PICPIs for class g. Let pgA_p_g denote the finite candidate interval family searched by Algorithm 2 for class g, and define Ncand:=∑g=1G|pg|.N_cand:= _g=1^G |A_p_g |. For each selected interval, define the empirical mass m^g,j:=1ncal∑k=1ncal(pg(Xk)∈I^g,j),π^g:=1ncal∑k=1ncal(Yk=g). m_g,j:= 1n_cal _k=1^n_cal 1 (p_g(X_k)∈ I_g,j ), π_g:= 1n_cal _k=1^n_cal 1(Y_k=g). For confidence levels δm,δπ∈(0,1) _m, _π∈(0,1), let εm:=log(2Ncand/δm)2ncal,επ:=log(2G/δπ)2ncal. _m:= (2N_cand/ _m)2n_cal, _π:= (2G/ _π)2n_cal. We then define m¯g,j:=(m^g,j−εm)+,π¯g:=(π^g−επ)+,π¯g:=(π^g+επ)∧1. m_g,j:= ( m_g,j- _m )_+, π_g:= ( π_g- _π )_+, π_g:= ( π_g+ _π ) 1. For any subset Sg⊆[J^g]S_g [ J_g], define Γ^gup(Sg):=1−∑j∈Sgm¯g,j−∑j∉Sg(1−b^g,j)m¯g,jπ¯g,π¯g>0,+∞,π¯g=0, _g^up(S_g):= cases 1- _j∈ S_g m_g,j- _j∉ S_g(1- b_g,j) m_g,j π_g,& π_g>0,\\[12.91663pt] +∞,& π_g=0, cases Γ^glow(Sg):=1−∑j∈Sga^g,jm¯g,jπ¯g,Γ^g(Sg):=minΓ^gup(Sg),Γ^glow(Sg). _g^low(S_g):=1- _j∈ S_g a_g,j m_g,j π_g, _g(S_g):= \ _g^up(S_g), _g^low(S_g) \. This yields the empirical relaxation S^g⋆(αg)∈argmin∑j∈S⊆[J^g]m^g,jsubject toΓ^g(S)≤αg, S_g ( _g)∈ _S [ J_g] _j∈ S m_g,j to _g(S)≤ _g, and the resulting feasible label set C^mc⋆(x,α):=g∈[G]:pg(x)∈⋃j∈S^g⋆(αg)I^g,j. C_mc (x;α):= \g∈[G]:p_g(x)∈ _j∈ S_g ( _g) I_g,j \. The key point is that the same calibration sample may be used both to construct the intervals and to estimate their masses: because Algorithm 2 only searches over a fixed finite candidate family, uniform concentration over that family automatically controls the masses of the finally selected intervals as well. Theorem 4.5 (Fully feasible same-sample certificate). Suppose Algorithm 2 is run for each class g∈[G]g∈[G] with classwise failure budgets summing to δPI _PI, so that with probability at least 1−δPI1- _PI all returned intervals in ℐ^pg=I^g,j=[a^g,j,b^g,j]j=1J^g,g∈[G], I_p_g= \ I_g,j=[ a_g,j, b_g,j] \_j=1 J_g, g∈[G], are genuine disjoint PICPIs simultaneously. For any S=(S1,…,SG)S=(S_1,…,S_G) with Sg⊆[J^g]S_g [ J_g], define C^mc(x,S)=g∈[G]:pg(x)∈⋃j∈SgI^g,j. C_mc(x;S)= \g∈[G]:p_g(x)∈ _j∈ S_g I_g,j \. Then, with probability at least 1−δPI−δm−δπ1- _PI- _m- _π over the draw of DcalD_cal, the following holds simultaneously for all classes g∈[G]g∈[G] and all subsets Sg⊆[J^g]S_g [ J_g]: ℙ(Y(X)∉C^mc(X,S)|Y(X)=g,Dcal)≤Γ^g(Sg).P (Y(X)∉ C_mc(X;S)\ |\ Y(X)=g,\ D_cal )≤ _g(S_g). In particular, whenever S^g⋆(αg) S_g ( _g) is feasible for every class, ℙ(Y(X)∉C^mc⋆(X,α)|Y(X)=g,Dcal)≤αgfor all g∈[G].P (Y(X)∉ C_mc (X;α)\ |\ Y(X)=g,\ D_cal )≤ _g all g∈[G]. The proof can be found in Appendix A.9. 5. Experiments We compare PICPI-based methods with other interval-prediction methods for predicting probabilities (Task 1) and multiclass classification (Task 2). In both tasks, point predictors are fitted on a training split, intervals (or label sets) are constructed using a separate calibration split, and the resulting intervals (or label sets) are evaluated on a test split. 5.1. Task 1: Interval prediction for probabilities. As described in Section 4.1, the goal is to produce an interval to predict ℙ(Y=1∣X=x)P(Y=1 X=x). Each baseline produces a set of intervals ℐI. That is, for a predictor p p, we obtain the interval I∈ℐI containing predicted p^(x) p(x) for an input x. The data-generating processes (DGPs) use 1-dimensional and 20-dimensional covariates with binary targets. Details of the target probabilities and label-generating mechanisms are provided in Appendix B.3. 5.2. Task 2: Label set prediction for multi-class classification. As described in Section 4.2, the goal is to produce a label set. For each class, interval-based baselines produce an interval set ℐ(g)I^(g), from which we select the interval containing p^g(x) p_g(x). For DGPs, target labels are drawn as Y|X=x∼Categorical(q(x))Y X=x (q(x)). The Linear-Gaussian DGP uses isotropic Gaussian covariates X∼(0,I10)X (0,I_10) with sparse linear softmax probabilities q(x)=softmax(Wx+b)q(x)=softmax(Wx+b). The Nonlinear-Mixture DGP generates non-Gaussian, multimodal covariates via a four-component latent-factor mixture. Its class probabilities scale the linear logits by a bounded covariate-dependent function s(x)s(x) such that q(x)=softmax((Wx+b)/s(x))q(x)=softmax ((Wx+b)/s(x) ). Details are provided in Appendix B.4. 5.3. Baselines We compare the following baselines: (i) PICPI constructs intervals from the calibration split using the same underlying procedure as our method, with an ϵε-expansion for Task 1 and classwise thresholding of upper endpoints for Task 2. (i) Simultaneous confidence intervals construct pointwise intervals by applying the Delta method to a Wald confidence region for the parameters of a parametric probability model. (i) Split conformal prediction forms symmetric intervals around the fitted probability using a calibration quantile of absolute residuals, providing marginal coverage guarantees. (iv) Calibration-based intervals discretize [0,1][0,1] into bins and retain bins whose mass-weighted calibration error is comparable to that of PICPI, assigning each prediction to the closest retained bin when necessary. (v) Fixed-width binning discretizes [0,1][0,1] into fixed-width bins and reports the bin containing the predicted probability. (vi) Conformal classification, used only for Task 2, constructs prediction sets by including classes whose predicted probabilities exceed a calibration-based threshold rather. This is not an interval-based method, but has marginal coverage of the true label. We extend their formulation details and metrics in Appendix B.1 and B.2, respectively. 5.4. PICPIs with empirical self-consistency Algorithm 3 PICPI Construction (empirical mode with self-consistency on observed data) 1: procedure Calibration(data points (xi,yi)i=1n,\(x_i,y_i)\_i=1^n, predictor p) ⊳ Obtain a set of PICPIs with respect to the empirical distribution 2: Compute predictions pi←p(xi)p_i← p(x_i) and sort pip_i in ascending order. Let pi′p_i and yi′y_i to denote the sorted prediction and outcome respectively. 3: Let t0←0t_0← 0. 4: for j=1,2,…,j=1,2,…, do 5: Find mintj∈(tj−1,1] t_j∈(t_j-1,1] subject to ⊳ find the next endpoint on a grid of [0,1] ∑i=1nyi′(pi′∈(tj−1,tj])∑i=1n(pi′∈(tj−1,tj])∈ _i=1^ny_i 1(p_i ∈(t_j-1,t_j]) _i=1^n 1(p_i ∈(t_j-1,t_j])∈ (tj−1,tj]. \ (t_j-1,t_j]. ∑i=1nyi′(pi′∈(tj,1])∑i=1n(pi′∈(tj,1])∈ _i=1^ny_i 1(p_i ∈(t_j,1]) _i=1^n 1(p_i ∈(t_j,1])∈ (tj,1]. \ (t_j,1]. 6: Let Ij←(tj−1,tj]I_j←(t_j-1,t_j]. 7: if tj=1t_j=1 then 8: Let J←jJ← j. Break. 9: end if 10: end for 11: Return Ijj=1J\I_j\_j=1^J 12: end procedure 13: 14: procedure Inference(covariate input x,x, predictor p) ⊳ Find the PICPI that contains p(x)p(x) 15: Ijj=1J←Calibration((xi,yi)i=1n,p)\I_j\_j=1^J← Calibration (\(x_i,y_i)\_i=1^n,p ) 16: Let I be the interval among Ijj=1J\I_j\_j=1^J that contains p(x)p(x) 17: Return I 18: end procedure We now present Algorithm 3, the empirical variant of algorithm used in experiments in this section. It returns PICPIs with respect to the empirical distribution of the calibration sample. We refer to this property as empirical self-consistency to distinguish it from the population validity guaranteed by Theorem 2.2. The key difference is that Algorithm 3 omits the concentration margin (Line 6, Algorithm 1) when evaluating intervals. While this omission yields narrower, more practical intervals, the theoretical population validity guaranteed by Theorem 2.2 no longer applies. Figure 3 illustrates the resulting tradeoff. Most proposed intervals also satisfy the PICPI condition on held-out data, whereas the shortest intervals (of length at most 0.030.03) are more prone to false acceptance. The algorithm proceeds as follows: starting at the left endpoint, it chooses the smallest next endpoint for which the empirical label mean of the new bin lies in that bin, and so does the remaining part in [0,1][0,1]. It then repeats this step on the remainder of [0,1][0,1]. At inference time, a new prediction is assigned to its unique bin in the resulting partition. Figure 3: Demonstration of empirical version (Algorithm 3). Left: Proportion of proposed intervals satisfying the PICPI property on a held-out dataset. Middle: With more calibration data, the empirical version admits narrower intervals. Right: While most intervals satisfy the PICPI property, the narrowest intervals (length ≤0.03≤ 0.03) at the bottom exhibit lower pass rates. Their pass rate follows a U-shaped pattern as the amount of calibration data increases: initially, additional data allow the procedure to admit increasingly narrow and more challenging intervals, leading to a lower pass rate; with sufficiently large calibration samples, however, the increased information improves calibration even for these narrow intervals, and the pass rate recovers. For a formal presentation of the setup studied here, see Appendix B.3. Figure 4: Visualization under mildly miscalibrated model. PICPI adapts its interval endpoints to local prediction bins to account for deviations between the fitted and target probabilities, while maintaining substantially shorter intervals than split conformal prediction. The simultaneous confidence interval is less responsive to bin-specific deviations from calibration, while the binning-based baselines have poor coverage of the target probability. Figure 5: Task 1 demonstration with a misspecified decision-tree point estimator. PICPIs adapt to biased, piecewise-constant predictions by adjusting interval endpoints and widening where calibration reveals instability or mismatch. The smooth curve shows the true conditional probability p∗(x)p^*(x), while the step function shows the misspecified decision-tree estimate p^(x) p(x). The dashed step curves mark the lower and upper endpoints of the PICPI interval, and the shaded boxes show the resulting interval assigned at each x. 5.5. Task 1 Results Figures 4 and 5, together with Table 1, show that PICPI achieves a favorable balance among length, calibration, and coverage measures. In the univariate visualization, PICPI adapts to the local score partition and tracks the nonlinear target probability p∗(x)p^*(x) more faithfully than the two binning baselines, while split conformal is visibly conservative. In the misspecified decision-tree example of Figure 5, the PICPIs remain responsive when the fitted predictor is structurally biased: the decision tree produces a piecewise-constant p^(x) p(x) that fails to track the nonlinear true probability curve, but PICPI adapts its interval endpoints around these score plateaus and widens in regions where the calibration data reveal instability or mismatch. Quantitatively, Table 1 shows that PICPI attains the lowest calibration errors while maintaining substantially shorter intervals than the simultaneous confidence interval and split conformal prediction, and it achieves much higher true-probability coverage than the calibration-based and fixed-width binning baselines. It avoids the over-conservatism of conformal-style wide intervals while retaining far better coverage and calibration behavior than binning alternatives. In particular, its calibration errors indicate that the reported intervals are well centered around the true probabilities and the its binning follows the calibration structure of predictor. Method Avg. Length Rel. ECEmid_mid Rel. ECEpred_pred Coverage over [0,1] Coverage of p∗p^* Simultaneous CI 0.41 ± 0.01 6.39 3.25 – 100.00% ± 0.00% PICPI 0.11 ± 0.01 1.00 1.00 100.00% ± 0.00% 79.47% ± 4.68% Split conformal prediction 0.87 ± 0.01 27.02 2.43 – 100.00% ± 0.00% Calibration-based interval 0.02 ± 0.00 1.16 1.01 60.78% ± 11.49% 28.07% ± 3.25% Fixed-width binning 0.02 ± 0.00 1.05 1.11 100.00% ± 0.00% 29.28% ± 3.34% Table 1: Comparison of interval length, calibration errors, and coverage across methods. PICPI achieves substantially shorter intervals than simultaneous and conformal methods while maintaining full coverage of the prediction space and improved coverage of the target probability. ECEmidECE_mid measures if interval midpoint is representative of the true probabilities of the predictions assigned to it. ECEpredECE_pred measures if average prediction agree with the true probability among the points receiving this interval. It evaluates how well the interval grouping preserves the calibration of p(X)p(X). Definitions of the metrics are supplemented in Appendix B.2. Covariate dimension d=20d=20. 5.6. Task 2 Results Figure 6: Task 2 results, averaged over 10001000 replications. First row: linear-Gaussian DGP. Second row: nonlinear-mixture softmax DGP. Each replication uses disjoint splits of size ntrain=ncal=neval=2000n_train=n_cal=n_eval=2000, with d=10d=10 covariates and G=4G=4 classes. PICPI maintains near-target worst-class coverage in both settings. This advantage is more pronounced under the more complex nonlinear-mixture DGP, where competing baselines show substantial degradation in worst-class coverage. Implementation details are given in Appendix B.4. Across both the linear-Gaussian and nonlinear-mixture softmax DGPs, PICPI delivers the most stable worst-class coverage relative to prediction-set size, as illustrated in the first column of Figure 6. This advantage becomes increasingly pronounced as the complexity of the DGP grows. The second and third columns of Figure 6 demonstrate that under the nonlinear-mixture DGP, the performance of competing methods deteriorates substantially. In contrast, PICPI continues to track the target coverage closely while maintaining the smallest prediction sets. Note that the average class coverage shown in the second column is distinct from marginal coverage–standard conformal classification, despite its marginal coverage guarantees, suffers from poor worst-class coverage. Conversely, PICPI consistently yields class-wise coverage that is lower-bounded by the target level. These empirical observations align exactly with Theorems 4.3, 4.4, and 4.5, suggesting that PICPI is robust to miscalibration resulting from nonlinear covariate structures and model misspecification. 6. Related Work Conformal prediction. Conformal prediction constructs set-valued predictors C(X)C(X) with finite-sample marginal validity, ℙ(Y∈C(X))≥1−α,P (Y∈ C(X) )≥ 1-α, under exchangeability (Vovk et al., 2005; Shafer and Vovk, 2008). Split, cross-, and jackknife-style variants—including the jackknife+, jackknife-after-bootstrap, conformalized quantile regression, and conformal predictive distributions—show that the framework can be layered on top of a broad range of base learners (Barber et al., 2021b; Kim et al., 2020; Romano et al., 2019; Vovk, 2015; Vovk et al., 2018). A stronger target is conditional coverage, ℙ(Y∈C(X)∣X=x)≥1−αfor a.e. x,P (Y∈ C(X) X=x )≥ 1-α a.e. x, but distribution-free versions of this guarantee are impossible except through vacuous sets in the worst case (Lei and Wasserman, 2014; Barber et al., 2021a). Recent work therefore studies weaker conditional notions obtained by coarsening the conditioning event: Mondrian or group-conditional validity conditions on a prespecified partition or label (Vovk et al., 2003a; Vovk, 2012; Martinez Gil et al., 2024); multivalid coverage asks for simultaneous approximate validity over a rich family of subpopulations (Jung et al., 2023; Gibbs et al., 2025); and recent analyses clarify sample- and training-conditional targets and when they can be attained under additional assumptions (Duchi, 2025; Liang and Barber, 2025). Related extensions handle covariate shift, label shift, online adaptation, and broader forms of nonexchangeability (Tibshirani et al., 2019; Gibbs and Candes, 2021; Podkopaev and Ramdas, 2021; Barber et al., 2022; Berta et al., 2024). The PICPI target, ℙ(Y=1∣p(X)∈I)∈I,P (Y=1 p(X)∈ I )∈ I, belongs to this “nearly conditional” regime, but the conditioning variable is the model score itself rather than the full feature vector. This distinction is central for our results: the score is exactly the quantity used downstream for thresholding, classification, and decision-making, so conditioning on p(X)p(X) avoids the impossibility of full covariate-conditional validity while remaining aligned with the deployed predictor. Theorem 4.1 converts this score-conditional statement into a conformal-style guarantee for the latent mean p∗(X)p (X), and Theorem 3.2 supplies an efficiency statement: under λ-regularity and L1L_1 score error ε , for all but a δ1 _1 fraction of score values there exists a PICPI of width of order (logn/n)1/3+ε/δ1( n/n)^1/3+ / _1. This complements a separate line of work that studies the length, volume, or oracle optimality of conformal regions. For density-level-set predictors, Lei et al. (2013) establish finite-sample validity together with asymptotic efficiency relative to oracle level sets; Lei et al. (2018) show that split conformal regression inherits consistency and can adapt interval length to heteroscedasticity; Chernozhukov et al. (2021) and Izbicki et al. (2020); Izbicki et al. (2022) obtain approximate conditional validity and convergence to oracle conditional intervals or highest-density regions when conditional distribution estimators are consistent; and Yang and Kuchibhotla (2025); Bars and Humbert (2025) quantify finite-sample width losses relative to oracle choices. Relative to these papers, Theorem 3.2 studies the shrinking object most relevant for our problem: not a response interval around Y, but a score-conditional calibration interval that certifies the meaning of the model’s own probability output. Calibration. For a binary score p(X)∈[0,1]p(X)∈[0,1], perfect calibration means ℙ(Y=1∣p(X)=v)=v,P (Y=1 p(X)=v )=v, for all relevant v (Dawid, 1982; Foster and Vohra, 1998). In practice this pointwise condition is usually assessed through coarser summaries: reliability diagrams compare empirical frequencies to predicted probabilities on bins of the score, and scalar criteria such as ECE aggregate these discrepancies across bins (Gneiting and Raftery, 2007; Vaicenavicius et al., 2019; Carrell et al., 2022; Blasiok and Nakkiran, 2024). A large post-hoc calibration literature then modifies the original score via parametric or nonparametric maps—Platt scaling, isotonic regression, Bayesian binning, temperature scaling, beta calibration, and Dirichlet calibration—to improve this global agreement (Platt et al., 1999; Zadrozny and Elkan, 2002; Niculescu-Mizil and Caruana, 2005; Naeini et al., 2015; Guo et al., 2017; Kull et al., 2017; Kull et al., 2019). Our objective is different: we do not alter p, but instead seek data-adaptive subsets of the score range on which the conditional mean outcome is certified to lie inside the same interval. In this sense, classical conformal prediction controls the random event Y∈C(X)\Y∈ C(X)\, whereas PICPI applies conformal-style uncertainty quantification to the conditional mean outcome along the score axis. This places PICPI between descriptive calibration diagnostics and conformal uncertainty quantification. Multicalibration asks for simultaneous approximate calibration over many subgroups and is closely tied to fairness constraints (Hébert-Johnson et al., 2018; Kleinberg et al., 2016; Pleiss et al., 2017). PICPI can be read as a one-dimensional, score-indexed analogue of this idea: the conditioning sets are subsets of the score axis rather than externally specified demographic groups, and the guarantee is interval-valued rather than an additive calibration error bound. Vovk and collaborators developed self-calibrating forecasters, Venn predictors, and Venn–Abers methods that output calibrated probabilities or probability intervals under exchangeability (Vovk et al., 2003b; Vovk and Petej, 2014). More recently, Gupta et al. (2020) give distribution-free confidence intervals for calibrated binary probabilities, Podkopaev and Ramdas (2021) study such statements under structured shift, and van der Laan and Alaa (2024) recalibrate the models with conformal guarantees. These works are the closest precedents to our calibration view. The key distinction is that PICPI produces the conditioning interval and the uncertainty statement simultaneously through the self-consistency condition, which is exactly the structure exploited in our theoretical results. 7. Conclusions In summary, we have introduced PICPIs as a framework for local, data-driven uncertainty quantification that bridges the gap between calibration and conformal prediction. By conditioning on the model’s predicted probabilities, PICPIs provide intervals that are both interpretable and adaptive to the precision of the underlying predictor. Unlike traditional calibration methods, which rely on fixed bins or post-hoc adjustments, our approach identifies regions of the prediction space where the model is self-certified to be reliable, yielding uncertainty statements that are directly actionable without modifying the base predictor. We show that, under suitable regularity conditions, the width of these intervals shrinks with increasing sample size, enabling increasingly precise risk characterization. Practically, PICPIs offer a tool for probabilistic prediction and multi-class classification with respective coverage guarantees. Overall, our results highlight that leveraging the model’s own risk stratification can produce rigorous, locally adaptive uncertainty quantification, providing an alternative to conventional recalibration or fully conditional prediction frameworks. This opens up directions on extending PICPIs to structured settings inherently require prediction-conditioned uncertainty quantification, such as representation-based inference (where the predictions acts as a low-dimensional summary), covariate shift, and distribution shift (where uncertainty guarantees may need to hold conditionally on prediction-induced strata that capture regions of varying model reliability). On the theoretical side, a direction for future work is to develop adaptive rates on the width bound that exploit additional properties or parametric structure of the prediction model. References Barber et al. (2021a) Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. The limits of distribution-free conditional predictive inference. Information and Inference: A Journal of the IMA, 10(2):455–482, 2021a. Barber et al. (2021b) 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, 2021b. Barber et al. (2022) Rina Foygel Barber, Emmanuel J. Candès, Aaditya Ramdas, and Ryan J. Tibshirani. Conformal prediction beyond exchangeability, 2022. Bars and Humbert (2025) Batiste Le Bars and Pierre Humbert. On volume minimization in conformal regression. In Proceedings Of The 42nd International Conference On Machine Learning, volume 267 of PMLR, pages 3030–3053. PMLR, 2025. Berta et al. (2024) Eugène Berta, Francis Bach, and Michael Jordan. Classifier calibration with ROC-regularized isotonic regression. In Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, volume 238 of Proceedings of Machine Learning Research, pages 1972–1980. PMLR, 02–04 May 2024. Blasiok and Nakkiran (2024) Jaroslaw Blasiok and Preetum Nakkiran. Smooth ece: Principled reliability diagrams via kernel smoothing. In International Conference on Learning Representations, volume 2024, pages 1679–1707, 2024. Carrell et al. (2022) A Michael Carrell, Neil Mallinar, James Lucas, and Preetum Nakkiran. The calibration generalization gap. arXiv preprint arXiv:2210.01964, 2022. Carriero et al. (2025) Alex Carriero, Kim Luijken, Anne de Hond, Karel GM Moons, Ben van Calster, and Maarten van Smeden. The harms of class imbalance corrections for machine learning based prediction models: a simulation study. Statistics in Medicine, 44(3-4):e10320, 2025. Chernozhukov et al. (2021) Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu. Distributional conformal prediction. Proceedings of the National Academy of Sciences, 118(48):e2107794118, 2021. Cox (1958) David Roxbee Cox. Some problems connected with statistical inference. The Annals of Mathematical Statistics, 29(2):357–372, 1958. Cox (1988) David Roxbee Cox. Some aspects of conditional and asymptotic inference: A review. Sankhyā: The Indian Journal of Statistics, Series A, pages 314–337, 1988. Dawid (1982) A Philip Dawid. The well-calibrated bayesian. Journal of the American Statistical Association, 77(379):605–610, 1982. Ding et al. (2023) Tiffany Ding, Anastasios N. Angelopoulos, Stephen Bates, Michael I. Jordan, and Ryan J. Tibshirani. Class-conditional conformal prediction with many classes. In Advances In Neural Information Processing Systems, volume 36, pages 64555–64576, 2023. doi: 10.52202/075280-2817. Duchi (2025) John C Duchi. A few observations on sample-conditional coverage in conformal prediction. arXiv preprint arXiv:2503.00220, 2025. Fisher (1934) Ronald Aylmer Fisher. Statistical methods for research workers. 1934. Foster and Vohra (1998) Dean P Foster and Rakesh V Vohra. Asymptotic calibration. Biometrika, 85(2):379–390, 1998. Gibbs and Candes (2021) Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems, 34:1660–1672, 2021. Gibbs et al. (2025) Isaac Gibbs, John J Cherian, and Emmanuel J Candès. Conformal prediction with conditional guarantees. Journal of the Royal Statistical Society Series B: Statistical Methodology, 87(4):1100–1126, 2025. Gneiting and Raftery (2007) Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477):359–378, 2007. Guo et al. (2017) Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger. On calibration of modern neural networks. In International Conference on Machine Learning, pages 1321–1330. PMLR, 2017. Gupta et al. (2020) Chirag Gupta, Aleksandr Podkopaev, and Aaditya Ramdas. Distribution-free binary classification: prediction sets, confidence intervals and calibration. Advances in Neural Information Processing Systems, 33:3711–3723, 2020. Hébert-Johnson et al. (2018) Ursula Hébert-Johnson, Michael Kim, Omer Reingold, and Guy Rothblum. Multicalibration: Calibration for the (computationally-identifiable) masses. In International Conference on Machine Learning, pages 1939–1948. PMLR, 2018. Izbicki et al. (2020) Rafael Izbicki, Gilson Shimizu, and Rafael Stern. Flexible distribution-free conditional predictive bands using density estimators. In Silvia Chiappa and Roberto Calandra, editors, Proceedings Of The Twenty Third International Conference On Artificial Intelligence And Statistics, volume 108 of PMLR, pages 3068–3077. PMLR, 2020. Izbicki et al. (2022) Rafael Izbicki, Gilson Shimizu, and Rafael B Stern. Cd-split and hpd-split: Efficient conformal regions in high dimensions. Journal of Machine Learning Research, 23(87):1–32, 2022. Jung et al. (2023) Christopher Jung, Georgy Noarov, Ramya Ramalingam, and Aaron Roth. Batch multivalid conformal prediction. In International Conference on Learning Representations, 2023. Kim et al. (2020) Byol Kim, Chen Xu, and Rina Barber. Predictive inference is free with the jackknife+-after-bootstrap. Advances in Neural Information Processing Systems, 33:4138–4149, 2020. Kleinberg et al. (2016) Jon Kleinberg, Sendhil Mullainathan, and Manish Raghavan. Inherent trade-offs in the fair determination of risk scores. arXiv preprint arXiv:1609.05807, 2016. Kull et al. (2017) Meelis Kull, Telmo Silva Filho, and Peter Flach. Beta calibration: a well-founded and easily implemented improvement on logistic calibration for binary classifiers. In Artificial Intelligence and Statistics, pages 623–631. PMLR, 2017. Kull et al. (2019) Meelis Kull, Miquel Perello Nieto, Markus Kängsepp, Telmo Silva Filho, Hao Song, and Peter Flach. Beyond temperature scaling: Obtaining well-calibrated multi-class probabilities with dirichlet calibration. Advances in Neural Information Processing Systems, 32, 2019. Lei and Wasserman (2014) 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. Lei et al. (2013) Jing Lei, James Robins, and Larry Wasserman. Distribution-free prediction sets. Journal of the American Statistical Association, 108(501):278–287, 2013. Lei et al. (2018) 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. Liang and Barber (2025) Ruiting Liang and Rina Foygel Barber. Algorithmic stability implies training-conditional coverage for distribution-free prediction methods. The Annals of Statistics, 53(4):1457–1482, 2025. Martinez Gil et al. (2024) Natalia Martinez Gil, Dhaval Patel, Chandra Reddy, Giri Ganapavarapu, Roman Vaculin, and Jayant Kalagnanam. Identifying homogeneous and interpretable groups for conformal prediction. In Proceedings of the Fortieth Conference on Uncertainty in Artificial Intelligence, volume 244 of PMLR, pages 2471–2485. PMLR, 2024. Naeini et al. (2015) Mahdi Pakdaman Naeini, Gregory Cooper, and Milos Hauskrecht. Obtaining well calibrated probabilities using bayesian binning. In Proceedings of the AAAI conference on Artificial Intelligence, volume 29, 2015. Niculescu-Mizil and Caruana (2005) Alexandru Niculescu-Mizil and Rich Caruana. Predicting good probabilities with supervised learning. In Proceedings of the 22nd International Conference on Machine Learning, pages 625–632, 2005. Papadopoulos et al. (2002) Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In Proceedings of the Thirteenth European Conference on Machine Learning, volume 2430 of Lecture Notes in Computer Science, pages 345–356. Springer, 2002. doi: 10.1007/3-540-36755-1_29. Platt et al. (1999) John Platt et al. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers, 10(3):61–74, 1999. Pleiss et al. (2017) Geoff Pleiss, Manish Raghavan, Felix Wu, Jon Kleinberg, and Kilian Q Weinberger. On fairness and calibration. Advances in Neural Information Processing Systems, 30, 2017. Podkopaev and Ramdas (2021) Aleksandr Podkopaev and Aaditya Ramdas. Distribution-free uncertainty quantification for classification under label shift. In Uncertainty in Artificial Intelligence, pages 844–853. PMLR, 2021. Reid (1995) Nancy Reid. The roles of conditioning in inference. Statistical Science, 10(2):138 – 157, 1995. doi: 10.1214/s/1177010027. URL https://doi.org/10.1214/s/1177010027. Romano et al. (2019) Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in Neural Information Processing Systems, 32, 2019. Shafer and Vovk (2008) Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008. Tibshirani et al. (2019) Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas. Conformal prediction under covariate shift. Advances in Neural Information Processing Systems, 32, 2019. Vaicenavicius et al. (2019) Juozas Vaicenavicius, David Widmann, Carl Andersson, Fredrik Lindsten, Jacob Roll, and Thomas Schön. Evaluating model calibration in classification. In Proceedings of the 22nd International Conference on Artificial Intelligence And Statistics, volume 89 of PMLR, pages 3459–3467. PMLR, 2019. van der Laan and Alaa (2024) Lars van der Laan and Ahmed M Alaa. Self-calibrating conformal prediction. Advances in Neural Information Processing Systems, 37:107138–107170, 2024. Vovk (2012) Vladimir Vovk. Conditional validity of inductive conformal predictors. In Asian Conference on Machine Learning, pages 475–490. PMLR, 2012. Vovk (2015) Vladimir Vovk. Cross-conformal predictors. Annals of Mathematics and Artificial Intelligence, 74(1):9–28, 2015. Vovk and Petej (2014) Vladimir Vovk and Ivan Petej. Venn–abers predictors. In Proceedings Of The Thirtieth Conference On Uncertainty In Artificial Intelligence, pages 829–838. AUAI Press, 2014. Vovk et al. (2003a) Vladimir Vovk, David Lindsay, Ilia Nouretdinov, Alex Gammerman, et al. Mondrian confidence machine. Technical Report, 2003a. Vovk et al. (2003b) Vladimir Vovk, Glenn Shafer, and Ilia Nouretdinov. Self-calibrating probability forecasting. Advances in Neural Information Processing Systems, 16, 2003b. Vovk et al. (2005) Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in a random world. Springer, 2005. Vovk et al. (2018) Vladimir Vovk, Ilia Nouretdinov, Valery Manokhin, and Alexander Gammerman. Cross-conformal predictive distributions. In Conformal and Probabilistic Prediction and Applications, pages 37–51. PMLR, 2018. Yang and Kuchibhotla (2025) Yachong Yang and Arun Kumar Kuchibhotla. Selection and aggregation of conformal prediction sets. Journal of the American Statistical Association, 120(549):435–447, 2025. Zadrozny and Elkan (2002) Bianca Zadrozny and Charles Elkan. Transforming classifier scores into accurate multiclass probability estimates. In Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 694–699, 2002. Appendix A Proofs A.1. Proof of Theorem 2.2 Proof. For each discretized interval [a,b][a,b] with empirical count N[a,b]=∑ip(xi)∈[a,b]>0N_[a,b]= _i 1\p(x_i)∈[a,b]\>0, we apply Hoeffding’s inequality and a union bound to derive that with probability at least 1−δ1-δ: |μ^[a,b]−μ[a,b]|≤2log(|ℐ|/δ)N[a,b], | μ_[a,b]- _[a,b]|≤ 2 (|I|/δ)N_[a,b], where μ[a,b]=ℙ(Y=1∣p(X)∈[a,b]) _[a,b]=P(Y=1 p(X)∈[a,b]), μ^[a,b]=∑i=1nyi⋅(p(xi)∈[a,b])∑i=1n(p(xi)∈[a,b]) μ_[a,b]= _i=1^ny_i· 1(p(x_i)∈[a,b]) _i=1^n 1(p(x_i)∈[a,b]) and |ℐ|=K2|I|=K^2 is the number of candidate intervals. On the event that all these bounds hold simultaneously, Line 6 implies that μ[a,b]∈[a,b]. _[a,b]∈[a,b]. This confirms that with probability at least 1−δ1-δ, ℐpI_p is a set of PICPIs. We note that a relaxed alternative for Line 6 is to add indicator functions at the boundary endpoints 0 and 1. That is, ∑i=1nyi⋅(p(xi)∈[a,b])∑i=1n(p(xi)∈[a,b]) _i=1^ny_i· 1(p(x_i)∈[a,b]) _i=1^n 1(p(x_i)∈[a,b]) ∈[a+2log(K2/δ)∑i=1n(p(xi)∈[a,b])a=0,b−2log(K2/δ)∑i=1n(p(xi)∈[a,b])b=1]. ∈ [a+2 (K^2/δ) _i=1^n 1(p(x_i)∈[a,b]) 1\a=0\,b-2 (K^2/δ) _i=1^n 1(p(x_i)∈[a,b]) 1\b=1\ ]. (21) This formulation may be preferable in edge cases where Y=0Y=0 or Y=1Y=1 almost surely. The proof still holds since μ[a,b]∈[0,1] _[a,b]∈[0,1]. For simplicity, we proceed with the simpler formulation in the main text. ∎ A.2. Proof of Theorem 3.2 Proof. Let ϵstat= _stat= 2(2log(3n/δ)λn)1/3, \ 2 ( 2 (3n/δ)λ n )^1/3, ϵ1= _1= ϵstat+4ϵδ1, \ _stat+ 4ε _1, K= K= max1,⌊14ϵ1⌋. \ \1, 14 _1 \. If K=1K=1, the interval [0,1][0,1] is returned by the algorithm and covers all prediction values, proving the result directly with C(n,ϵ,δ,δ1)=1C(n,ε,δ, _1)=1. For K≥2K≥ 2, the proof consists of three steps. Step 1. We show that there exists J∈ℤ+J _+ and disjoint intervals [aj,bj],j∈[J][a_j,b_j],j∈[J] such that: 1. For any j∈[J]j∈[J], aj+ϵstat⋅(aj≠0)≤ℙ[Y=1∣p(X)∈[aj,bj]]≤bj−ϵstat⋅(bj≠1); a_j+ _stat· 1(a_j≠ 0) [Y=1 p(X)∈[a_j,b_j]]≤ b_j- _stat· 1(b_j≠ 1); (22) 2. ℙ(p(X)∉∪j∈[J][aj,bj])≤δ1; (p(X)∉ _j∈[J][a_j,b_j] )≤ _1; (23) 3. For any j∈[J]j∈[J], 4ϵ1≤|bj−aj|≤56ϵ1⋅log(6/δ1). 4 _1≤ |b_j-a_j |≤ 56 _1· (6/ _1). (24) First, applying Lemma 3.4 with ν being the distribution of p(X)p(X), uncovered mass budget δ1 _1, and accuracy parameter ϵ1 _1, we see that there exist disjoint intervals [ai,bi],i∈[L][a_i,b_i],i∈[L] that satisfy ℙ(p(X)∉∪i∈[L][ai,bi])≤δ1/2, (p(X)∉ _i∈[L][a_i,b_i] )≤ _1/2, Eq. (24), and ai+ϵ1⋅(ai≠0)≤[p(X)∣p(X)∈[ai,bi]]≤bi−ϵ1⋅(bi≠1). a_i+ _1· 1(a_i≠ 0) [p(X) p(X)∈[a_i,b_i]]≤ b_i- _1· 1(b_i≠ 1). (25) By Lemma 3.5, there exist i1,…,iJ∈[L]i_1,…,i_J∈[L] such that for any j∈[J]j∈[J], |ℙ(Y=1∣p(X)∈[aij,bij])−[p(X)∣p(X)∈[aij,bij]]|≤4ϵ/δ1 |P(Y=1 p(X)∈[a_i_j,b_i_j])-E[p(X) p(X)∈[a_i_j,b_i_j]] |≤ 4ε/ _1 (26) and ℙ(p(X)∈∪j′∉i1,…,iJ[aj′,bj′])≤δ1/2. (p(X)∈ _j ∉\i_1,…,i_J\[a_j ,b_j ] )≤ _1/2. Combining Eq. (26) and Eq. (25), we have ℙ(Y=1∣p(X)∈[aij,bij]) P(Y=1 p(X)∈[a_i_j,b_i_j]) ≤[p(X)∣p(X)∈[aij,bij]]+4ϵ/δ1 [p(X) p(X)∈[a_i_j,b_i_j]]+4ε/ _1 ≤bij−ϵ1⋅(bij≠1)+4ϵ/δ1 ≤ b_i_j- _1· 1(b_i_j≠ 1)+4ε/ _1 =bij−ϵstat⋅(bij≠1), =b_i_j- _stat· 1(b_i_j≠ 1), where the last step follows from the definition of ϵ1 _1, and ℙ(Y=1∣p(X)∈[aij,bij]) P(Y=1 p(X)∈[a_i_j,b_i_j]) ≥[p(X)∣p(X)∈[aij,bij]]−4ϵ/δ1 [p(X) p(X)∈[a_i_j,b_i_j]]-4ε/ _1 ≥aij+ϵ1⋅(aij≠0)−4ϵ/δ1 ≥ a_i_j+ _1· 1(a_i_j≠ 0)-4ε/ _1 =aij+ϵstat⋅(aij≠0). =a_i_j+ _stat· 1(a_i_j≠ 0). (if aij=0a_i_j=0 then trivially aij≤ℙ(Y=1∣p(X)∈[aij,bij])a_i_j (Y=1 p(X)∈[a_i_j,b_i_j])) This establishes aij+ϵstat⋅(aij≠0)≤ℙ(Y=1∣p(X)∈[aij,bij])≤bij−ϵstat⋅(bij≠1),∀j∈[J]. a_i_j+ _stat· 1(a_i_j≠ 0) (Y=1 p(X)∈[a_i_j,b_i_j])≤ b_i_j- _stat· 1(b_i_j≠ 1), \ ∀ j∈[J]. Furthermore, ℙ(p(X)∉∪j∈[J][aij,bij])≤ (p(X)∉ _j∈[J][a_i_j,b_i_j] )≤ ℙ(p(X)∉∪i∈[L][ai,bi])+ℙ(p(X)∈∪j′∉i1,…,iJ[aj′,bj′]) \ P (p(X)∉ _i∈[L][a_i,b_i] )+P (p(X)∈ _j ∉\i_1,…,i_J\[a_j ,b_j ] ) ≤ ≤ δ1. \ _1. In sum, the intervals [ai1,bi1],…,[aiJ,biJ][a_i_1,b_i_1],…,[a_i_J,b_i_J] satisfy Eq. (22)-(24) (where we view [aij,bij][a_i_j,b_i_j] as [aj,bj][a_j,b_j]). Step 2. We show that with probability at least 1−δ1-δ over the randomness of (xi,yi)i=1n\(x_i,y_i)\_i=1^n, the intervals in Step 1 will be included in the set of intervals ℐI output by Algorithm 1. For simplicity of notation, we let [a1,b1],…,[aJ,bJ][a_1,b_1],…,[a_J,b_J] denote the intervals found in Step 1; i.e. [a1,b1],…,[aJ,bJ]=[ai1,bi1],…,[aiJ,biJ][a_1,b_1],…,[a_J,b_J]=[a_i_1,b_i_1],…,[a_i_J,b_i_J]. Let zi=p(xm+i)z_i=p(x_m+i) be a shorthand for the predicted probability. Define events A1= A_1= ∑i=1n(zi∈[aj,bj])≥2λϵ1n,∀j∈[J] \ \ _i=1^n 1(z_i∈[a_j,b_j])≥ 2λ _1n,∀ j∈[J] \ A2= A_2= |∑i=1nyi⋅(zi∈[aj,bj])∑i=1n(zi∈[aj,bj])−ℙ(Y=1∣p(X)∈[aj,bj])|≤ϵstat/2,∀j∈[J]. \ \ | _i=1^ny_i· 1(z_i∈[a_j,b_j]) _i=1^n 1(z_i∈[a_j,b_j])-P(Y=1 p(X)∈[a_j,b_j]) |≤ _stat/2,∀ j∈[J] \. Next, we show that under events A1,A2A_1,A_2, the intervals [a1,b1],…,[aJ,bJ][a_1,b_1],…,[a_J,b_J] are included in ℐI by Algorithm 1. Indeed, we have ∑i=1nyi⋅(zi∈[aj,bj])∑i=1n(zi∈[aj,bj])≤ _i=1^ny_i· 1(z_i∈[a_j,b_j]) _i=1^n 1(z_i∈[a_j,b_j])≤ ℙ(Y=1∣p(X)∈[aj,bj])+ϵstat/2 \ P(Y=1 p(X)∈[a_j,b_j])+ _stat/2 ≤ ≤ bj−ϵstat⋅(bj≠1)+ϵstat/2 \ b_j- _stat· 1(b_j≠ 1)+ _stat/2 ≤ ≤ bj−ϵstat/2⋅(bj≠1) \ b_j- _stat/2· 1(b_j≠ 1) ≤ ≤ bj−2log(2K2/δ)∑i=1n(zi∈[aj,bj]), \ b_j-2 (2K^2/δ) _i=1^n 1(z_i∈[a_j,b_j]), where the last inequality is due to 2log(2K2/δ)∑i=1n(zi∈[aj,bj])≤ 2 (2K^2/δ) _i=1^n 1(z_i∈[a_j,b_j])≤ 22log(3n/δ)2λϵ1n \ 2 2 (3n/δ)2λ _1n ≤ ≤ (2log(3n/δ)λn)1/3=ϵstat/2, \ ( 2 (3n/δ)λ n )^1/3= _stat/2, under event A1A_1. Similarly, ∑i=1nyi⋅(zi∈[aj,bj])∑i=1n(zi∈[aj,bj])≥aj+2log(2K2/δ)∑i=1n(zi∈[aj,bj]). _i=1^ny_i· 1(z_i∈[a_j,b_j]) _i=1^n 1(z_i∈[a_j,b_j])≥ a_j+2 (2K^2/δ) _i=1^n 1(z_i∈[a_j,b_j]). Therefore, the interval [aj,bj][a_j,b_j] will be included in ℐI under A1∩A2A_1∩ A_2. Next, we show ℙ(A1∩A2)≥1−δP(A_1∩ A_2)≥ 1-δ. Since |bj−aj|≥4ϵ1, |b_j-a_j |≥ 4 _1, applying Lemma 3.3 with c=4ϵ1c=4 _1, we have ℙ(A1)≥ (A_1)≥ 1−Jexp(−λϵ1n2) \ 1-J (- λ _1n2 ) ≥ ≥ 1−nexp(−λϵstatn2) \ 1-n (- λ _statn2 ) ≥ ≥ 1−nexp(−log(3n/δ)) \ 1-n (- (3n/δ) ) (n≥λ−1log(3n/δ)/2n≥λ^-1 (3n/δ)/ 2) ≥ ≥ 1−δ/3. \ 1-δ/3. Under A1A_1, by Hoeffding’s inequality and a union bound, we have ℙ(A2∣A1)≥ (A_2 A_1)≥ 1−2Jexp(−ϵstat2⋅∑i=1n(zi∈[aj,bj])/2) \ 1-2J (- _stat^2· _i=1^n 1(z_i∈[a_j,b_j])/2 ) ≥ ≥ 1−2nexp(−ϵstat2⋅λϵ1n) \ 1-2n (- _stat^2·λ _1n ) ≥ ≥ 1−2nexp(−ϵstat3⋅λn) \ 1-2n (- _stat^3·λ n ) ≥ ≥ 1−2nexp(−16log(3n/δ)λn⋅λn) \ 1-2n (- 16 (3n/δ)λ n·λ n ) ≥ ≥ 1−2nexp(−16log(3n/δ)) \ 1-2n (-16 (3n/δ) ) ≥ ≥ 1−2n(3n/δ)−16 \ 1-2n(3n/δ)^-16 ≥ ≥ 1−2δ16316n15 \ 1- 2δ^163^16n^15 ≥ ≥ 1−δ/6. \ 1-δ/6. Combining, we have ℙ(A1∩A2)≥1−δ. (A_1∩ A_2)≥ 1-δ. Step 3. We show that the intervals [a1,b1],…,[aJ,bJ][a_1,b_1],…,[a_J,b_J] in Step 2 satisfy the property ℙx∼μ(∃[aj,bj]s.t.|bj−aj|≤C(n,ϵ,δ,δ1)andp(x)∈[aj,bj])≥1−δ1, _x μ (∃[a_j,b_j] \ s.t. \ |b_j-a_j|≤ C(n,ε,δ, _1) \ and \ p(x)∈[a_j,b_j] )≥ 1- _1, where C(n,ϵ,δ,δ1) C(n,ε,δ, _1) =log(6/δ1)⋅56(2(2log(3n/δ)λn)1/3+4ϵδ1). = (6/ _1)· 56 (2 ( 2 (3n/δ)λ n )^1/3+ 4ε _1 ). Define event A3=∃[aj,bj]s.t.|b−a|≤C(n,ϵ,δ,δ1)andp(x)∈[aj,bj]. A_3= \∃[a_j,b_j] \ s.t. \ |b-a|≤ C(n,ε,δ, _1) \ and \ p(x)∈[a_j,b_j] \. By Eq. (24), for any j∈[J]j∈[J], |bj−aj|≤C(n,ϵ,δ,δ1). |b_j-a_j|≤ C(n,ε,δ, _1). It follows from Eq. (23) that ℙx∼μ(∃[aj,bj]s.t.|bj−aj|≤C(n,ϵ,δ,δ1)andp(x)∈[aj,bj])= _x μ (∃[a_j,b_j] \ s.t. \ |b_j-a_j|≤ C(n,ε,δ, _1) \ and \ p(x)∈[a_j,b_j] )= ℙ(p(X)∈∪j∈[J][aj,bj]) \ P (p(X)∈ _j∈[J][a_j,b_j] ) ≥ ≥ 1−δ1. \ 1- _1. ∎ A.3. Proof of Lemma 3.3 Proof. Fix an interval [aj,bj][a_j,b_j] for some j∈[J]j∈[J]. Define wi:=(zi∈[aj,bj])w_i:= 1(z_i∈[a_j,b_j]) for i∈[n]i∈[n]. Then wii=1n\w_i\_i=1^n are i.i.d. Bernoulli random variables with mean [wi]=ℙ(zi∈[aj,bj])≥λ(bj−aj)≥λc. [w_i]=P(z_i∈[a_j,b_j])≥λ(b_j-a_j)≥λ c. Let Sj:=∑i=1nwiS_j:= _i=1^nw_i. The total expected sample count satisfies [Sj]≥λcnE[S_j]≥λ cn. By the lower-tail multiplicative Chernoff bound, ℙ(Sj≤(1−δ)[Sj])≤exp(−δ2[Sj]/2).P(S_j≤(1-δ)E[S_j])≤ (-δ^2E[S_j]/2). Applying it with δ=1/2δ=1/2, ℙ(∑i=1nwi<λcn/2)≤ ( _i=1^nw_i<λ cn/2 )≤ ℙ(Sj≤[Sj]2) \ P (S_j≤ E[S_j]2 ) ≤ ≤ exp(−[Sj]8) \ (- E[S_j]8 ) ≤ ≤ exp(−λcn8). \ (- λ cn8 ). By a union bound over all j∈[J]j∈[J], with probability at least 1−Jexp(−λcn8)1-J (- λ cn8 ), ∑i=1n(zi∈[aj,bj])≥λcn/2 _i=1^n 1(z_i∈[a_j,b_j])≥λ cn/2 holds for all j∈[J]j∈[J]. This completes the proof. ∎ A.4. Proof of Lemma 3.4 Proof. We present a constructive proof based on discretization of [0,1][0,1]. To simplify notation, we abbreviate Z∼νE_Z ν as E and ℙZ∼νP_Z ν as ℙP. Step 1: Construct the graph. Consider [(i−1)/K,i/K],i=1,2,…,K[(i-1)/K,i/K],i=1,2,…,K with K=⌊1/(4ϵ)⌋K= 1/(4ε) . Without loss of generality assume K≥2K≥ 2. Construct a directed graph G=(V,E)G=(V,E) where V=[K]V=[K] and let the edges be defined as follows: 1→2,K−1←K1→ 2,K-1← K, where there is a directed edge i→ji→ j iff both of the following conditions hold: • |i−j|=1|i-j|=1, and • Either ℙ(Z∈[(i−1)/K,i/K]) and ([Z∣Z∈[(i−1)/K,i/K]]−i−0.5K)⋅(j−i)≥0, (Z∈[(i-1)/K,i/K]) and (E [Z Z∈[(i-1)/K,i/K] ]- i-0.5K )·(j-i)≥ 0, (27) or ℙ(Z∈[(i−1)/K,i/K])≤3⋅ℙ(Z∈[(j−1)/K,j/K]). (Z∈[(i-1)/K,i/K] )≤ 3·P (Z∈[(j-1)/K,j/K] ). (28) Eq. (27) can be interpreted as: if [Z∣Z∈[(i−1)/K,i/K]]E [Z Z∈[(i-1)/K,i/K] ] is closer to (i−1)/K(i-1)/K, then draw a directed edge i−1←i-1← i; if [Z∣Z∈[(i−1)/K,i/K]]E [Z Z∈[(i-1)/K,i/K] ] is closer to i/Ki/K, then draw a directed edge i→i+1i→ i+1. It directs an edge toward the side of the bin on which its conditional mean lies. Eq. (28) can be interpreted as: we draw a directed edge to i→ji→ j (for j=i−1j=i-1 or i+1i+1) if the probability mass on [(i−1)/K,i/K][(i-1)/K,i/K] is not substantially greater than that of [(j−1)/K,j/K][(j-1)/K,j/K]. It directs an edge toward a neighboring bin whose mass is not substantially smaller. By construction, it is obvious that every vertex has at least one outgoing edge since for positive-mass vertex Eq. (27) holds either for j=i−1j=i-1 or j=i+1j=i+1, and for zero-mass vertex Eq. (28) holds for both neighbors. There is at least one bidirectional adjacent pair. The existence of such pair follows from a continuity argument: since 1→2,K−1←K1→ 2,K-1← K, and every vertex has at least one outgoing edge, let k∗∈2,…,Kk^*∈\2,…,K\ be the smallest integer such that there exists an edge k∗−1←k∗k^*-1← k^*. Then k∗−1→k∗k^*-1→ k^* because every vertex has at least one outgoing edge, and if the outgoing edge for k∗−1k^*-1 is k∗−1→k∗−2k^*-1→ k^*-2, then k∗k^* would not be the smallest vertex as defined. This means that k∗−1⇆k∗k^*-1 k^* forms a bidirectional pair, thus confirming the existence. Now, denote the collection of H maximal blocks of consecutive vertices connected by bidirectional edges as Di:Di=li,li+1,…,ri, where li⇆li+1⇆…⇆ri,i∈[H] \D_i:D_i=\l_i,l_i+1,…,r_i\, where l_i l_i+1 … r_i, \ i∈[H]\ Then we have 1≤l1<r1<l2<r2<…<lH≤K. 1≤ l_1<r_1<l_2<r_2<…<l_H≤ K. Step 2: Construct intervals from the graph. We now construct intervals, where we first subdivide each block into short groups, and then assign the vertices between consecutive blocks to one of the two neighboring groups. We now divide each bidirectional block DiD_i into mi:=⌊(ri−li+1)/2⌋m_i:= (r_i-l_i+1)/2 consecutive groups, whereas every group has two vertices, except that the last has three when ri−li+1r_i-l_i+1 is odd. That is, denote the j-th group as αi,j,…,βi,j\ _i,j,…, _i,j\, with endpoints αi,j _i,j =li+2(j−1), =l_i+2(j-1), βi,j _i,j =li+2j−1,j<mi,ri,j=mi. = casesl_i+2j-1,&j<m_i,\\ r_i,&j=m_i. cases (29) Notice that for any i∈[H−1]i∈[H-1], all integer vertices between two consecutive bidirectional blocks must be connected separately by a leftwards chain and a rightwards chain (upon merging duplicated edges) as one move from left to right. That is, ri←⋯←wi→⋯→li+1, r_i←·s← w_i→·s→ l_i+1, for some wi∈ri,…,li+1w_i∈\r_i,…,l_i+1\. (Otherwise, there exists another bidirectional block between rir_i and li+1l_i+1.) We call wiw_i as an intervening source vertex, where we assign it the block on its left when constructing intervals. Denote the endpoints of constructed intervals as ui,vi,i∈[H]\u_i,v_i\,i∈[H]. Define the cut points and the ranges assigned to the blocks by ω0 _0 =0,ωH=K,ωi=minwi,li+1−1(i∈[H−1]), =0, _H=K, _i= \w_i,l_i+1-1\ (i∈[H-1]), ui u_i =ωi−1+1,vi=ωi(i∈[H]). = _i-1+1, v_i= _i (i∈[H]). (30) They satisfy ui≤li<ri≤vi,ui+1=vi+1(i∈[H−1]), u_i≤ l_i<r_i≤ v_i, u_i+1=v_i+1 (i∈[H-1]), and all edges in each assigned range have the form ui→⋯→li⇆⋯⇆ri←⋯←vi. u_i→·s→ l_i ·s r_i←·s← v_i. (31) Truncate the chains in Eq. (31) after at most ⌈log(3/δ)⌉ (3/δ) vertices by setting si=ui∨(li−⌈log(3/δ)⌉),ti=vi∧(ri+⌈log(3/δ)⌉). s_i=u_i (l_i- (3/δ) ), t_i=v_i (r_i+ (3/δ) ). (32) Attach the retained left chain to the first group and the retained right chain to the last group. That is, set si,j s_i,j =si,j=1,αi,j,j>1, = casess_i,&j=1,\\ _i,j,&j>1, cases ti,j t_i,j =βi,j,j<mi,ti,j=mi. = cases _i,j,&j<m_i,\\ t_i,&j=m_i. cases Define intervals Ii,j=[ai,j,bi,j]=[(si,j−1)/K,ti,j/K],i∈[H],j∈[mi]. I_i,j=[a_i,j,b_i,j]=[(s_i,j-1)/K,t_i,j/K], i∈[H],\ j∈[m_i]. (33) An example illustration of the intervals (in boxes, without truncation) as: ⋯←vi−1→ui→⋯→li⇆li+1⏟attach chain to first block pair⇆li+2⇆li+3⏟grouped by 2⇆⋯⇆ri−1⇆ri←⋯←vi⏟attach chain to last block pair→⋯. ·s← v_i-1→ u_i→·s→ l_i l_i+1_attach chain to first block pair l_i+2 l_i+3_grouped by 2 ·s r_i-1 r_i←·s← v_i_attach chain to last block pair→·s. (34) Truncation could be illustrated as: ⋯←vi−1→ui→⋯→li−⌈log(3/δ)⌉→⋯→li⇆li+1⏟attach truncated chain⇆li+2⇆li+3⏟grouped by 2⇆⋯. ·s← v_i-1→ u_i→·s→ l_i- (3/δ) →·s→ l_i l_i+1_attach truncated chain l_i+2 l_i+3_grouped by 2 ·s. (35) Remove any intervals of zero probability. Relabel those remaining as Ii=[ai,bi]I_i=[a_i,b_i], i∈[L]i∈[L], in increasing order. Step 3: Establish theoretical properties. We show that the intervals [ai,bi],i∈[L][a_i,b_i],i∈[L] satisfy the desired properties. By construction each interval contains between 11 and 2⌈log(3/δ)⌉+32 (3/δ) +3 vertices, thus Eq. (10) is satisfied. To show Eq. (8), Fix i∈[L]i∈[L]. By si≥ui,ti≤vis_i≥ u_i,t_i≤ v_i and the definition of ui,viu_i,v_i (after relabeling), there exist edges si→si+1s_i→ s_i+1 and ti−1←tit_i-1← t_i. We will show that si→si+1s_i→ s_i+1 implies that [Z∣Z∈[ai,bi]]≥ai+ϵ⋅(ai≠0)E[Z Z∈[a_i,b_i]]≥ a_i+ε· 1(a_i≠ 0) by considering the following three cases of si→si+1s_i→ s_i+1: Case 1: si=1s_i=1. Then ai=0a_i=0 and trivially [Z∣Z∈[ai,bi]]≥ai+ϵ⋅(ai≠0).E[Z Z∈[a_i,b_i]]≥ a_i+ε· 1(a_i≠ 0). Case 2: The edge is drawn by Eq. (27), which is [Z∣Z∈[(si−1)/K,si/K]]≥si−0.5K. [Z Z∈[(s_i-1)/K,s_i/K] ]≥ s_i-0.5K. Then we have [Z∣Z∈[ai,bi]]= [Z Z∈[a_i,b_i]]= ∑j=sitiℙ(Z∈[(j−1)/K,j/K])ℙ(Z∈[ai,bi])⋅[Z∣Z∈[(j−1)/K,j/K]] \ _j=s_i^t_i P(Z∈[(j-1)/K,j/K])P(Z∈[a_i,b_i])·E[Z Z∈[(j-1)/K,j/K]] ≥ ≥ [Z∣Z∈[(si−1)/K,si/K]] \ E [Z Z∈[(s_i-1)/K,s_i/K] ] ≥ ≥ si−0.5K \ s_i-0.5K ≥ ≥ ai+ϵ. \ a_i+ε. The first step comes from the law of total expectation and the fact that [(j−1)/K,j/K],j=si,…,ti[(j-1)/K,j/K],j=s_i,…,t_i forms a partition of [ai,bi][a_i,b_i]; the second step comes from [Z∣Z∈[(j−1)/K,j/K]]≥[Z∣Z∈[(si−1)/K,si/K]]E[Z Z∈[(j-1)/K,j/K]] [Z Z∈[(s_i-1)/K,s_i/K] ] for all j≥sij≥ s_i; the third step comes from Eq. (27); and the final step comes from ai=(si−1)/K,a_i=(s_i-1)/K, ϵ≤1/(4K)ε≤ 1/(4K), so ai+ϵ≤(si−3/4)/K≤(si−0.5)/Ka_i+ε≤(s_i-3/4)/K≤(s_i-0.5)/K. Case 3: The edge is drawn by Eq. (28), which is ℙ(Z∈[(si−1)/K,si/K])≤3⋅ℙ(Z∈[si/K,(si+1)/K]). (Z∈[(s_i-1)/K,s_i/K] )≤ 3·P (Z∈[s_i/K,(s_i+1)/K] ). Under Case 3, we have [Z∣Z∈[ai,bi]]= [Z Z∈[a_i,b_i]]= ∑j=sitiℙ(Z∈[(j−1)/K,j/K])ℙ(Z∈[ai,bi])⋅[Z∣Z∈[(j−1)/K,j/K]] \ _j=s_i^t_i P(Z∈[(j-1)/K,j/K])P(Z∈[a_i,b_i])·E[Z Z∈[(j-1)/K,j/K]] ≥ ≥ ℙ(Z∈[(si−1)/K,si/K])ℙ(Z∈[ai,bi])⋅[Z∣Z∈[(si−1)/K,si/K]] \ P(Z∈[(s_i-1)/K,s_i/K])P(Z∈[a_i,b_i])·E[Z Z∈[(s_i-1)/K,s_i/K]] +(1−ℙ(Z∈[(si−1)/K,si/K])ℙ(Z∈[ai,bi]))⋅[Z∣Z∈[si/K,(si+1)/K]] \ + (1- P(Z∈[(s_i-1)/K,s_i/K])P(Z∈[a_i,b_i]) )·E[Z Z∈[s_i/K,(s_i+1)/K]] ≥ ≥ 3/4⋅(si−1)/K+1/4⋅si/K \ 3/4·(s_i-1)/K+1/4· s_i/K ≥ ≥ ai+ϵ. \ a_i+ε. The first step comes from the law of total expectation and [(j−1)/K,j/K],j=si,…,ti[(j-1)/K,j/K],j=s_i,…,t_i forms a partition of [ai,bi][a_i,b_i]; the second step comes from pulling the first term out j=sij=s_i and then using [Z∣Z∈[(j−1)/K,j/K]]≥[Z∣Z∈[si/K,(si+1)/K]]E[Z Z∈[(j-1)/K,j/K]] [Z Z∈[s_i/K,(s_i+1)/K] ] for all j≥si+1j≥ s_i+1; the third step comes from Eq. (28); the last step comes from ai=(si−1)/K,a_i=(s_i-1)/K, ϵ≤1/(4K)ε≤ 1/(4K), so ai+ϵ≤(si−3/4)/Ka_i+ε≤(s_i-3/4)/K. Repeating similar arguments we can show ti−1←tit_i-1← t_i implies that [Z∣Z∈[ai,bi]]≤bi−ϵ⋅(bi≠1)E[Z Z∈[a_i,b_i]]≤ b_i-ε· 1(b_i≠ 1). Therefore, Eq. (8) is confirmed. Finally, we establish Eq. (9) by bounding the probability of the cells removed when the chains were truncated. Define qi:=ℙ(Z∈[(i−1)/K,i/K])q_i:=P(Z∈[(i-1)/K,i/K]). Recall that as demonstrated, on the left chain assigned to block DiD_i, the reverse edge is absent at every step. In particular, for each d∈1,…,li−uid∈\1,…,l_i-u_i\, the edge li−d+1→li−dl_i-d+1→ l_i-d does not exist and thus Eq.(28) is violated. Therefore, qli−d+1>3qli−d. q_l_i-d+1>3q_l_i-d. Iterating this inequality along the chain gives qli−d<3−dqli. q_l_i-d<3^-dq_l_i. We then have ∑d=⌈log(3/δ)⌉+1li−uiqli−d<∑d=⌈log(3/δ)⌉+1li−ui3−dqli. _d= (3/δ) +1^l_i-u_iq_l_i-d< _d= (3/δ) +1^l_i-u_i3^-dq_l_i. The analogous argument on the right chain (d∈1,…,vi−rid∈\1,…,v_i-r_i\) gives qri+d q_r_i+d <3−dqri <3^-dq_r_i ∑d=⌈log(3/δ)⌉+1vi−riqri+d _d= (3/δ) +1^v_i-r_iq_r_i+d <∑d=⌈log(3/δ)⌉+1vi−ri3−dqri. < _d= (3/δ) +1^v_i-r_i3^-dq_r_i. Combining, the total mass sum of discarded cells on two sides is at most ∑d=⌈log(3/δ)⌉+1li−uiqli−d+∑d=⌈log(3/δ)⌉+1vi−riqri+d _d= (3/δ) +1^l_i-u_iq_l_i-d+ _d= (3/δ) +1^v_i-r_iq_r_i+d ≤(qli+qri)∑d=⌈log(3/δ)⌉+1∞3−d ≤(q_l_i+q_r_i) _d= (3/δ) +1^∞3^-d =12 3−⌈log(3/δ)⌉(qli+qri). = 12\,3^- (3/δ) (q_l_i+q_r_i). (sum of geometric series) Recall that each point can at most belong to two of the intervals, ∑i=1H(qli+qri)≤2, _i=1^H(q_l_i+q_r_i)≤ 2, by union bound, ℙ(Z∉∪i∈[L][ai,bi]) (Z∉ _i∈[L][a_i,b_i] ) ≤12 3−⌈log(3/δ)⌉∑i=1H(qli+qri) ≤ 12\,3^- (3/δ) _i=1^H(q_l_i+q_r_i) ≤3−⌈log(3/δ)⌉≤e−⌈log(3/δ)⌉≤δ/3≤δ. ≤ 3^- (3/δ) ≤ e^- (3/δ) ≤δ/3≤δ. This proves Eq. (9) and completes the proof. ∎ A.5. Proof of Lemma 3.5 Proof. Let j1,j2,…,jMj_1,j_2,…,j_M be all integers in [L][L] such that Eq. (11) is not satisfied. If M=0M=0, then all intervals satisfy Eq. (11), and choosing i1,…,iJ=[L]\i_1,…,i_J\=[L] completes the proof. Otherwise, for any w∈j1,j2,…,jMw∈\j_1,j_2,…,j_M\, we have ℙ(V∈[aw,bw])>0P(V∈[a_w,b_w])>0 and 2ϵδ< 2εδ< |[U∣V∈[aw,bw]]−[V∣V∈[aw,bw]]| \ |E [U V∈[a_w,b_w] ]-E [V V∈[a_w,b_w] ] | ≤ ≤ [|U−V|∣V∈[aw,bw]] \ E [|U-V| V∈[a_w,b_w] ] (by the triangle inequality) = = 1ℙ(V∈[aw,bw])⋅[|U−V|⋅(V∈[aw,bw])]. \ 1P(V∈[a_w,b_w])·E [|U-V|· 1(V∈[a_w,b_w]) ]. Each point belongs to at most two intervals, so ∑w∈j1,…,jM(V∈[aw,bw])≤2 _w∈\j_1,…,j_M\ 1(V∈[a_w,b_w])≤ 2. It follows that 2ϵ≥ 2ε≥ 2[|U−V|] \ 2E[|U-V|] ≥ ≥ [|U−V|∑w∈j1,…,jM(V∈[aw,bw])] \ E [|U-V| _w∈\j_1,…,j_M\ 1(V∈[a_w,b_w]) ] = = ∑w∈j1,j2,…,jM[|U−V|⋅(V∈[aw,bw])] \ _w∈\j_1,j_2,…,j_M\E [|U-V|· 1(V∈[a_w,b_w]) ] > > ∑w∈j1,j2,…,jMℙ(V∈[aw,bw])⋅2ϵδ \ _w∈\j_1,j_2,…,j_M\P(V∈[a_w,b_w])· 2εδ ≥ ≥ ℙ(V∈∪w∈j1,j2,…,jM[aw,bw])⋅2ϵδ. \ P (V∈ _w∈\j_1,j_2,…,j_M\[a_w,b_w] )· 2εδ. Dividing both sides by 2ϵ/δ2ε/δ yields ℙ(V∈∪w∈j1,j2,…,jM[aw,bw])<δ. (V∈ _w∈\j_1,j_2,…,j_M\[a_w,b_w] )<δ. Setting i1,…,iJ:=[L]∖j1,…,jM\i_1,…,i_J\:=[L] \j_1,…,j_M\, the remaining intervals satisfy Eq. (11) by construction and Eq. (12). ∎ A.6. Proof of Theorem 4.1 Proof. Let Aj=p(X)∈IjA_j=\p(X)∈ I_j\, j∈[J]j∈[J], and A0=p(X)∉∪j=1JIjA_0=\p(X)∉ _j=1^JI_j\. These events form a partition. For each j with mj>0m_j>0, Eq. (13) gives Ctp(X,ϵ)=Eϵ(Ij)C_tp(X;ε)=E^ε(I_j) on AjA_j. The law of total probability therefore implies ℙ(p∗(X)∉Ctp(X;ϵ))≤∑j:mj>0mjℙ(p∗(X)∉Eϵ(Ij)∣Aj)+ℙ(A0). (p^*(X)∉ C_tp(X;ε) ) \ ≤ _j:m_j>0m_jP (p^*(X)∉ E^ε(I_j) A_j )+P(A_0). (36) For a positive-mass interval, let μj=[p∗(X)∣Aj] _j=E[p^*(X) A_j]. By the tower property and the PICPI condition, μj=[[Y∣X]∣Aj]=[Y∣Aj]∈Ij. _j=E[E[Y X] A_j]=E[Y A_j]∈ I_j. If p∗(X)∉Eϵ(Ij)p^*(X)∉ E^ε(I_j), its distance from every point of IjI_j, including μj _j, exceeds ϵε. Chebyshev’s inequality conditionally on AjA_j gives ℙ(p∗(X)∉Eϵ(Ij)∣Aj) (p^*(X)∉ E^ε(I_j) A_j ) ≤ℙ(|p∗(X)−μj|>ϵ∣Aj) (|p^*(X)- _j|>ε A_j ) ≤Var(p∗(X)∣Aj)ϵ2. ≤ (p^*(X) A_j)ε^2. Since ℙ(A0)=1−∑j=1JmjP(A_0)=1- _j=1^Jm_j, substituting them into Eq. (36) yields ℙ(p∗(X)∉Ctp(X;ϵ))≤ϵ−2∑j:mj>0mjVar(p∗(X)∣p(X)∈Ij)+(1−∑j=1Jmj). (p^*(X)∉ C_tp(X;ε) ) \ ≤ε^-2 _j:m_j>0m_j (p^*(X) p(X)∈ I_j )+ (1- _j=1^Jm_j ). This proves the theorem. ∎ A.7. Proof of Theorem 4.3 Proof. Fix a class g∈[G]g∈[G] and abbreviate π:=πg,t:=tg,J:=Jg,S:=Sg,π:= _g, t:=t_g, J:=J_g, S:=S_g, mj:=mg,j,aj:=ag,j,bj:=bg,j,j=1,…,J.m_j:=m_g,j, a_j:=a_g,j, b_j:=b_g,j, j=1,…,J. Define the disjoint events A0:=pg(X)∉⋃j=1J[aj,bj],Aj:=pg(X)∈[aj,bj],j=1,…,J.A_0:= \p_g(X)∉ _j=1^J[a_j,b_j] \, A_j:= \p_g(X)∈[a_j,b_j] \, j=1,…,J. Let μ0:=ℙ(Y(X)=g∣A0)∈[0,1],μj:=ℙ(Y(X)=g∣Aj)∈[aj,bj],j=1,…,J, _0:=P (Y(X)=g A_0 )∈[0,1], _j:=P (Y(X)=g A_j )∈[a_j,b_j], j=1,…,J, where the inclusion μj∈[aj,bj] _j∈[a_j,b_j] follows from the definition of PICPIs. Since the events A0,A1,…,AJA_0,A_1,…,A_J partition the sample space, π=ℙ(Y(X)=g)=tμ0+∑j=1Jmjμj. π=P (Y(X)=g )=t _0+ _j=1^Jm_j _j. (37) Now, ℙ(Y(X)∉Cmc(X,S),Y(X)=g) \ P (Y(X)∉ C_mc(X;S),\,Y(X)=g ) = = tμ0+∑j∉Smjμj, \ t _0+ _j∉ Sm_j _j, (38) because class g is omitted precisely on the uncovered region A0A_0 and on those intervals AjA_j with j∉Sj∉ S. Dividing Eq. (38) by π gives ℙ(Y(X)∉Cmc(X,S)∣Y(X)=g)=tμ0+∑j∉Smjμjπ. (Y(X)∉ C_mc(X;S) Y(X)=g )= t _0+ _j∉ Sm_j _jπ. (39) We first prove the upper-endpoint bound. Since μ0≤1 _0≤ 1 and μj≤bj _j≤ b_j for every j∉Sj∉ S, tμ0+∑j∉Smjμj≤t+∑j∉Smjbj.t _0+ _j∉ Sm_j _j≤ t+ _j∉ Sm_jb_j. Combining with Eq. (39) yields ℙ(Y(X)∉Cmc(X,S)∣Y(X)=g)≤Γgup(Sg).P (Y(X)∉ C_mc(X;S) Y(X)=g )≤ _g^up(S_g). Next we prove the lower-endpoint bound. By Eqs. (37) and (38), tμ0+∑j∉Smjμj=π−∑j∈Smjμj.t _0+ _j∉ Sm_j _j=π- _j∈ Sm_j _j. Since μj≥aj _j≥ a_j for all j∈Sj∈ S, tμ0+∑j∉Smjμj≤π−∑j∈Smjaj.t _0+ _j∉ Sm_j _j≤π- _j∈ Sm_ja_j. Combining again with Eq. (39) gives ℙ(Y(X)∉Cmc(X,S)∣Y(X)=g)≤Γglow(Sg).P (Y(X)∉ C_mc(X;S) Y(X)=g )≤ _g^low(S_g). Therefore ℙ(Y(X)∉Cmc(X,S)∣Y(X)=g)≤Γg(Sg).P (Y(X)∉ C_mc(X;S) Y(X)=g )≤ _g(S_g). It remains to prove exactness. Define U:=t+∑j∉Smjbj,L:=π−∑j∈Smjaj.U:=t+ _j∉ Sm_jb_j, L:=π- _j∈ Sm_ja_j. We show that the miscoverage numerator in Eq. (38) can attain minU,L \U,L\. Case 1: U≤LU≤ L. Set μ0=1,μj=bj for all j∉S. _0=1, _j=b_j\ for all j∉ S. Then the miscoverage numerator equals U. To make Eq. (37) hold, we need to choose μj∈[aj,bj] _j∈[a_j,b_j] for j∈Sj∈ S so that ∑j∈Smjμj=π−U. _j∈ Sm_j _j=π-U. This is possible because U≤L=π−∑j∈Smjaj⟹π−U≥∑j∈Smjaj,U≤ L=π- _j∈ Sm_ja_j π-U≥ _j∈ Sm_ja_j, and, since any admissible law must satisfy π≤t+∑j∉Smjbj+∑j∈Smjbj=U+∑j∈Smjbj,π≤ t+ _j∉ Sm_jb_j+ _j∈ Sm_jb_j=U+ _j∈ Sm_jb_j, we also have π−U≤∑j∈Smjbj.π-U≤ _j∈ Sm_jb_j. Hence π−Uπ-U lies in the interval [∑j∈Smjaj,∑j∈Smjbj], [ _j∈ Sm_ja_j,\ _j∈ Sm_jb_j ], so such a choice of μjj∈S\ _j\_j∈ S exists. Therefore equality is attained with value Uπ=Γgup(Sg)=Γg(Sg). Uπ= _g^up(S_g)= _g(S_g). Case 2: L<UL<U. Set μj=aj for all j∈S. _j=a_j\ for all j∈ S. Then the retained intervals contribute exactly ∑j∈Smjaj _j∈ Sm_ja_j, so by Eq. (37) the miscoverage numerator must equal L. It remains to show that one can choose μ0∈[0,1] _0∈[0,1] and μj∈[aj,bj] _j∈[a_j,b_j] for j∉Sj∉ S so that tμ0+∑j∉Smjμj=L.t _0+ _j∉ Sm_j _j=L. The attainable range of the left-hand side is [∑j∉Smjaj,t+∑j∉Smjbj]. [ _j∉ Sm_ja_j,\ t+ _j∉ Sm_jb_j ]. Because any admissible law must satisfy π≥∑j=1Jmjaj,π≥ _j=1^Jm_ja_j, we have L=π−∑j∈Smjaj≥∑j∉Smjaj.L=π- _j∈ Sm_ja_j≥ _j∉ Sm_ja_j. And by assumption, L<U=t+∑j∉Smjbj.L<U=t+ _j∉ Sm_jb_j. Hence L lies in the attainable range, so equality is attained with value Lπ=Γglow(Sg)=Γg(Sg). Lπ= _g^low(S_g)= _g(S_g). This proves exactness in both cases. ∎ A.8. Proof of Theorem 4.4 Proof. By construction, Γg(Sg⋆(αg))≤αg _g (S_g ( _g) )≤ _g for every g∈[G]g∈[G]. Applying Theorem 4.3 with S=S⋆(α)S=S (α) yields ℙ(Y(X)∉Cmc⋆(X,α)∣Y(X)=g)≤αg,P (Y(X)∉ C_mc (X;α) Y(X)=g )≤ _g, which proves the coverage claim. We now prove size-optimality. For any S=(S1,…,SG)S=(S_1,…,S_G), [|Cmc(X,S)|]=∑g=1Gℙ(g∈Cmc(X,S)). [|C_mc(X;S)| ]= _g=1^GP (g∈ C_mc(X;S) ). Since the intervals in ℐpgI_p_g are disjoint for each fixed class g, ℙ(g∈Cmc(X,S))=∑j∈Sgℙ(pg(X)∈[ag,j,bg,j])=∑j∈Sgmg,j.P (g∈ C_mc(X;S) )= _j∈ S_gP (p_g(X)∈[a_g,j,b_g,j] )= _j∈ S_gm_g,j. Therefore [|Cmc(X,S)|]=∑g=1G∑j∈Sgmg,j. [|C_mc(X;S)| ]= _g=1^G _j∈ S_gm_g,j. (40) The objective (40) separates over classes, and the feasibility constraints Γg(Sg)≤αg,g∈[G], _g(S_g)≤ _g, g∈[G], also separate over classes. Hence the global minimization problem decomposes into G independent classwise problems, and each classwise optimizer is exactly Sg⋆(αg)S_g ( _g) by (18). Summing the classwise minima proves the result. ∎ A.9. Proof of Theorem 4.5 Proof. Let ℰPIE_PI denote the event that all intervals returned by Algorithm 1 are genuine disjoint PICPIs simultaneously. By assumption, ℙ(ℰPI)≥1−δPI.P(E_PI)≥ 1- _PI. For each class g∈[G]g∈[G] and each candidate interval I∈pgI _p_g, define Zk(g,I):=(pg(Xk)∈I),k=1,…,ncal.Z_k^(g,I):= 1 (p_g(X_k)∈ I ), k=1,…,n_cal. Conditioned on the data-generating distribution, these are i.i.d. Bernoulli random variables with mean mg(I):=ℙ(pg(X)∈I).m_g(I):=P (p_g(X)∈ I ). Hence Hoeffding’s inequality gives, for every fixed (g,I)(g,I), ℙ(|1ncal∑k=1ncalZk(g,I)−mg(I)|>εm)≤2e−2ncalεm2.P ( | 1n_cal _k=1^n_calZ_k^(g,I)-m_g(I) |> _m )≤ 2e^-2n_cal _m^2. Applying a union bound over all candidate intervals in all classes yields ℙ(supg∈[G]supI∈pg|1ncal∑k=1ncalZk(g,I)−mg(I)|>εm)≤2Ncande−2ncalεm2=δm.P ( _g∈[G] _I _p_g | 1n_cal _k=1^n_calZ_k^(g,I)-m_g(I) |> _m )≤ 2N_cande^-2n_cal _m^2= _m. Therefore, with probability at least 1−δm1- _m, we have simultaneously for every class g and every candidate interval I∈pgI _p_g, |m^g(I)−mg(I)|≤εm. | m_g(I)-m_g(I) |≤ _m. Since each selected interval I^g,j I_g,j belongs to pgA_p_g by construction, it follows that mg,j:=ℙ(pg(X)∈I^g,j)≥m¯g,jfor all g and j.m_g,j:=P (p_g(X)∈ I_g,j )≥ m_g,j all g and j. Let ℰmE_m denote this event. Next, for each class g, the random variables (Yk=g) 1(Y_k=g) are i.i.d. Bernoulli with mean πg _g. Another Hoeffding bound and a union bound over g∈[G]g∈[G] give ℙ(supg∈[G]|π^g−πg|>επ)≤2Ge−2ncalεπ2=δπ.P ( _g∈[G] | π_g- _g |> _π )≤ 2G\,e^-2n_cal _π^2= _π. Hence, with probability at least 1−δπ1- _π, we simultaneously have π¯g≤πg≤π¯gfor all g∈[G]. π_g≤ _g≤ π_g all g∈[G]. Let ℰπE_π denote this event. By the union bound, ℙ(ℰPI∩ℰm∩ℰπ)≥1−δPI−δm−δπ.P (E_PI _m _π )≥ 1- _PI- _m- _π. Fix a realization of DcalD_cal in the event ℰ:=ℰPI∩ℰm∩ℰπ.E:=E_PI _m _π. Conditioned on this realization, the intervals I^g,j=[a^g,j,b^g,j] I_g,j=[ a_g,j, b_g,j] are now deterministic and genuine PICPIs, so Theorem 4.3 applies to them. Fix a class g∈[G]g∈[G] and a subset Sg⊆[J^g]S_g [ J_g]. By Theorem 4.3, ℙ(Y(X)∉C^mc(X,S)|Y(X)=g,Dcal) \ P (Y(X)∉ C_mc(X;S)\ |\ Y(X)=g,\ D_cal ) ≤ ≤ mint^g+∑j∉Sgm^g,jpopb^g,jπg, 1−∑j∈Sgm^g,jpopa^g,jπg, \ \ t_g+ _j∉ S_g m_g,j^pop b_g,j _g,\ 1- _j∈ S_g m_g,j^pop a_g,j _g \, (41) where m^g,jpop:=ℙ(pg(X)∈I^g,j∣Dcal),t^g:=1−∑j=1J^gm^g,jpop. m_g,j^pop:=P (p_g(X)∈ I_g,j D_cal ), t_g:=1- _j=1 J_g m_g,j^pop. Using t^g+∑j∉Sgm^g,jpopb^g,j=1−∑j∈Sgm^g,jpop−∑j∉Sg(1−b^g,j)m^g,jpop, t_g+ _j∉ S_g m_g,j^pop b_g,j=1- _j∈ S_g m_g,j^pop- _j∉ S_g(1- b_g,j) m_g,j^pop, together with m^g,jpop≥m¯g,j m_g,j^pop≥ m_g,j on ℰmE_m, we obtain t^g+∑j∉Sgm^g,jpopb^g,j≤1−∑j∈Sgm¯g,j−∑j∉Sg(1−b^g,j)m¯g,j. t_g+ _j∉ S_g m_g,j^pop b_g,j≤ 1- _j∈ S_g m_g,j- _j∉ S_g(1- b_g,j) m_g,j. Since πg≥π¯g _g≥ π_g on ℰπE_π, it follows that t^g+∑j∉Sgm^g,jpopb^g,jπg≤Γ^gup(Sg). t_g+ _j∉ S_g m_g,j^pop b_g,j _g≤ _g^up(S_g). Likewise, on ℰmE_m and ℰπE_π, ∑j∈Sgm^g,jpopa^g,j≥∑j∈Sgm¯g,ja^g,j,πg≤π¯g, _j∈ S_g m_g,j^pop a_g,j≥ _j∈ S_g m_g,j a_g,j, _g≤ π_g, so 1−∑j∈Sgm^g,jpopa^g,jπg≤1−∑j∈Sgm¯g,ja^g,jπ¯g=Γ^glow(Sg).1- _j∈ S_g m_g,j^pop a_g,j _g≤ 1- _j∈ S_g m_g,j a_g,j π_g= _g^low(S_g). Combining these inequalities with Ineq. (41) yields ℙ(Y(X)∉C^mc(X,S)|Y(X)=g,Dcal)≤Γ^g(Sg).P (Y(X)∉ C_mc(X;S)\ |\ Y(X)=g,\ D_cal )≤ _g(S_g). Since the argument holds simultaneously for every class and every subset on the event ℰE, the theorem follows. Finally, if S^g⋆(αg) S_g ( _g) is feasible for every class, then Γ^g(S^g⋆(αg))≤αg. _g ( S_g ( _g) )≤ _g. Applying the displayed bound above with Sg=S^g⋆(αg)S_g= S_g ( _g) gives ℙ(Y(X)∉C^mc⋆(X,α)|Y(X)=g,Dcal)≤αgP (Y(X)∉ C_mc (X;α)\ |\ Y(X)=g,\ D_cal )≤ _g for every class g∈[G]g∈[G]. ∎ Appendix B Implementation Details B.1. Baselines 1. PICPI. We construct PICPIs from the calibration split. For Task 1, We adapt the procedure of Eq. (13) for an ϵε-expansion of PICPIs. For Task 2, we use Eq. (14). That is, include class if its predicted probability lies in the an acceptance region by the subset solver in Algorithm 2. 2. Simultaneous confidence interval. Suppose a point estimator pθp_θ is parametrized by θ∈ℝdθ ^d. Given covariate value x, denote pθ∗(x)p_θ^*(x) as the true conditional probability when correctly specified. The simultaneous confidence interval for pθ∗(x)p_θ^*(x) is computed from a Wald confidence region for θ∗θ^* (using the chi-square critical value). Let θ θ be the fitted parameter and let V V be an estimate of its covariance matrix. Let pθ(x)=g−1(z(x)⊤θ),p_θ(x)=g^-1\! (z(x) θ ), where z(x)z(x) is the feature vector (for example, z(x)=(1x⊤)⊤z(x)=(1 \ x ) ) and g−1g^-1 is the inverse link function. The simultaneous confidence interval for pθ∗(x)p_θ^*(x) is SCI1−α(x)=[infθ∈1−αpθ(x),supθ∈1−αpθ(x)],SCI_1-α(x)= [ _θ _1-αp_θ(x),\; _θ _1-αp_θ(x) ], where 1−α=θ∈ℝd:(θ−θ^)⊤V^−1(θ−θ^)≤χd, 1−α2.C_1-α= \θ ^d:(θ- θ) V^-1(θ- θ)≤χ^2_d,\,1-α \. 3. Split conformal prediction. For each test covariate input x, the split conformal interval is C^1−α(x)=[p^(x)−q^1−α,p^(x)+q^1−α]∩[0,1], C_1-α(x)= [ p(x)- q_1-α,\ p(x)+ q_1-α ]∩[0,1], where q^1−α q_1-α is an empirical 1−α1-α quantile of the residual scores |Yi−p^(Xi)||Y_i- p(X_i)| on ncaln_cal calibration data points. The coverage satisfies ℙ(Y∈C^1−α(x)∣Xi,Yii∈[ncal])∈[1−α,1−α+(1+ncal)−1) (Y∈ C_1-α(x) \X_i,Y_i\_i∈[n_cal] )∈ [1-α,1-α+(1+n_cal)^-1 ) (42) (Papadopoulos et al., 2002). 4. Calibration-based interval. We discretize [0,1][0,1] into bins Bk=(bk−1,bk]B_k=(b_k-1,b_k] with midpoint mk=(bk−1+bk)/2m_k=(b_k-1+b_k)/2. We retain a subset of bins such that their mass-weighted absolute calibration error is close to that of PICPI; if that bin is not retained, we use the retained bin whose midpoint is closest to p^(x) p(x). 5. Fixed-width binning. We discretize [0,1][0,1] into fixed bins and use the bin containing p^(x) p(x) as the reported interval. 6. Conformal classification. This an additional baseline only for Task 2 (which differs from the others by not relying on intervals of point predictions). Let p^g(x) p_g(x) denote the predicted probability of class g. Given a test covariate input x, the label set includes classes whose predicted probability exceeds a threshold τ^α τ_α, Γ^1−α(x)=g∈[K]:p^g(x)≥τ^α, _1-α(x)=\g∈[K]: p_g(x)≥ τ_α\, where τ^α τ_α is obtained from the lower tail of scores on calibration data: τ^α=⌈α(ncal+1)⌉ smallest of p^Yi(Xi)i∈[ncal]. τ_α= α(n_cal+1) smallest of $\ p_Y_i(X_i)\_i∈[n_cal]$. Analogous to Eq. (42), it is marginally guaranteed that on the test data ℙ(Y∈Γ^1−α(x)∣Xi,Yii∈[ncal])∈[1−α,1−α+(1+ncal)−1).P (Y∈ _1-α(x) \X_i,Y_i\_i∈[n_cal] )∈ [1-α,1-α+(1+n_cal)^-1 ). Interval subset solver for Task 2 For PICPI, fixed-width binning, calibration-based intervals, and split conformal prediction, the classwise candidate intervals from the calibration set are passed to Algorithm 2. For each class the solver keeps a min-mass subset whose plugin two-sided certificate meets the target α. Simultaneous confidence intervals and conformal classification bypass the solver. For simultaneous CIs we include class g if the upper endpoint satisfies Ug(x)≥τU_g(x)≥τ, with τ the ⌈α(ncal+1)⌉ α(n_cal+1) -th smallest of UYi(Xi)\U_Y_i(X_i)\. For conformal classification, we use the calibration set to select a threshold on point predictions. B.2. Metrics For Task 1, we report the average interval length, calibration errors, the coverage of the prediction space [0,1][0,1], and the coverage of the target probability p∗(x)p^*(x) on the test data. For a method that outputs intervals C^(Xi)=[Li,Ui] C(X_i)=[L_i,U_i] on evaluation points Xii=1neval\X_i\_i=1^n_eval, let ℐ=I1,…,IMI=\I_1,…,I_M\ denote the set of distinct reported intervals, where Ij=[ℓj,uj]I_j=[ _j,u_j]. For each interval IjI_j, define its midpoint as mj=(ℓj+uj)/2m_j=( _j+u_j)/2 and let Aj=i:C^(Xi)=IjA_j=\i: C(X_i)=I_j\ be the set of evaluation points assigned to that interval. We estimate the expected calibration error of the intervals relative to their midpoints by: ECEmid=1M∑j=1M|1|Aj|∑i∈Ajp∗(Xi)−mj|.ECE_mid= 1M _j=1^M | 1|A_j| _i∈ A_jp^*(X_i)-m_j |. We also estimate the expected calibration error of the intervals relative to the predictions p(Xi)p(X_i) by: ECEpred=1M∑j=1M|1|Aj|∑i∈Ajp∗(Xi)−1|Aj|∑i∈Ajp(Xi)|.ECE_pred= 1M _j=1^M | 1|A_j| _i∈ A_jp^*(X_i)- 1|A_j| _i∈ A_jp(X_i) |. Intuitively, ECEmidECE_mid assesses the discrepancy between an interval’s midpoint and the average target probability p∗(X)p^*(X) among points assigned to that interval, capturing how well the interval is centered around the target probabilities. ECEpredECE_pred measures the discrepancy between the average prediction p(X)p(X) and the average target probability p∗(X)p^*(X) among points assigned to each interval, representing if grouping align with the underlying predictions. The reported relative ECE is the mean ECE of each method divided by the smallest mean ECE across all methods. For Task 2, we measure class-conditional error on label set C(X,τ)C(X;τ), ∑i=1nevalqg(Xi)g∉C(Xi;τ)∑i=1nevalqg(Xi). _i=1^n_evalq_g(X_i)1\g∉ C(X_i;τ)\ _i=1^n_evalq_g(X_i). We also report the average prediction-set size. B.3. Setup for Task 1 Univariate interval visualization (Figure 4). The covariate is sampled uniformly from a discrete grid X_GRID=−2,…,2X\_GRID=\-2,…,2\ with 10001000 equally spaced points, and the noise-free target probability used for evaluation is p∗(x)=σ(2x),p^*(x)=σ(2x), where σ(t)=1/(1+e−t)σ(t)=1/(1+e^-t). Observed labels are generated from the noisy logistic model Y∣X=x,ξ∼Bernoulli(σ(2x+ξ)),ξ∼N(0,0.12).Y X=x,ξ \! (σ(2x+ξ) ), ξ N(0,0.1^2). Thus the marginal conditional label probability is ℙ(Y=1∣X=x)=ξ[σ(2x+ξ)]P(Y=1 X=x)=E_ξ[σ(2x+ξ)], whereas the reported target-probability coverage uses the noise-free p∗(x)=σ(2x)p^*(x)=σ(2x) above. We set ntrain=2000,ncal=2000,neval=10000.n_train=2000, n_cal=2000, n_eval=10000. The common method settings are α=0.1α=0.1, PICPI with num_bin=50=50 in empirical mode, PICPI expansion parameter ϵ=0.02ε=0.02, and 5050 bins for the calibration-based and fixed-width baselines. The shared point predictor is scikit-learn logistic regression with its default ℓ2 _2 penalty and regularization strength 11. Misspecified decision-tree illustration (Figure 5). For the qualitative misspecification figure, the covariate is drawn from X∼Unif[−2,2],X [-2,2], with true probability p∗(x)=σ(1.2sin(2.5x)+0.8x),p^*(x)=σ\! (1.2 (2.5x)+0.8x ), and labels Y|X=x∼Bernoulli(p∗(x)).Y X=x (p^*(x)). The point predictor is a deliberately misspecified decision tree with max_depth = 5 and min_samples_leaf = 25, fitted on ntrain=20000,ncal=20000n_train=20000, n_cal=20000 samples. PICPI is built with num_bin=100=100 candidate grid points in empirical mode and expanded by ϵ=0.02ε=0.02. The curves are plotted on a dense evaluation grid of 600600 evenly spaced points in [−2,2][-2,2]. Multivariate result table (Table 1). The covariate is X∼N(0,I20),X N(0,I_20), and the true target probability is p∗(x)=σ(β0+x⊤β),β0=0,β=(1.6,−1.1, 0.9, 0,…,0)⊤∈ℝ20.p^*(x)=σ( _0+x β), _0=0, β=(1.6,\,-1.1,\,0.9,\,0,…,0) ^20. Labels are again generated from a noisy logistic model, Y∣X=x,ξ∼Bernoulli(σ(β0+x⊤β+ξ)),ξ∼N(0,0.12).Y X=x,ξ \! (σ( _0+x β+ξ) ), ξ N(0,0.1^2). As in the univariate experiment, the reported target-probability coverage uses the noise-free p∗(x)p^*(x) above, not the Gaussian-averaged marginal conditional label probability. Each replication uses ntrain=2000,ncal=2000,neval=10000,n_train=2000, n_cal=2000, n_eval=10000, and the table is averaged over 100100 Monte Carlo replications. The shared point predictor is scikit-learn logistic regression with its default ℓ2 _2 penalty. Unless otherwise stated, the method hyperparameters are the same as in the univariate setting: α=0.1α=0.1, PICPI with num_bin=50=50 candidate grid points and εPICPI=0.02 _PICPI=0.02, and 5050 bins for the two binning baselines. Split conformal implementation. For the Task 1 split conformal baseline, the implementation computes residual scores Ri=|Yi−p^(Xi)|R_i=|Y_i- p(X_i)| on the calibration split and sets u=min⌈(1−α)(ncal+1)⌉ncal,1,q^1−α=Quantileulinear(Rii∈[ncal]).u= \! \ (1-α)(n_cal+1) n_cal,1 \, q_1-α=Quantile^linear_u (\R_i\_i∈[n_cal] ). Here QuantilelinearQuantile^linear denotes the default linear-interpolation convention of numpy.quantile. This numerical implementation uses an interpolated empirical quantile rather than the exact conformal order statistic. Empirical mode diagnostics (Figure 3). The covariate and noisy-logit DGP, training size ntrain=2000n_train=2000, and fitted logistic model are the same as for Table 1. We compare empirical-mode and population-mode PICPI using 100100 equal-width score bins, with δ=0.1δ=0.1 for population mode. We vary the calibration size over 103,3×103,104,3×104,105,3×105,106,3×106,107\10^3,3× 10^3,10^4,3× 10^4,10^5,3× 10^5,10^6,3× 10^6,10^7\, use neval=106n_eval=10^6, and run 5050 replications for each calibration size. In these diagnostics, the held-out true probability is the Gaussian-averaged marginal probability ξ[σ(β0+x⊤β+ξ)]E_ξ[σ( _0+x β+ξ)]. B.4. Setup for Task 2 Coverage sweep for two DGPs (Figure 6) Both use a d=10d=10 dimensional covariate and G=4G=4 classes. Given the true class-probability vector q(x)∈[0,1]Gq(x)∈[0,1]^G, labels are drawn as Y|X=x∼Categorical(q(x)).Y X=x (q(x) ). Each replication uses disjoint splits of sizes ntrain=ncal=neval=2000,n_train=n_cal=n_eval=2000, and all curves are averages over 10001000 independent replications. PICPI uses num_bin=100=100 equally spaced candidate grid points in empirical mode, while the calibration-based and fixed-width baselines use 100100 bins. Candidate intervals from PICPI, split conformal prediction, calibration-based intervals, and fixed-width binning are passed to the classwise interval-selection solver (Algorithm 2). For split conformal prediction, each class uses at most 200200 interval centers sampled without replacement from its calibration scores. Simultaneous confidence intervals and conformal classification instead use a calibration threshold, as described in Section 5. The coverage sweep evaluates all methods at 1616 equally spaced target coverage levels 1−α∈0,1/15,…,11-α∈\0,1/15,…,1\. DGP 1: Linear-Gaussian. The covariate is drawn from X∼N(0,I10),X N(0,I_10), and the true class probabilities follow a linear softmax model, q(x)=softmax(Wx+b),q(x)=softmax(Wx+b), with a sparse weight matrix W=[1.4−1.00.800000000−0.701.10.800000000.900−1.21.0000−0.60000000.8−0.91.2]W= bmatrix1.4&-1.0&0.8&0&0&0&0&0&0&0\\ 0&-0.7&0&1.1&0.8&0&0&0&0&0\\ 0&0&0.9&0&0&-1.2&1.0&0&0&0\\ -0.6&0&0&0&0&0&0&0.8&-0.9&1.2 bmatrix and an intercept b=(0.25,−0.20, 0.15,−0.20)⊤b=(0.25,\,-0.20,\,0.15,\,-0.20) . The multinomial logistic working model is therefore correctly specified. DGP 2: Nonlinear-mixture. The covariate is drawn from a four-component latent-factor mixture with component weights (0.32, 0.24, 0.27, 0.17).(0.32,\,0.24,\,0.27,\,0.17). Equivalently, a latent component label is first drawn as k∼Categorical(0.32, 0.24, 0.27, 0.17).k (0.32,\,0.24,\,0.27,\,0.17). Conditional on component k∈1,…,4k∈\1,…,4\, an initial draw X~∈ℝ10 X ^10 is generated as X~=μk+Λ⊤Z+diag(σk)h(Z)Sε,Z∼N(0,I2),ε∼N(0,I10), X= _k+ Z+diag( _k)\,h(Z)\,S\, , Z N(0,I_2), N(0,I_10), where S=7/χ72S= 7/χ^2_7 induces Student-t (heavy) tails and h(z)=1+0.25tanh(z1+0.45z2)h(z)=1+0.25 (z_1+0.45z_2) induces heteroscedasticity in the latent factors. The component-specific location and scale vectors μk,σk∈ℝ10 _k, _k ^10 are the rows of M=[−1.40.5−0.90.2−0.30.6−0.50.4−0.70.3−0.2−1.00.7−0.60.8−0.40.3−0.50.6−0.80.90.6−0.21.1−0.70.40.90.7−0.10.51.6−0.31.00.50.4−0.8−1.0−0.60.9−0.4], M= bmatrix-1.4&0.5&-0.9&0.2&-0.3&0.6&-0.5&0.4&-0.7&0.3\\ -0.2&-1.0&0.7&-0.6&0.8&-0.4&0.3&-0.5&0.6&-0.8\\ 0.9&0.6&-0.2&1.1&-0.7&0.4&0.9&0.7&-0.1&0.5\\ 1.6&-0.3&1.0&0.5&0.4&-0.8&-1.0&-0.6&0.9&-0.4 bmatrix, V=[0.600.450.550.400.500.350.450.420.500.380.450.650.500.600.550.400.350.480.440.570.550.500.700.450.600.500.550.520.460.500.500.400.600.550.450.650.600.580.540.47], V= bmatrix0.60&0.45&0.55&0.40&0.50&0.35&0.45&0.42&0.50&0.38\\ 0.45&0.65&0.50&0.60&0.55&0.40&0.35&0.48&0.44&0.57\\ 0.55&0.50&0.70&0.45&0.60&0.50&0.55&0.52&0.46&0.50\\ 0.50&0.40&0.60&0.55&0.45&0.65&0.60&0.58&0.54&0.47 bmatrix, and the shared factor-loading matrix is Λ=[0.95−0.550.350.45−0.300.40−0.200.32−0.280.240.300.60−0.450.200.55−0.350.65−0.250.48−0.42]. = bmatrix0.95&-0.55&0.35&0.45&-0.30&0.40&-0.20&0.32&-0.28&0.24\\ 0.30&0.60&-0.45&0.20&0.55&-0.35&0.65&-0.25&0.48&-0.42 bmatrix. The covariate X is then obtained from X~ X by fixed deterministic warps. Starting from x=x~+0.22sin(x~1)δkx= x+0.22 ( x_1)\, _k, the coordinates are updated sequentially (each line using the current values) via x3 x_3 ←x3+0.28(x12−1), ← x_3+0.28\,(x_1^2-1), x5 x_5 ←x5+0.16x2x4, ← x_5+0.16\,x_2x_4, x7 x_7 ←x7+0.14x3x6, ← x_7+0.14\,x_3x_6, x8 x_8 ←x8+0.18sin(x1x3), ← x_8+0.18 (x_1x_3), x9 x_9 ←x9+0.15x5x7, ← x_9+0.15\,x_5x_7, x10 x_10 ←x10+0.12(x82−1), ← x_10+0.12\,(x_8^2-1), followed by the component offset x←x+okx← x+o_k and the final warps x6 x_6 ←x6+0.20sgn(x1)|x1|+10−6, ← x_6+0.20\,sgn(x_1) |x_1|+10^-6, x2 x_2 ←x2−0.15tanh(x7−x5), ← x_2-0.15 (x_7-x_5), x9 x_9 ←x9+0.18sgn(x8)|x8|+10−6, ← x_9+0.18\,sgn(x_8) |x_8|+10^-6, x10 x_10 ←x10−0.12tanh(x9−x6). ← x_10-0.12 (x_9-x_6). The drift and offset vectors δk,ok∈ℝ10 _k,o_k ^10 are the rows of D=[0.30−0.200.100.000.12−0.080.050.07−0.060.04−0.100.28−0.120.08−0.060.14−0.04−0.090.10−0.050.08−0.120.26−0.160.10−0.060.180.11−0.040.09−0.140.06−0.080.24−0.100.20−0.22−0.050.12−0.13], D= bmatrix0.30&-0.20&0.10&0.00&0.12&-0.08&0.05&0.07&-0.06&0.04\\ -0.10&0.28&-0.12&0.08&-0.06&0.14&-0.04&-0.09&0.10&-0.05\\ 0.08&-0.12&0.26&-0.16&0.10&-0.06&0.18&0.11&-0.04&0.09\\ -0.14&0.06&-0.08&0.24&-0.10&0.20&-0.22&-0.05&0.12&-0.13 bmatrix, O=[−0.120.06−0.040.02−0.050.03−0.020.04−0.030.010.05−0.100.08−0.060.09−0.040.03−0.020.05−0.070.100.04−0.060.11−0.080.060.120.03−0.090.080.16−0.020.100.070.05−0.11−0.14−0.080.07−0.05]. O= bmatrix-0.12&0.06&-0.04&0.02&-0.05&0.03&-0.02&0.04&-0.03&0.01\\ 0.05&-0.10&0.08&-0.06&0.09&-0.04&0.03&-0.02&0.05&-0.07\\ 0.10&0.04&-0.06&0.11&-0.08&0.06&0.12&0.03&-0.09&0.08\\ 0.16&-0.02&0.10&0.07&0.05&-0.11&-0.14&-0.08&0.07&-0.05 bmatrix. The resulting marginal covariate law is non-Gaussian, multimodal, and dependent across coordinates. The true class probabilities follow a softmax model whose linear logits are rescaled by a covariate-dependent scale function, q(x)=softmax(Wx+bs(x)), q(x)=softmax\! ( Wx+bs(x) ), s(x)=clip(0.86+0.22tanh(0.60x1−0.35x4+0.20x9), 0.65, 1.40), s(x)=clip (0.86+0.22 (0.60x_1-0.35x_4+0.20x_9),\ 0.65,\ 1.40 ), where the weight matrix W∈ℝ4×10W ^4× 10 and intercept b∈ℝ4b ^4 were generated once (by combining i.i.d. Gaussian entries with structured rank-one terms drawn from a fixed seed) and held constant across all replications; rounded to two decimals, W=[−1.491.57−0.40−1.121.37−0.670.81−0.920.960.07−0.100.32−0.75−0.450.29−0.140.54−0.37−0.30−0.100.32−0.650.47−0.050.11−0.12−0.440.320.160.571.53−1.481.061.230.051.37−0.500.910.180.41]W= bmatrix-1.49&1.57&-0.40&-1.12&1.37&-0.67&0.81&-0.92&0.96&0.07\\ -0.10&0.32&-0.75&-0.45&0.29&-0.14&0.54&-0.37&-0.30&-0.10\\ 0.32&-0.65&0.47&-0.05&0.11&-0.12&-0.44&0.32&0.16&0.57\\ 1.53&-1.48&1.06&1.23&0.05&1.37&-0.50&0.91&0.18&0.41 bmatrix and b=(0.30, 0.16,−0.05,−0.41)⊤b=(0.30,\,0.16,\,-0.05,\,-0.41) . Because the log-odds are nonlinear in x through s(x)s(x), the multinomial logistic working model is misspecified and thus is not expected to be well calibrated under this DGP. B.5. Setup for Figure 2 The covariate and noisy-label DGP are the same as in the univariate interval visualization of Section B.3. We use Algorithm 1 with 100100 bins and δ=0.1δ=0.1 and 100100 repetitions. The noise-free target is p∗(x)=σ(2x)p^*(x)=σ(2x). We fit an unpenalized logistic-regression estimator once using a training set of size 10,00010,000, fix an evaluation set of size 10,00010,000, and vary the calibration-set size over 103,3×103,104,3×104,105,3×105,106,3×106,107\10^3,3× 10^3,10^4,3× 10^4,10^5,3× 10^5,10^6,3× 10^6,10^7\. 1. fit a logistic-regression estimator p p on the training fold; 2. apply the calibration procedure on the calibration fold to obtain an interval collection ℐ=[aj,bj]jI=\[a_j,b_j]\_j; 3. evaluate the interval system on the independent evaluation fold by computing, for each point, the minimal width among all intervals that contain the target probability p∗(Xi)p^*(X_i), with width 1.01.0 assigned if no interval contains p∗(Xi)p^*(X_i). Appendix C Additional Algorithms Algorithm 4 selects a disjoint subset for intervals returned by Algorithm 1. The DisjointCalibration procedure first imposes a maximum length and then uses weighted interval scheduling to select a non-overlapping subcollection with maximum covered length. The MaxCoverageIntervals procedure solves a weighted interval scheduling problem to select a set of non-overlapping intervals that maximizes total coverage, where each interval is weighted by its length. It sorts intervals by end time, precomputes the last compatible interval for each candidate, and then uses dynamic programming to recursively compute the optimal value, followed by a backtracking step to recover the selected intervals. Algorithm 4 Maximum-coverage non-overlapping intervals (weighted interval scheduling) 1: procedure DisjointCalibration(length bound l, PICPIs ℐpI_p) 2: ⊳ Restrict PICPIs to intervals of length ≤l≤ l and maximize covered length 3: ℐp(≤l)←∅I_p^(≤ l)← 4: for each interval [a,b]∈ℐp[a,b] _p do 5: if b−a≤lb-a≤ l then 6: ℐp(≤l)←ℐp(≤l)∪[a,b]I_p^(≤ l) _p^(≤ l)∪\[a,b]\ 7: end if 8: end for 9: if ℐp(≤l)=∅I_p^(≤ l)= then 10: Return ∅ ⊳ No interval satisfies the length constraint 11: end if 12: Relabel ℐp(≤l)I_p^(≤ l) as [si,ei]i=1m\[s_i,e_i]\_i=1^m 13: ⋆←MaxCoverageIntervals([si,ei]i=1m)S ← MaxCoverageIntervals (\[s_i,e_i]\_i=1^m ) ⊳ Select a disjoint set intervals of maximal coverage 14: Return ⋆S 15: end procedure 16: 17: procedure MaxCoverageIntervals(Intervals [si,ei]i=1m\[s_i,e_i]\_i=1^m) ⊳ Find a maximum-length set of non-overlapping intervals 18: Sort the intervals so that e1≤e2≤⋯≤eme_1≤ e_2≤·s≤ e_m 19: for j←1j← 1 to m do 20: wj←ej−sjw_j← e_j-s_j ⊳ Weight = interval length 21: end for 22: for j←1j← 1 to m do 23: p(j)←maxi<j:ei≤sjp(j)← \\,i<j:e_i≤ s_j\,\, or 00 if no such i exists ⊳ Binary search over e1,…,ej−1\e_1,…,e_j-1\ 24: end for 25: Initialize array OPT[0..m]OPT[0..m] ⊳ Dynamic Programming 26: OPT[0]←0OPT[0]← 0 27: for j←1j← 1 to m do 28: OPT[j]←max(wj+OPT[p(j)],OPT[j−1])OPT[j]← (w_j+OPT[p(j)], \ OPT[j-1] )⊳ Maximum coverage using 1,2,…,j\1,2,…,j\ 29: end for 30: Initialize selected set ⋆←∅S ← 31: j←mj← m 32: while j>0j>0 do 33: if wj+OPT[p(j)]>OPT[j−1]w_j+OPT[p(j)]>OPT[j-1] then 34: ⋆←⋆∪[sj,ej]S ∪\[s_j,e_j]\ 35: j←p(j)j← p(j) 36: else 37: j←j−1j← j-1 38: end if 39: end while 40: Return ⋆S 41: end procedure