Scaling-Score Conformal Prediction for Multi-Target Regression Sylvain Rousseau Soundouss Messoudi
arXiv:2609.17091v1 [cs.AI] 15 Sep 2026
Université de Technologie de Compiègne, CNRS, Heudiasyc, Compiègne, France [email protected], [email protected]
Abstract Multi-target regression requires a model to simultaneously predict several related outputs. Conformal prediction provides distribution-free, finite-sample marginal coverage guarantees, but extending these to joint multi-dimensional regions in a model-agnostic, sample-efficient manner remains challenging: max-aggregation ignores scale differences, copula-based methods are only asymptotically valid, rectangular methods typically split the calibration set, and quantile- or density-based methods require training a specialised model beyond a plain point predictor. We propose the scaling-score conformal method, which is model-agnostic (requires only componentwise absolute residuals), uses a single calibration set, and yields four nested output types: an outer rectangle (SCO) with valid joint coverage, the exact set Rα , a staircase (SC2 ) overapproximation of Rα , and an inner rectangle (SCI). A single hyperparameter γ ∈ (0, 1) controls the base-rectangle quantile level independently of α. We prove downward-closedness and a rectangular sandwich bound and derive a closed-form outer rectangle. Experiments on 29 realworld datasets confirm valid joint coverage; SC2 with γ = 1−α consistently achieves competitive volume relative to baselines, with the advantage growing with output dimension d. Keywords: conformal prediction, multi-target regression, joint prediction regions, uncertainty quantification, nonconformity score
1
Introduction
Modern applications increasingly require the simultaneous prediction of several correlated response variables. Such settings arise whenever the response is intrinsically multivariate: energy systems, environmental monitoring, finance, and supply chains all produce outputs whose components must be predicted jointly. In all these settings, downstream decisions depend on the joint distribution of the outputs, making a marginal error bar for each target individually insufficient. The relevant object is a prediction region C(x) ⊆ Rd that simultaneously covers all d components with a specified probability. Conformal prediction (Vovk et al., 2005) provides prediction regions with a finite-sample, distribution-free coverage guarantee: under exchangeability of the calibration and test data, P[Yn+1 ∈ C(Xn+1 )] ≥ 1 − α for any miscoverage level α ∈ (0, 1), irrespective of the data-generating distribution. The method is model-agnostic: it wraps any point predictor, and the guarantee requires no distributional assumption beyond exchangeability. Extending these guarantees to d-dimensional outputs in a simultaneously valid, model-agnostic, and sample-efficient way is non-trivial, and current methods each suffer from at least one key weakness. The simplest multi-target approach, max aggregation, takes the maximum of per-target absolute residuals; it yields valid coverage but produces axis-aligned hypercubes that are over-conservative when targets have heterogeneous scales. Copula-based methods (Messoudi et al., 2021) model the 1
joint dependence structure through a copula and directly construct a prediction region; coverage guarantees are at best asymptotic. CHR (Sampson and Chan, 2024) achieves valid rectangular regions with per-target width adaptation, but requires splitting the calibration set in two, reducing the effective sample size for the conformal quantile. Transductively Standardized Conformal Prediction (Fan and Sesia, 2025) is valid and model-agnostic with a single calibration set, but standardizes by the empirical mean and standard deviation of calibration residuals, which implicitly assumes these moments are well-defined and well-estimated; normalizing by calibration quantiles instead avoids this assumption, since a quantile depends only on the relative order of the observations rather than on every calibration point. Methods based on conformalized quantile regression (Romano et al., 2019), conditional density estimation (Izbicki et al., 2020), or normalizing flows require a trained model beyond a plain point predictor. We introduce a scaling score that, for a given residual vector, measures the minimum factor by which a per-target base vector must be scaled to contain it. The base vector is derived from per-target calibration quantiles at a tunable level, encoding the relative scale of each output coordinate. The base vector is defined transductively: it depends on the test residual through a clamping formula, so that all calibration and test scores share the same base and exchangeability is preserved. From a single calibration set, our method simultaneously produces four nested output types (Figure 2(b)): an outer rectangle (SCO), a single axis-aligned hyperrectangle that is valid but conservative; the exact prediction set Rα , which is the tightest valid prediction set but defined implicitly via a cell decomposition; a staircase region (SC2 ), a union of at most 2d axis-aligned hyperrectangles that is a valid over-approximation of Rα , computable at calibration time; and an inner rectangle (SCI), the most compact of the four, tighter even than Rα , but without a formal coverage guarantee. Our contributions are threefold. First, we propose a model-agnostic, single-calibration-set conformal score for multi-target regression that simultaneously produces four nested output types, three of which carry valid joint coverage guarantees. Second, we derive an exact cell-based analytical characterisation of the prediction region, rectangular sandwich bounds with closed-form expressions, and a proof of downward-closedness. Third, experiments on 29 real-world datasets confirm valid joint coverage and competitive efficiency across a range of target dimensions. Section 2 reviews conformal prediction and existing multi-target methods. Section 3 establishes notation and problem setup. Section 4 develops the scaling-score method. Section 5 reports experiments on 29 real-world datasets. Section 6 concludes. Proofs are collected in Appendix A; dataset details appear in Appendix B; full experimental figures appear in Appendix C.
2
Background and Related Work
2.1
Conformal prediction for regression
Let (X1 , Y1 ), . . . , (Xn , Yn ), (Xn+1 , Yn+1 ) be exchangeable random variables taking values in X × Rd . Given a pre-trained predictor fˆ : X → Rd , a nonconformity score s(x, y) measures how poorly fˆ(x) predicts y. In split (inductive) conformal prediction (Papadopoulos et al., 2002; Vovk et al., 2005), the predictor is fitted on training data and the score evaluated on a held-out calibration set Dcal = {(Xi , Yi )}ni=1 . Setting si = s(Xi , Yi ), the conformal quantile is q̂1−α = ⌈(1 − α)(n + 1)⌉-th smallest value in {s1 , . . . , sn , +∞}, and the prediction region C(x) = {y : s(x, y) ≤ q̂1−α } satisfies P[Yn+1 ∈ C(Xn+1 )] ≥ 1 − α by a simple rank argument (Lei et al., 2018). 2
Standard nonconformity scores for univariate regression include the absolute residual |y − fˆ(x)| and the conformalized quantile regression score (Romano et al., 2019). Normalized scores (Papadopoulos et al., 2002; Lei et al., 2018) divide the residual by a local scale estimate to improve conditional coverage.
2.2
Conformal prediction for multi-target regression
With d-dimensional outputs, the joint coverage event {Yn+1 ∈ C(Xn+1 )} requires the entire response vector to lie within the prediction region C(Xn+1 ) ⊆ Rd . Achieving at least 1 − α joint coverage while keeping the region small is harder than marginal coverage of each component individually. The simplest approach uses the Max (ℓ∞ ) score s(x, y) = maxk∈[d] |y (k) − fˆ(k) (x)|, which is valid by the standard conformal argument and produces a hypercube of radius q̂1−α centered at fˆ(x). However, the hypercube is driven by the coordinate with the largest residuals; when targets have heterogeneous scales, the resulting region is over-conservative for the other targets. The Bonferroni correction applied coordinate-wise achieves 1 − dα′ joint coverage at level α′ = α/d per dimension, but this is even more conservative. Copula-based methods (Messoudi et al., 2021) model the joint dependence structure of the residuals directly. Parametric copulas assume a specific family for the joint distribution of errors; (k) the empirical copula instead uses the rank transform, mapping each error coordinate ei through its marginal ECDF F̂k to place all targets on a common [0, 1] scale. Both approaches produce axisaligned rectangular regions. However, neither defines a nonconformity score in the split-conformal sense, so finite-sample coverage is not guaranteed. The empirical copula achieves only asymptotic validity as n → ∞; parametric copulas additionally require that the residuals follow the assumed copula distribution. CHR (Sampson and Chan, 2024) (Conformal Multi-Target Hyperrectangles) is a valid splitconformal method. It uses a two-stage calibration: a first portion of the calibration set estimates per-target prediction intervals (capturing scale differences across targets), and a second portion computes a scaled nonconformity score on top of the estimated intervals. The method achieves finite-sample joint coverage and asymptotic balance (equal marginal coverage across targets). The main limitation is that the required calibration split reduces the effective sample size available for the conformal quantile, which can inflate the prediction region when n is moderate or d is large. TSCP and TSCP-R (Fan and Sesia, 2025) (Transductively Standardized Conformal Prediction) standardize the residual vector by the empirical mean and standard deviation of calibration residuals in a transductive fashion, preserving exchangeability. Both are valid, model-agnostic, and use a single calibration set. The TSCP-R variant restricts to a rectangular output region. These methods are the closest antecedent to our work; the key difference is that TSCP standardizes by location and scale (mean ± std), while our method normalizes purely by calibration quantiles, giving more direct control over the coverage level of the base rectangle. Density-based conformal methods (Izbicki et al., 2020) can produce non-rectangular regions that capture correlation geometry, and may therefore achieve smaller volume than rectangular methods. However, they require learning a conditional density estimator in addition to a point predictor, which adds modeling complexity. Non-rectangular regions are also harder to interpret and communicate than simple axis-aligned boxes. Dheur et al. (2025) provide a unified comparative study of several such methods.
3
3
Problem Setup and Notation
Let X be a feature space and Y ∈ Rd the response, and let fˆ : X → Rd be a pre-trained point predictor. We observe a calibration set Dcal = {(Xi , Yi )}ni=1 and assume that {(Xi , Yi )}n+1 i=1 are exchangeable, the sole assumption underlying all validity guarantees. Define the component-wise absolute calibration errors ei = |Yi − fˆ(Xi )| ∈ Rd+ ,
i = 1, . . . , n.
For a candidate test response y ∈ Rd define en+1 = en+1 (y) = |y − fˆ(Xn+1 )|. The goal is to construct a prediction region C(Xn+1 ) ⊆ Rd such that P[Yn+1 ∈ C(Xn+1 )] ≥ 1 − α. Write [d] = {1, . . . , d}. For a vector a ∈ Rd , we write a(k) for its k-th coordinate; the same convention applies to b, e, Y , and all other vectors throughout. For a, b ∈ Rd , write a ⪯ b for the (k) component-wise partial order: a(k) ≤ b(k) for all k ∈ [d]. Let e(j) denote the j-th order statistic of (k)
(k)
{e1 , . . . , en }.
4
The Scaling-Score Method
Our method constructs a rectangular prediction region from a per-target base rectangle whose side lengths encode the relative scale of each output. The base rectangle is then inflated by the smallest scaling factor (Figure 1(a)) that covers a 1 − α fraction of the calibration points. The main subtlety is that the base rectangle depends on the test error through a clamping formula, and handling this dependency correctly while preserving exchangeability is the central challenge of the construction. We describe the construction in four steps: the scaling score (Section 4.1), the lower and upper base rectangles (Section 4.2), the cell decomposition (Section 4.3), and the prediction region characterisation and its properties (Section 4.4).
4.1
Scaling score
Given a base vector b ∈ Rd++ , the scaling score measures the minimum factor by which b must be scaled to contain a residual e (Figure 1(a)). Definition 1 (Scaling score) For a base vector b ∈ Rd++ and any e ∈ Rd+ , the scaling score is Q the minimum dilation factor needed to fit e inside the base rectangle dk=1 [0, b(k) ]: e(k) . k∈[d] b(k)
σ(e; b) = max
(1)
Geometrically, bounding σ(e; b) by a threshold λ is equivalent to requiring that e lies in the axis-aligned rectangle obtained by scaling b coordinate-wise by λ: σ(e; b) ≤ λ ⇐⇒ e ⪯ λb ⇐⇒ e ∈
d Y
0, λ b(k) .
(2)
k=1
4.2
Lower and upper base rectangles
Fix a level γ ∈ (0, 1) and set m′ = ⌈γ(n + 1)⌉. Define the lower and upper base rectangles via their per-coordinate bounds (k) (k) Lk = e(m′ −1) , Uk = e(m′ ) , k ∈ [d]. (3) 4
(k)
(k)
We require 2 ≤ m′ ≤ n so that Lk = e(m′ −1) is a well-defined order statistic (m′ ≥ 2) and Uk = e(m′ ) is finite (m′ ≤ n); in practice we clip m′ = clip(⌈γ(n + 1)⌉, 2, n). By construction 0 ≤ Lk ≤ Uk for all k. The base rectangle is constructed transductively: its side lengths (each in [Lk , Uk ]) depend on the test residual en+1 through the rank it would occupy in the augmented calibration set, keeping all scores exchangeable. Definition 2 (Transductive base vector) statistic of the pooled set in each coordinate:
The base vector bn+1 ∈ Rd++ is the m′ -th order
(k)
(k)
bn+1 = m′ -th order statistic of {ei }n+1 i=1 , The associated base rectangle is Bn+1 =
k ∈ [d].
(k) k=1 [0, bn+1 ].
Qd
(k)
The base vector depends on en+1 only through the rank of en+1 relative to the calibration values. Specifically, (k) if en+1 ≤ Lk , Lk (k) (k) bn+1 = e(k) (4) n+1 if Lk < en+1 ≤ Uk , (k) Uk if en+1 > Uk , (k)
so Lk ≤ bn+1 ≤ Uk for all en+1 . Figure 1(b) illustrates the thresholds Lk and Uk , the resulting cell structure for d = 2, and the clamping relationship between en+1 and bn+1 . In the transductive setting, all n + 1 scores share the same base bn+1 : σi (en+1 ) = σ(ei ; bn+1 ),
i = 1, . . . , n + 1.
The following theorem is the main coverage guarantee of the method. It states that the prediction region defined via a conformal quantile achieves the desired joint marginal coverage, under the sole assumption of exchangeability. Theorem 3 (Validity) Let q = ⌈(1 − α)(n + 1)⌉ and let λ̂(en+1 ) denote the q-th smallest value in {σ1 (en+1 ), . . . , σn (en+1 ), +∞}. Define the prediction region in error space Rα = en+1 ∈ Rd+ : σn+1 (en+1 ) ≤ λ̂(en+1 ) . If (e1 , . . . , en , en+1 ) are exchangeable, then P(en+1 ∈ Rα ) ≥ 1 − α. Proof See the proof in Appendix A. The validity guarantee in Proposition 3 establishes that Rα achieves the desired coverage, but the definition Rα = {en+1 : σn+1 (en+1 ) ≤ λ̂(en+1 )} is an implicit equation: the construction is transductive, meaning the base vector bn+1 depends on en+1 through the clamping formula (4), so every candidate test error induces a different calibration ranking and, potentially, a different quantile threshold. Computing Rα explicitly therefore requires a dedicated analysis. To handle this, we decompose Rd+ into a finite collection of cells and treat each cell separately: conditioned on the test error en+1 belonging to a given cell, the base vector bn+1 is constant throughout that cell, so the calibration scores and the conformal quantile can be computed analytically.
5
e(2)
high
e(2)
C{1},∅,{2}
C∅,{1},{2}
mid
C{1},{2},∅
L2
b2
en+1
bn+1
U2 σ = 1.33
C∅,∅,{1,2}
C∅,{1,2},∅
C∅,{2},{1}
bn+1
low
en+1
C{1,2},∅,∅
σ = 1.5 0
σ = 0.6
0
b1
low
C{2},{1},∅
L1
mid
corner cell (Imid
e(1)
C{2},∅,{1}
U1
e(1) high
= ∅, 22 = 4 cells)
mixed cell (Imid 6= ∅, 32 − 22 = 5 cells)
(a) Scaling score for three residuals (d = 2, (b) The 32 = 9 cells of the decomposition for base vector (b1 , b2 )). The score σ is the d = 2. The thresholds L1 , U1 and L2 , U2 smallest factor by which the base rectandivide R2+ into a 3 × 3 grid. Corner cells (blue) have Imid = ∅; mixed cells (red) have gle must be scaled to contain the point. |Imid | ≥ 1. The clamping relationship beAll scaled-rectangle corners lie on the ray tween en+1 and bn+1 is also visible: each cothrough (b1 , b2 ) (gray line). ordinate of bn+1 equals the clamp of the corresponding coordinate of en+1 to the interval [Lk , Uk ]. Figure 1: Scaling score and cell decomposition (d = 2).
4.3
Cell decomposition and scores
The thresholds Lk and Uk partition the k-th coordinate axis into three intervals: below Lk (lower region), between Lk and Uk (middle region), and above Uk (upper region). Taking the product over all d coordinates yields a partition of Rd+ into 3d cells. Each cell is characterised by specifying, for every coordinate k ∈ [d], whether that coordinate falls in the lower, middle, or upper interval. We encode this by a triple (Ilow , Imid , Ihigh ), a partition of [d], where Ilow collects the coordinates in the lower region, Imid those in the middle region, and Ihigh those in the upper region. Figure 1(b) illustrates the 32 = 9 cells for d = 2. This is captured formally by the following definition. Definition 4 (Cells) For a partition (Ilow , Imid , Ihigh ) of [d], the corresponding cell is e(k) ≤ Lk , ∀k ∈ Ilow , CIlow ,Imid ,Ihigh = e ∈ Rd+ : Lk < e(k) ≤ Uk , ∀k ∈ Imid , . (k) e > Uk , ∀k ∈ Ihigh Assuming the test error en+1 belongs to cell CIlow ,Imid ,Ihigh , the base vector (4) is fully determined (k)
(k)
(k)
(k)
by the partition: bn+1 = Lk for k ∈ Ilow , bn+1 = en+1 for k ∈ Imid , and bn+1 = Uk for k ∈ Ihigh . Figure 1(b) illustrates how bn+1 depends on en+1 in the three-region partition for each coordinate. 6
Cells with Imid = ∅ are called corner cells: the base vector does not depend on the test error throughout the cell, and we denote its value bIlow ,Ihigh . Cells with Imid ̸= ∅ are mixed cells. There are 2d corner cells and 3d − 2d mixed cells. Figure 1(b) illustrates both types for d = 2. Given the scaling score σ(e; b) = maxk e(k) /b(k) and the cell formula for bn+1 , the scores σn+1 and σi take a simple explicit form on each cell, as summarised in the following lemma. Lemma 5 (Scores in any cell) Let en+1 ∈ CIlow ,Imid ,Ihigh . With the convention max∅ := 0: (k) (k) en+1 en+1 σn+1 (en+1 ) = max max , max 1, max , k∈Ilow Lk k∈Ihigh Uk k∈Imid (k) (k) (k) e ei e σi (en+1 ) = max max i , max (k) , max i , k∈Ilow Lk k∈Imid e k∈Ihigh Uk
(5)
i ≤ n.
(6)
n+1
When en+1 belongs to a corner cell (Imid = ∅), σi (en+1 ) is constant on the cell. When en+1 belongs to a mixed cell (Imid ̸= ∅), σn+1 (en+1 ) ≥ 1, and the mid terms in σi (en+1 ) are non-increasing in (k) each en+1 , k ∈ Imid . Proof See the proof in Appendix A.
4.4
Prediction region characterisation
Having partitioned Rd+ into cells, we now characterise Rα by working one cell at a time. The key observation is the following: assuming the test error en+1 belongs to a given cell, the base vector bn+1 is fully determined and remains constant for all test errors confined to that cell (in corner cells) or depends on en+1 only through an explicit formula (in mixed cells). As a consequence, all calibration scores {σi (en+1 )}ni=1 are determined, the conformal quantile λ̂(bn+1 ) can be computed analytically. For corner cells, the base vector is fixed and the region boundary has an explicit closed form. For mixed cells, the base vector still depends on en+1 through its mid coordinates, so the boundary retains an implicit form; we characterise it analytically and provide sandwich bounds. We treat corner cells and mixed cells in turn. 4.4.1
Corner cells
When the test error en+1 belongs to a corner cell CIlow ,∅,Ihigh , the base vector bIlow ,Ihigh does not depend on the test error, so all calibration scores (k) (k) e e I ,I σi low high = max max i , max i k∈Ilow Lk k∈Ihigh Uk
are independent of en+1 , and the conformal quantile λ̂Ilow ,Ihigh is a fixed number. Proposition 6 (Exact region on corner cells) Let λ̂ = λ̂Ilow ,Ihigh . The exact prediction region restricted to the corner cell is Y Y Rα ∩ CIlow ,∅,Ihigh = 0, min(λ̂, 1) Lk × Uk , λ̂ Uk , (7) k∈Ilow
k∈Ihigh
with the convention that (Uk , λ̂ Uk ] = ∅ when λ̂ ≤ 1. 7
Proof See the proof in Appendix A. Two boundary cases of Proposition 6 are worth noting. When Ihigh ̸= ∅, the intersection (7) is non-empty if and only if λ̂ > 1: the conformal region extends beyond Uk only when the quantile exceeds 1. When Ihigh = ∅ (all-low corner cell), the region always contains the origin and equals Q the full cell dk=1 [0, Lk ] if and only if λ̂ ≥ 1. 4.4.2
Mixed cells
When the test error en+1 belongs to a mixed cell CIlow ,Imid ,Ihigh (with Imid ̸= ∅), the base vector bn+1 (k)
(k)
depends on en+1 through the mid coordinates (since bn+1 = en+1 for k ∈ Imid ). As a consequence, the calibration scores {σi (en+1 )} and the conformal quantile also depend on en+1 , but only through (Imid ) its Imid coordinates; we write λ̂(en+1 ) to make this dependence explicit. Proposition 7 (Exact region on mixed cells) Assuming the test error en+1 belongs to a mixed cell CIlow ,Imid ,Ihigh : 1. If Ihigh = ∅: σn+1 (en+1 ) = 1, and (Imid ) Rα ∩ CIlow ,Imid ,∅ = en+1 ∈ CIlow ,Imid ,∅ : λ̂(en+1 )≥1 . (k)
2. If Ihigh ̸= ∅: σn+1 (en+1 ) = maxk∈Ihigh en+1 /Uk , and (k) (Imid ) Rα ∩ CIlow ,Imid ,Ihigh = en+1 ∈ CIlow ,Imid ,Ihigh : en+1 ≤ λ̂(en+1 ) Uk ∀k ∈ Ihigh . In both cases, the low coordinates are unconstrained within the cell: for k ∈ Ilow , cell membership (k) (k) (k) (k) already forces en+1 ≤ Lk = bn+1 , so the low terms en+1 /bn+1 ≤ 1 never drive σn+1 , and σn+1 ≤ λ̂ (k)
imposes no further constraint on them beyond en+1 ≤ Lk . Proof See the proof in Appendix A. (I
)
mid Item 1 of Proposition 7 (Ihigh = ∅) reduces to the single condition λ̂(en+1 ) ≥ 1: since σn+1 = 1 throughout this cell, inclusion in Rα depends only on whether the conformal quantile reaches 1. (Imid ) In item 2 (Ihigh ̸= ∅), the condition additionally constrains the high coordinates via λ̂(en+1 ).
(I
)
mid Both cases are implicit: λ̂(en+1 ) depends on en+1 through the mid coordinates, and there is in
(I
)
mid general no closed-form description of Rα on mixed cells. Nevertheless, λ̂(en+1 ) can be bracketed analytically, as the following proposition shows.
Proposition 8 (Monotonicity and interpolation) On a mixed cell CIlow ,Imid ,Ihigh : (I
)
mid 1. λ̂(en+1 ) is non-increasing in each mid coordinate.
(k)
2. At the boundary en+1 → L+ k for all k ∈ Imid : λ̂ → λ̂Ilow ∪Imid , Ihigh . (k)
3. At en+1 = Uk for all k ∈ Imid : λ̂ = λ̂Ilow , Imid ∪Ihigh . (I
)
mid Hence λ̂Ilow , Imid ∪Ihigh ≤ λ̂(en+1 ) ≤ λ̂Ilow ∪Imid , Ihigh throughout the cell.
8
Proof See the proof in Appendix A. Propositions 7 and 8 together imply that the boundary of Rα is continuous across cell boundaries: items 2 and 3 of Proposition 8 show that λ̂ matches the corner-cell values at each face, so no jump can occur when crossing from one cell to an adjacent one. Together, Propositions 6, 7 and 8 give a complete analytical picture of Rα : on corner cells it reduces to an explicit rectangle, and on mixed cells it is characterised by the two-item structure of Proposition 7. While Rα is not itself a hyperrectangle, it can be sandwiched between two axisaligned hyperrectangles (the inner and outer rectangles of Section 4.5), making it both principled and tractable. 4.4.3
Structural properties
S Assembling the cell-by-cell characterisations above gives the full prediction region Rα = (Ilow ,Imid ,Ihigh ) Rα ∩ CIlow ,Imid ,Ihigh . A key global property of this union is that it is downward-closed with respect to the component-wise order. Intuitively, if a residual en+1 is already within the conformal region, then any smaller residual e′ ⪯ en+1 must also lie within it. This follows from the monotonicity of the score ratio x 7→ x/ clamp(x, L, U ), which is non-decreasing on R+ for any fixed 0 < L ≤ U (the proof is immediate by case analysis on the three branches of the clamp). Proposition 9 (Downward-closedness) Rα is downward-closed with respect to ⪯: if en+1 ∈ Rα and e′ ⪯ en+1 , then e′ ∈ Rα . Proof See the proof in Appendix A.
4.5
Nested prediction regions
From a single calibration run, our method produces four related prediction regions with different trade-offs between validity, tightness, and computational cost (Proposition 12); see Figure 2(b) for an illustration. The inner rectangle is the most compact prediction region our method produces. It is contained inside Rα but carries no formal coverage guarantee. The exact region Rα is the tightest valid prediction set, but it is defined implicitly and costly to characterise on mixed cells. The staircase is a valid over-approximation of Rα given as an explicit union of at most 2d hyperrectangles, each computable at calibration time using the corner-cell quantiles. The outer rectangle is the most conservative: a single valid hyperrectangle, but the largest of the four. Outer rectangle (SCO) — valid. The outer rectangle is obtained by being doubly conservative: scores are computed relative to the lower base rectangle (Lk sides), which yields the largest possible calibration scores and therefore the largest possible quantile λ̂ub ; it then builds the confidence region using the largest possible base rectangle (Uk sides). The resulting hyperrectangle is guaranteed to contain Rα and hence achieves valid joint coverage, at the cost of being conservative. Assume Lk > 0 for all k. Formally, define upper-bound scores (k)
ei , k∈[d] Lk
σiub = max
λ̂ub = q-th order statistic of {σiub }ni=1 .
9
e(2)
e(2)
e(2)
e(2)
U2 U2
L2
U2 U2 λ̂ > 1
e(1)
U1
γ < γrect
λ̂ < 1
λ̂ = 1 U1
e(1)
γ = γrect
U1
0
e(1)
γ > γrect
L1 Inner Staircase
e(1)
U1 Exact Outer
(a) Effect of γ on Rα for Q fixed 1 − α (d = 2). The dashed (b) Prediction boundaries for d = 2. box is the rectangle k [0, Uk ], which grows with γ. The exact set Rα (red, zigzag) and its staircase approximation (green) Left (γ < γrect ): base rectangle too small, λ̂ > 1, Rα are both valid. The outer rectanextends beyond it. Centre (γ = γrect ): λ̂ = 1, Rα gle (blue) is conservative. The incoincides with the base rectangle. Right (γ > γrect ): ner rectangle (violet, dashed) lacks λ̂ < 1, Rα strictly inside the base rectangle. a coverage guarantee. Boundaries are slightly offset for visual clarity. Figure 2: Effect of γ on the prediction region shape, and the four nested region types (d = 2). Proposition 10 (Rectangular outer bound) Assume Uk > 0 for all k ∈ [d]. The exact prediction region is contained in the following closed-form hyperrectangle: Rα ⊆
d Y
0, λ̂ub Uk .
(8)
k=1
The region
Q
k [0, λ̂ub Uk ] achieves joint coverage ≥ 1 − α.
Proof See the proof in Appendix A. Algorithm 1 gives the calibration and prediction procedures for the outer rectangle. Inner rectangle (SCI) — no formal validity guarantee. The inner rectangle is doubly optimistic: scores are computed relative to the upper base rectangle (Uk sides), which yields the smallest possible calibration scores, giving the lowest possible quantile λ̂lb ; it then uses the smallest possible base rectangle (Lk sides) to build the confidence region. The result is tight but not guaranteed to be valid. Formally, define lower-bound scores (k)
ei , k∈[d] Uk
σilb = max
λ̂lb = q-th order statistic of {σilb }ni=1 .
Algorithm 2 gives the calibration and prediction procedures for the inner rectangle.
10
Algorithm 1: SCO: Scaling-Score Outer-Bound Predictor Input: Calibration errors {ei }ni=1 ⊂ Rd+ , miscoverage level α ∈ (0, 1), base-rectangle level γ ∈ (0, 1) Output: Half-widths h1 , . . . , hd ≥ 0 — Calibration phase (run once) — m′ ← clip(⌈γ(n + 1)⌉, 2, n); for k = 1, . . . , d do (k) (k) Lk ← e(m′ −1) , Uk ← e(m′ ) ;
end for i = 1, . . . , n do (k) σiub ← maxk∈[d] ei /Lk ; end q ← ⌈(1 − α)(n + 1)⌉; λ̂ub ← q-th order statistic of {σiub }ni=1 ; for k = 1, . . . , d do hk ← λ̂ub · Uk ; end — Prediction phase (per test point Xn+1 ) — Compute ŷ ← fˆ(Xn+1 ); Q return dk=1 [ŷ (k) − hk , ŷ (k) + hk ]
Exact region Rα (valid, tightest but implicit). The exact prediction region Rα is characterised cell by cell via Propositions 6 and 7: on corner cells and all-low-mid mixed cells it has an explicit form, but on mixed cells with Ihigh ̸= ∅ the boundary is only implicitly defined through (Imid ) λ̂(en+1 ), which depends on the test error. Although Rα is the tightest possible valid region, its boundary forms broken decision lines that are difficult to characterise explicitly and whose computation requires evaluating all calibration scores for every query. Staircase (SC2 ) — valid, union of at most 2d rectangles. The staircase is a valid overapproximation of Rα that replaces the implicit mixed-cell boundary by the upper-bounding cornercell quantile from Proposition 8. This turns every cell boundary into a fixed rectangle that is independent of the test error. The result is an explicit union of 2d hyperrectangles, one per subset S ⊆ [d], all computable at calibration time (Algorithm 3). The staircase is valid since it contains Rα , and it is tight in the sense that the 2d rectangles each touch the boundary of Rα at their corner points, as the following proposition shows. Proposition 11 (Staircase: union of hyperrectangles) Assume Uk > 0 for all k ∈ [d]. For each S ⊆ [d], define the corner point ( λ̂S̄,S Uk k ∈ S, ∗(k) eS = Uk k∈ / S, where S̄ = [d] \ S and λ̂S̄,S is the conformal quantile of the corner cell CS̄,∅,S . Then Rα ⊆
d [ Y
S⊆[d] k=1
11
∗(k)
[0, eS ].
Algorithm 2: SCI: Scaling-Score Inner-Bound Predictor Input: Calibration errors {ei }ni=1 ⊂ Rd+ , miscoverage level α ∈ (0, 1), base-rectangle level γ ∈ (0, 1) lb Output: Half-widths hlb 1 , . . . , hd ≥ 0 (not a valid prediction region) — Calibration phase — Compute m′ , Lk , Uk as in Algorithm 1; for i = 1, . . . , n do (k) σilb ← maxk∈[d] ei /Uk ; end q ← ⌈(1 − α)(n + 1)⌉; λ̂lb ← q-th order statistic of {σilb }ni=1 ; for k = 1, . . . , d do hlb k ← λ̂lb · Lk ; end
Proof See the proof in Appendix A. S Q ∗(k) The staircase S k [0, eS ] is valid (it contains Rα ), and all 2d corner rectangles are computable at calibration time: each λ̂S̄,S requires one pass over calibration scores. In practice, a simple post-processing step retains only the Pareto-optimal rectangles (those not dominated componentwise by any other), which often reduces the number of active rectangles well below 2d (empirically quantified in Section 5.2, Table 3). The four prediction regions are related by the following chain of containments. Proposition 12 (Rectangular sandwich) Assume Uk > 0 for all k ∈ [d]. Combining Proposition 10 with the cell-based analysis, whenever λ̂ub ≥ 1 (the regime that arises in practice), the four prediction regions form the nested chain d Y
0, λ̂lb Lk
⊆
k=1
|
Rα |{z}
⊆
exact region
{z
}
inner rectangle
d [ Y
∗(k) [0, eS ] S⊆[d] k=1
|
{z
staircase
}
⊆
d Y
0, λ̂ub Uk .
(9)
k=1
|
{z
outer rectangle
}
In the edge case λ̂ub < 1, the staircase need not be contained in the outer rectangle, but the inner sandwich d d Y [ Y ∗(k) 0, λ̂lb Lk ⊆ Rα ⊆ [0, eS ] (10) k=1
still holds, and Proposition 10 provides Rα ⊆
S⊆[d] k=1
Q
k [0, λ̂ub Uk ] as a separate over-approximation.
Proof See the proof in Appendix A. The inner rectangle is the most compact but offers no coverage guarantee; it may under-cover since it is strictly contained in Rα . Both the staircase and the outer rectangle are valid prediction regions (they contain Rα ). The full chain (9) holds in the typical regime λ̂ub ≥ 1. When λ̂ub < 1, which can occur at large γ, the outer rectangle shrinks below [0, U ], while the staircase still contains [0, U ] via its empty-subset corner e∗∅ = U , so the two regions become incomparable. These four objects are illustrated in Figure 2(b). 12
Algorithm 3: SC2 : Staircase Construction (Calibration Phase) Input: Calibration errors {ei }ni=1 ⊂ Rd+ , miscoverage level α ∈ (0, 1), base-rectangle level γ ∈ (0, 1) Output: Pareto-optimal staircase rectangles {[0, e∗S ]}S∈P Compute m′ , Uk as in Algorithm 1; q ← ⌈(1 − α)(n + 1)⌉; for each S ⊆ [d] do (k) Compute scores σiS̄,S ← maxk∈S ei /Uk for i = 1, . . . , n; λ̂S̄,S ← q-th order statistic of {σiS̄,S }ni=1 ; ( λ̂S̄,S Uk k ∈ S, ∗(k) Set corner point eS ← ; Uk k∈ / S. end Retain only Pareto-optimal rectangles: ∗(k) ∗(k) P ← {S ⊆ [d] : ∄ S ′ s.t. eS ′ ≤ eS ∀k, strictly for some k}; ∗ return {[0, eS ]}S∈P
4.6
Choice of γ
The hyperparameter γ ∈ (0, 1) determines the marginal quantile level at which the thresholds (k) Uk = e(m′ ) are set. Choosing γ too small or too large both degrade performance, though for different reasons. (k) When γ is small, the thresholds Uk = e(m′ ) sit at low marginal quantiles of the calibration errors. The base rectangle then under-represents the scale of large errors on each target: the outer quantile λ̂ must compensate by expanding the region uniformly, which dilutes the per-target scale adaptation that is the key advantage of the method. When γ is large, the thresholds Uk are set by the most extreme calibration errors, making them sensitive to outliers. Moreover, once γ > γrect the base rectangle already covers more than 1 − α of calibration points, so λ̂ < 1 and the prediction region shrinks to a rectangle strictly inside the base rectangle, losing the non-rectangular staircase structure. We define γrect as the value of γ for Q which the base rectangle covers exactly a fraction 1 − α of calibration errors. The base rectangle k [0, Uk ] has the equal-marginal-quantile property: each side sits at ECDF level γ. Because joint domination is strictly harder than marginal domination when targets are not perfectly co-monotone, γrect ≥ 1 − α in general, with equality only under perfect positive dependence. When γ = γrect , λ̂lb ≈ 1 and the inner rectangle approximates the empirical copula rectangle (with equality in the limit Lk → Uk , i.e. when the m′ -th and (m′ − 1)-th calibration order statistics coincide), while the outer rectangle provides a finitely valid conformal guarantee that the empirical copula lacks. We recommend targeting γ < γrect , so that the staircase retains its non-rectangular structure and λ̂ > 1 provides a meaningful scale signal. The choice γ = 1 − α is a natural default: it is tied to the coverage level, always satisfies γ ≤ γrect in practice, and our experiments confirm that it consistently matches or outperforms other values of γ across datasets and dimensions (see Section 5.2).
13
5
Experiments
5.1
Protocol
We use a random forest base model (100 trees) throughout. Each dataset is split into 60% train / 15% calibration / 25% test. Results are averaged over 20 random seeds. We evaluate at α ∈ {0.05, 0.10}. Baselines. • Max: max-aggregation score producing an axis-aligned hypercube. • CHR (Sampson and Chan, 2024): valid rectangular regions with per-target width adaptation, at the cost of splitting the calibration set in two. • TSCP (Fan and Sesia, 2025): Transductively Standardized Conformal Prediction, valid and model-agnostic with a single calibration set. Our method. We evaluate the SC2 prediction set under three choices of the hyperparameter γ, and SCO with γ = 1 − α: • SC2 , γ = 1 − α: a natural default tied to the coverage level. • SC2 , γ = 0.8: a fixed, coverage-level-agnostic choice, and a lower bound on all 1 − α levels considered in our experiments (0.80 ≤ 0.90, 0.95, 0.99). • SC2 , γrect : chosen on the calibration set so that exactly a fraction 1 − α of calibration points are dominated by the base rectangle (i.e. all their residuals fall within the base half-widths); see Section 4.6. • SCO, γ = 1 − α: the closed-form valid outer rectangle. Datasets. We evaluate on 29 real-world datasets spanning d ∈ {2, 3, 4, 6, 8, 12, 14, 16} targets, covering a wide range of sample sizes and output dimensions. See Appendix B for the full list with sample sizes and sources. Metrics. We evaluate each method on two metrics. Joint coverage: the fraction of test points whose full d-dimensional response falls inside the prediction region (target: at least 1 − α). Volume: the log10 volume of the prediction region, lower being better for a fixed coverage level.
5.2
Results
Coverage. Table 1 reports joint coverage at α = 0.10. All methods achieve valid coverage close to the nominal 1 − α = 0.90 level across all datasets, confirming the finite-sample guarantee. The one exception is SC2 with γ = 0.8 on Births 2 (d = 4): the random forest predicts one target near-perfectly, so more than 80% of calibration residuals on that target are exactly zero, yielding Uk = 0 and violating the assumption of Propositions 11 and 10. The same degenerate case arises for Solar Flare 1 at α = 0.05 (see Appendix C). Choosing γ ≥ 1 − α (all other SC2 and SCO variants) avoids this failure: if γ ≥ 1 − α and Uk = 0, then more than 1 − α of calibration residuals are zero, the transductive quantile remains finite, and coverage is maintained.
14
Dataset
d
Max
CHR TSCP
SC2
SCO
γ=1−α γ=0.8 γ=1−α γrect ANSUR-II Bio Births-1 Blog CalCOFI EDM ENB House Taxi Wage
2 2 2 2 2 2 2 2 2 2
0.903 0.899 0.899 0.901 0.899 0.905 0.915 0.900 0.901 0.899
0.897 0.900 0.897 0.901 0.900 0.929 0.926 0.902 0.902 0.902
0.904 0.900 0.897 0.902 0.900 0.929 0.908 0.900 0.901 0.903
0.908 0.900 0.896 0.901 0.901 0.936 0.917 0.900 0.901 0.902
0.906 0.900 0.895 0.901 0.900 0.927 0.910 0.899 0.901 0.899
0.906 0.900 0.895 0.901 0.901 0.922 0.911 0.900 0.901 0.901
0.907 0.900 0.896 0.902 0.901 0.906 0.922 0.901 0.902 0.902
Jura SCPF SF1 SF2 Student
3 3 3 3 3
0.918 0.902 0.932 0.905 0.909
0.937 0.912 0.946 0.889 0.904
0.912 0.909 0.932 0.906 0.916
0.914 0.913 0.935 0.872 0.917
0.911 0.909 0.854 0.868 0.911
0.910 0.908 0.925 0.870 0.910
0.945 0.911 0.962 0.907 0.930
Births-2 Households Stock
4 0.896 0.906 4 0.899 0.899 4 0.900 0.902
0.902 0.901 0.911
0.905 0.902 0.912
0.789 0.901 0.908
0.905 0.903 0.901 0.902 0.910 0.922
Air ATP1d ATP7d
6 0.899 0.899 6 0.902 0.908 6 0.915 0.929
0.900 0.890 0.910
0.898 0.908 0.921
0.899 0.896 0.915
0.898 0.902 0.899 0.916 0.910 0.923
RF1 RF2
8 0.899 0.902 8 0.903 0.903
0.899 0.904
0.900 0.905
0.900 0.904
0.900 0.899 0.905 0.906
OSales
12 0.902 0.907
0.891
0.901
0.905
0.899 0.918
WQ
14 0.887 0.898
0.899
0.902
0.895
0.897 0.895
OES10 OES97 SCM1d SCM20d
16 16 16 16
0.897 0.916 0.900 0.905
0.904 0.914 0.900 0.905
0.900 0.907 0.901 0.904
0.900 0.907 0.900 0.905
0.914 0.914 0.902 0.903
0.922 0.923 0.902 0.905
0.910 0.920 0.903 0.907
Table 1: Joint coverage on real-world datasets (random-forest base model, α = 0.10). Mean over 20 seeds. Target: ≥ 1 − α = 0.90.
15
Volume. Table 2 reports log10 volume. SC2 with γ = 1 − α achieves the best mean rank across the seven evaluated methods (mean rank 2.03) and ranks first on 12 of 29 datasets; the volume advantage over Max grows consistently with d. Among hyperrectangular methods, our SCO variant achieves a mean rank strictly better than all baselines (Max, CHR, TSCP), making it the strongest fully rectangular option. TSCP is competitive at low d but falls behind at high d. The SC2 γrect variant, which corresponds to the valid version of the copula-based construction (see Section 4.6), does not improve consistently over the default γ = 1−α, suggesting it offers no systematic advantage. Dataset
d
Max
CHR
TSCP
SC2
SCO γ=1−α
γ=0.8
γ=1−α
γrect
ANSUR-II Bio Births-1 Blog CalCOFI EDM ENB House Taxi Wage
2 3.40±0.01 2 6.08±0.00 2 6.42±0.00 2 3.33±0.01 2 1.59±0.00 2 0.47±0.03 2 1.74±0.02 2 11.17±0.00 2 −1.64±0.00 2 0.54±0.00
Jura SCPF SF1 SF2 Student
3 3 3 3 3
4.50±0.10 3.67±0.11 5.10±0.07 2.88±0.08 1.07±0.07 1.90±0.27 1.25±0.05 −0.72±0.43 3.46±0.03 3.48±0.05
3.38±0.06 2.78±0.05 1.04±0.08 0.11±0.09 3.53±0.03
3.35±0.07 2.83±0.06 1.44±0.13 0.51±0.04 3.52±0.03
3.33±0.08 2.94±0.07 1.28±0.19 0.45±0.05 3.49±0.03
3.31±0.07 2.76±0.06 1.11±0.12 0.47±0.04 3.49±0.03
3.66±0.07 2.79±0.05 1.41±0.09 0.03±0.15 3.61±0.04
Births-2 Households Stock
4 4 4
12.45±0.01 19.32±0.02 2.53±0.03
3.14±0.04 17.69±0.02 2.42±0.04
3.51±0.03 17.72±0.01 2.44±0.03
3.07±0.04 17.71±0.01 2.42±0.03
3.99±0.05 17.71±0.01 2.41±0.03
3.07±0.04 17.70±0.01 2.40±0.03
2.77±0.02 17.72±0.01 2.52±0.02
Air ATP1d ATP7d
6 6 6
9.47±0.02 15.88±0.17 15.24±0.19
4.59±0.03 4.65±0.01 15.97±0.18 15.52±0.15 15.85±0.25 15.34±0.23
4.57±0.01 15.79±0.17 15.27±0.28
4.64±0.02 4.57±0.01 15.65±0.15 15.61±0.16 15.14±0.26 15.10±0.27
4.65±0.02 15.97±0.12 15.69±0.17
RF1 RF2
8 8
5.56±0.04 5.86±0.05
3.54±0.07 3.94±0.08
4.51±0.14 4.87±0.16
3.47±0.05 3.96±0.06
3.62±0.04 4.08±0.06
3.46±0.05 3.95±0.06
3.54±0.05 4.03±0.05
OSales
12
58.14±0.40
52.73±0.50
51.94±0.35
51.90±0.37
52.08±0.36 51.69±0.35
53.36±0.28
WQ
14
12.76±0.04
13.61±0.14 12.47±0.04
13.34±0.07
13.28±0.06
13.21±0.07
12.59±0.04
OES10 OES97 SCM1d SCM20d
16 16 16 16
56.80±0.50 59.82±0.26 44.28±0.06 45.22±0.06
55.90±0.85 58.50±0.62 43.93±0.06 44.93±0.07
54.25±0.42 57.48±0.42 43.89±0.05 44.91±0.06
53.72±0.47 56.89±0.39 43.88±0.05 44.89±0.05
54.26±0.47 53.41±0.45 56.89±0.31 56.51±0.38 43.88±0.06 43.87±0.05 44.86±0.06 44.88±0.05
56.44±0.63 58.66±0.51 44.01±0.05 45.01±0.06
5.66
4.24
4.24
3.83
Mean rank
3.33±0.01 3.34±0.01 3.35±0.01 3.34±0.01 3.34±0.01 3.35±0.01 4.01±0.00 4.24±0.01 4.01±0.00 4.01±0.00 4.01±0.00 4.01±0.00 4.05±0.01 4.05±0.01 4.04±0.00 4.04±0.00 4.04±0.00 4.04±0.00 3.34±0.01 3.57±0.01 3.33±0.01 3.34±0.01 3.33±0.01 3.42±0.01 0.94±0.00 0.94±0.00 0.94±0.00 0.94±0.00 0.94±0.00 0.94±0.00 0.56±0.07 0.41±0.05 0.56±0.07 0.51±0.06 0.45±0.06 0.38±0.04 1.38±0.04 1.32±0.02 1.34±0.02 1.33±0.02 1.32±0.02 1.36±0.02 5.23±0.01 5.27±0.01 5.22±0.01 5.22±0.01 5.22±0.01 5.24±0.01 −1.63±0.01 −1.63±0.00 −1.64±0.00 −1.64±0.00 −1.64±0.00 −1.64±0.00 0.55±0.01 0.56±0.00 0.55±0.00 0.54±0.00 0.55±0.00 0.56±0.00
3.31
2.03
4.69
Table 2: log10 volume on real-world datasets (random-forest base model, α = 0.10). Mean over 20 seeds. Bold: lowest value per row. Staircase compactness. Table 3 reports the fraction of the 2d staircase corner rectangles that survive Pareto pruning (Algorithm 3), i.e. the proportion of the SC2 candidate set that is not 16
dominated and therefore stays part of the region. For the default γ = 1 − α, this fraction falls from a mean of 0.57 at d = 2 to 0.28 at d ∈ {3, 4}, 0.05 at d ∈ {6, 8}, and below 0.001 at d ≥ 12: although the number of candidate rectangles grows as 2d , pruning keeps the active set small in practice, so SC2 remains cheap to evaluate at the dimensions considered here. The γrect column illustrates the boundary case discussed in Section 4.6: when the base rectangle covers a fraction 1 − α of calibration points, λ̂S̄,S ≤ 1 for every S ̸= ∅ and the staircase collapses to its single empty-subset corner, 1/2d of the candidate set. This holds on 25 of the 29 datasets. It fails on ATP1d, ATP7d, OES10, and OES97, most severely on OES10 and OES97, where d = 16 but the calibration set has only n = 50 to 60 points, too few for any rectangle, even the loosest one, to jointly cover 1 − α of calibration points across all 16 targets. Up to several thousand of the 65,536 candidate rectangles then survive pruning. Dataset
SC2
d
γ=0.8 γ=1−α γrect ANSUR-II Bio Births-1 Blog CalCOFI EDM ENB House Taxi Wage
2 2 2 2 2 2 2 2 2 2
0.537 0.562 0.450 0.537 0.600 0.500 0.613 0.550 0.575 0.400
0.613 0.662 0.512 0.525 0.613 0.550 0.613 0.588 0.650 0.362
0.250 0.250 0.250 0.250 0.250 0.250 0.250 0.250 0.250 0.250
Jura SCPF SF1 SF2 Student
3 3 3 3 3
0.338 0.287 0.250 0.231 0.350
0.381 0.381 0.344 0.256 0.338
0.125 0.125 0.125 0.125 0.125
Births-2 Households Stock
4 4 4
0.078 0.228 0.212
0.081 0.062 0.231 0.062 0.169 0.062
Air ATP1d ATP7d
6 6 6
0.048 0.059 0.052
0.070 0.016 0.080 0.017 0.055 0.054
RF1 RF2
8 8
0.011 0.011
0.015 0.004 0.019 0.004
OSales
12
0.001
0.001 0.001
WQ
14
0.000
0.000 0.000
OES10 OES97 SCM1d SCM20d
16 16 16 16
0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000
0.004 0.007 0.000 0.000
Table 3: Fraction of the 2d staircase corner rectangles that survive Pareto pruning (random-forest base model, α = 0.10). Mean over 20 seeds. Lower values mean the staircase collapses to fewer active rectangles.
17
Runtime. Table 4 reports mean calibration time and formal complexity, for the methods compared in Table 2, on synthetic calibration errors at d = 8, ncalib = 300 (prediction is negligible for every method, well under one microsecond per point, and is omitted). Max and our SCO variant reduce to a handful of order-statistic computations and take a fraction of a millisecond; CHR has the same O(nd log n) complexity but a larger constant from its two-stage calibration split. The staircase variants with fixed γ (SC2 , γ ∈ {1 − α, 0.8}) enumerate up to 2d candidate corner rectangles and prune the dominated ones with a pairwise comparison, costing under 3 ms at d = 8; the calibration-adaptive γrect variant additionally scans O(n) candidate base-rectangle levels, each requiring an O(nd) pass over the calibration set, which raises its calibration time to about 7 ms. TSCP enumerates all partition cells in principle (O(nd )), but falls back to an O(dn log n) coordinate-wise search once d exceeds a small cutoff, which is the regime at d = 8 and explains why its time is comparable to the staircase variants despite the worst-case exponent. All calibration costs remain on the order of milliseconds or less at this scale, well within the cost of fitting the base regressor. Method Max CHR TSCP SCO (γ=1-α) SC2 (γ=0.8) SC2 (γ=1-α) SC2 (γrect )
Calib. time (ms) 0.034 0.096 4.654 0.059 2.756 2.770 7.291
Complexity O(nd + n log n) O(nd log n) O(min(nd , dn log n)) O(nd log n) O(2d n log n + 4d ) O(2d n log n + 4d ) O(n2 d + 2d n log n + 4d )
Table 4: Mean calibration time and formal complexity (n = ncalib , d = number of targets) for the methods compared in Table 2. Timings measured at d = 8, ncalib = 300, α = 0.1, synthetic errors, 10 seeds.
6
Conclusion
We introduced the scaling-score method for conformal multi-target regression. The method wraps any point predictor, requires only component-wise absolute residuals, and uses a single calibration set to produce four nested prediction regions: an inner rectangle (SCI, tight but without formal coverage guarantee), the exact set Rα (valid and tightest, characterised cell by cell), the staircase (SC2 , a valid union of at most 2d hyperrectangles, computable at calibration time), and the outer rectangle (SCO, a single closed-form valid hyperrectangle). All output types share a single calibration run and a single hyperparameter γ, for which γ = 1 − α is a natural and empirically reliable default. On the theoretical side, we gave an exact cell-by-cell characterisation of Rα , proved downward-closedness, and derived closed-form sandwich bounds. Experiments on 29 real-world datasets spanning d ∈ {2, . . . , 16} targets confirm valid joint coverage throughout. SC2 consistently matches or outperforms all baselines in volume, with the advantage over max-aggregation growing systematically with d, and the closed-form SCO outperforms all hyperrectangular baselines in mean rank. Future directions include extending the method to conditional coverage guarantees, to non-exchangeable settings such as time series, and to datadriven selection of γ.
18
References Annika Camehl, Dennis Fok, and Kathrin Gruber. On superlevel sets of conditional densities and multivariate quantile regression. Journal of Econometrics, 249:105807, May 2025. ISSN 0304-4076. doi: 10.1016/j.jeconom.2024.105807. Domagoj Cevid, Loris Michel, Jeffrey Näf, Peter Bühlmann, and Nicolai Meinshausen. Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1–79, 2022. Eustasio del Barrio, Alberto Gonzalez Sanz, and Marc Hallin. Nonparametric Multiple-Output Center-Outward Quantile Regression, April 2022. Victor Dheur, Matteo Fontana, Yorick Estievenart, Naomi Desobry, and Souhaib Ben Taieb. MultiOutput Conformal Regression: A Unified Comparative Study with New Conformity Scores, January 2025. Sašo Džeroski, Damjan Demšar, and Jasna Grbović. Predicting chemical parameters of river water quality from bioindicator data. Applied Intelligence, 13(1):7–17, 2000. Yunjie Fan and Matteo Sesia. Interpretable Multivariate Conformal Prediction with Fast Transductive Standardization, December 2025. Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated multiple-output quantile regression with representation learning. Journal of Machine Learning Research, 24(24):1–48, 2023. Pierre Goovaerts. Geostatistics for Natural Resources Evaluation. Oxford university press, 1997. Rafael Izbicki, Gilson T. Shimizu, and R. Stern. CD-split and HPD-split: Efficient conformal regions in high dimensions. Journal of Machine Learning Research, 23:87:1–87:32, 2020. Aram Karalič and Ivan Bratko. First Order Regression. Machine Learning, 26(2-3):147–176, February 1997. ISSN 0885-6125, 1573-0565. doi: 10.1023/A:1007365207130. Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman. DistributionFree Predictive Inference for Regression. Journal of the American Statistical Association, 113 (523):1094–1111, July 2018. ISSN 0162-1459. doi: 10.1080/01621459.2017.1307116. Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Copula-based conformal prediction for Multi-Target Regression. arXiv:2101.12002 [cs, stat], January 2021. Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In European Conference on Machine Learning, pages 345–356. Springer, 2002. Yaniv Romano, Evan Patterson, and Emmanuel J. Candès. Conformalized Quantile Regression. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 3538–3548, 2019. Max Sampson and Kung-Sik Chan. Conformal Multi-Target Hyperrectangles. Statistical Analysis and Data Mining: The ASA Data Science Journal, 17(5):e11710, October 2024. ISSN 1932-1864, 1932-1872. doi: 10.1002/sam.11710. 19
Athanasios Tsanas and Angeliki Xifara. Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools. Energy and buildings, 49:560–567, 2012. Grigorios Tsoumakas, Eleftherios Spyromitros-Xioufis, Jozef Vilcek, and Ioannis Vlahavas. Mulan: A java library for multi-label learning. The Journal of Machine Learning Research, 12:2411–2414, 2011. Vladimir Vovk, Alex Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. Springer Science & Business Media, 2005. Zhendong Wang, Ruijiang Gao, Mingzhang Yin, Mingyuan Zhou, and David M. Blei. Probabilistic Conformal Prediction Using Conditional Random Samples, June 2022.
20
A
Proofs
Proof [Proof of Theorem 3] Both bn+1 and the scores σi (en+1 ) depend only on the multiset {e1 , . . . , en , en+1 } and are symmetric in the indices {1, . . . , n + 1}. Under exchangeability, the rank of σn+1 among {σ1 , . . . , σn+1 } is uniformly distributed on {1, . . . , n + 1}. Hence P(σn+1 (en+1 ) ≤ λ̂(en+1 )) ≥ 1 − α, which is exactly P(en+1 ∈ Rα ) ≥ 1 − α. (k)
Proof [Proof of Lemma 5] Substitute the cell formula for bn+1 into (1): each coordinate contributes (k)
(k)
(k)
(k)
(k)
(k)
(k)
(k)
en+1 /Lk , en+1 /en+1 = 1, or en+1 /Uk to σn+1 , and ei /Lk , ei /en+1 , or ei /Uk to σi . Non-increase (k)
(k)
(k)
of ei /en+1 in en+1 is immediate. Proof [Proof of Proposition 6] The intersection requires, for each coordinate: (k)
(k)
(k)
(k)
(k)
• k ∈ Ilow : en+1 ≤ Lk (cell) and en+1 ≤ λ̂Lk (rectangle), giving en+1 ≤ min(λ̂, 1)Lk . (k)
• k ∈ Ihigh : en+1 > Uk (cell) and en+1 ≤ λ̂Uk (rectangle), giving en+1 ∈ (Uk , λ̂Uk ], non-empty iff λ̂ > 1.
(k)
(k)
(k)
Proof [Proof of Proposition 8] Each σi (en+1 ) = max(Di , maxk∈Imid ei /en+1 ), where Di = max(maxk∈Ilow ei /Lk , (k)
(k)
(k)
is constant on the cell. Since ei /en+1 is decreasing in en+1 , σi is non-increasing in each mid coordi(k)
nate, and so is the q-th order statistic. The boundary values follow by substituting en+1 = Lk (resp. (k)
(k)
Uk ) in the mid terms, which gives ei /Lk (resp. ei /Uk ): the scores then match those of the corner cell (Ilow ∪ Imid , ∅, Ihigh ) (resp. (Ilow , ∅, Imid ∪ Ihigh )). Since λ̂ is the quantile of these scores, it too matches the adjacent corner-cell value at each face, so the boundary of Rα is continuous across cell boundaries. Proof [Proof of Proposition 7] The test score formula follows from Lemma 5: on the cell, the low (k) (k) terms en+1 /Lk ≤ 1 and the mid terms equal 1, so σn+1 = max(1, maxk∈Ihigh en+1 /Uk ). (k)
For item 1: σn+1 = 1, so σn+1 ≤ λ̂ iff λ̂ ≥ 1. The low coordinates satisfy en+1 ≤ Lk automatically. (k) (k) For item 2: σn+1 ≤ λ̂ iff en+1 ≤ λ̂ Uk for all k ∈ Ihigh . Combined with en+1 > Uk , this requires λ̂ > 1. Proof [Proof of Proposition 9] Base vector. From (4), b(k) (e) = clamp(e(k) , Lk , Uk ) is nondecreasing in e(k) , so e′ ⪯ en+1 ⇒ b(e′ ) ⪯ b(en+1 ). Test score. Since x/ clamp(x, L, U ) is non-decreasing (each branch x/L, 1, x/U is non-decreasing), (k) for each k: e′(k) /b(k) (e′ ) ≤ en+1 /b(k) (en+1 ). Taking the max: σn+1 (e′ ) ≤ σn+1 (en+1 ). Calibration quantile. Since b(e′ ) ⪯ b(en+1 ), each calibration score satisfies σi (e′ ) ≥ σi (en+1 ) for all i ≤ n, so λ̂(b(e′ )) ≥ λ̂(b(en+1 )). Conclusion. σn+1 (e′ ) ≤ σn+1 (en+1 ) ≤ λ̂(b(en+1 )) ≤ λ̂(b(e′ )), where the middle inequality uses en+1 ∈ Rα . Hence e′ ∈ Rα . (k)
(k)
(k)
Proof [Proof of Proposition 10] Since bn+1 ≥ Lk , for each i ≤ n: σi (en+1 ) = maxk ei /bn+1 ≤ (k)
maxk ei /Lk = σiub . Hence λ̂(bn+1 ) ≤ λ̂ub . For any en+1 ∈ Rα , by definition σn+1 (en+1 ) ≤ λ̂(en+1 ), (k) (k) so: σn+1 (en+1 ) ≤ λ̂(bn+1 ) ≤ λ̂ub . By equation (2): en+1 ≤ λ̂ub bn+1 ≤ λ̂ub Uk for all k. 21
(k)
(k)
Proof [Proof of Proposition 12] Inner ⊆ exact. Let en+1 ≤ λ̂lb Lk for all k. Since bn+1 ≥ Lk : (k)
σn+1 (en+1 ) = max k
en+1 (k)
bn+1
(k)
≤ max
(k)
k
en+1 ≤ λ̂lb . Lk
(k)
Since bn+1 ≤ Uk , for each i ≤ n: σi (en+1 ) ≥ maxk ei /Uk = σilb . Taking the q-th order statistic: λ̂(bn+1 ) ≥ λ̂lb . Combining: σn+1 (en+1 ) ≤ λ̂lb ≤ λ̂(bn+1 ), i.e. en+1 ∈ Rα . Exact ⊆ staircase. This is Proposition 11. This proves (10). ∗(k) Staircase ⊆ outer rectangle (under λ̂ub ≥ 1). Each corner rectangle satisfies eS ≤ λ̂ub Uk : for k ∈ S, λ̂S̄,S ≤ λ̂ub since σiub ≥ σiS̄,S ; for k ∈ / S, Uk ≤ λ̂ub Uk because λ̂ub ≥ 1. Hence the union Q is contained in k [0, λ̂ub Uk ], yielding (9) together with Proposition 10. The hypothesis λ̂ub ≥ 1 ∗(k) is necessary: if λ̂ub < 1, then λ̂ub Uk < Uk = e∅ in every coordinate, so the empty-subset corner e∗∅ = U of the staircase already lies strictly outside the outer rectangle. Proof [Proof of Proposition 11] Let e ∈ Rα and set S = {k : e(k) > Uk }. For k ∈ / S: e(k) ≤ ∗(k) Uk = eS by definition. For k ∈ S: e belongs to a cell with Ihigh = S and Ilow ∪ Imid = S̄. When Imid = ∅, Proposition 6 gives e(k) ≤ λ̂S̄,S Uk . When Imid ̸= ∅, Proposition 8 gives λ̂(e(Imid ) ) ≤ λ̂S̄,S ; the exact-region formula for mixed cells then gives e(k) ≤ λ̂(e(Imid ) )Uk ≤ λ̂S̄,S Uk . In either case ∗(k)
e(k) ≤ λ̂S̄,S Uk = eS . Hence e ⪯ e∗S and e lies in the box indexed by S.
B
Dataset Details Table 5: Real-world datasets used in the experiments. Split: 60 % train / 15 % calibration / 25 % test averaged over 20 random seeds. Dataset
d
n
ntrain
ncalib
ntest
EDM ENB ANSUR II Births 1 Wage Blog Feedback House Prices NYC Taxi Bio CalCOFI
2 2 2 2 2 2 2 2 2 2
154 768 1,986 10,000 10,000 20,000 21,613 30,000 45,730 50,000
92 460 1,191 6,000 6,000 12,000 12,967 18,000 27,438 30,000
23 115 297 1,500 1,500 3,000 3,241 4,500 6,859 7,500
39 193 498 2,500 2,500 5,000 5,405 7,500 11,433 12,500
Solar Flare 1 Jura Student Perf. Solar Flare 2
3 3 3 3
323 359 395 1,066
193 215 237 639
48 53 59 159
82 91 99 268
Source Karalič and Bratko (1997) Tsanas and Xifara (2012) del Barrio et al. (2022) Cevid et al. (2022) Cevid et al. (2022) Feldman et al. (2023) Feldman et al. (2023) Wang et al. (2022) Feldman et al. (2023) del Barrio et al. (2022) UCI ML Repository Goovaerts (1997) UCI ML Repository UCI ML Repository (continued on next page)
22
C
Dataset
d
n
ntrain
ncalib
ntest
Source
SCPF
3
1,137
682
170
285
Kaggle
Stock Households Births 2
4 4 4
950 7,207 10,000
570 4,324 6,000
142 1,081 1,500
238 1,802 2,500
Wang et al. (2022) Camehl et al. (2025) Cevid et al. (2022)
ATP 7D ATP 1D Air Quality
6 6 6
296 337 10,000
177 202 6,000
44 50 1,500
75 85 2,500
Tsoumakas et al. (2011) Tsoumakas et al. (2011) Cevid et al. (2022)
River Flow 2 River Flow 1
8 8
7,679 9,005
4,607 5,403
1,151 1,350
1,921 2,252
Tsoumakas et al. (2011) Tsoumakas et al. (2011)
Online Sales
12
556
333
83
140
Kaggle
Water Quality
14
1,060
636
159
265
Džeroski et al. (2000)
OES 1997 OES 2010 SCM 20-day SCM 1-day
16 16 16 16
334 403 8,966 9,803
200 241 5,379 5,881
50 60 1,344 1,470
84 102 2,243 2,452
Tsoumakas et al. (2011) Tsoumakas et al. (2011) Tsoumakas et al. (2011) Tsoumakas et al. (2011)
Full Experimental Results
Each figure row below corresponds to a group of datasets sharing the same number of targets d. The title of each panel follows the format Dataset name (d targets, ntrain train, ncalib calib, ntest test). The left plot shows joint coverage (dashed line: nominal 1 − α) and the right plot shows log10 volume, both as a function of the method. Four groups of figures are included: random forest with α = 0.10, random forest with α = 0.05, linear model with α = 0.10, and linear model with α = 0.05. Within each panel, methods are listed in decreasing order of performance, so the method at the top is the best performer for that particular dataset.
Random forest, α = 0.10 ANSUR-II (d=2) 1191/297/498
Bio (d=2) 27438/6859/11433
CHR
Births-1 (d=2) 6000/1500/2500
CHR
TSCP SC2 =1
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 =1
SC2 rect
SC2 =0.8
SC2 rect
CHR
SC2 rect
TSCP
TSCP
Max
Max 3.35 log10 vol.
3.40
0.7
0.8 0.9 Joint cov.
1.0
Max 4
5 log10 vol.
Blog (d=2) 12000/3000/5000
6
0.7
0.8 0.9 Joint cov.
1.0
Max
SC2 rect
SC2 =0.8
SC2 =0.8
CHR
1.0
0.8 0.9 Joint cov.
1.0
TSCP SC2 =1 Max SC2 =0.8
Max 0.8 0.9 Joint cov.
0.7
EDM (d=2) 92/23/39
SC2 =1
TSCP
6
SC2 rect
TSCP
SC2 rect 0.7
5 log10 vol.
CalCOFI (d=2) 30000/7500/12500 CHR
SC2 =1
3.4 3.5 log10 vol.
4
CHR 1.0
1.2 1.4 log10 vol.
23
1.6 0.7
0.8 0.9 Joint cov.
1.0
0.4
0.5 log10 vol.
0.6
0.7
0.8 0.9 Joint cov.
1.0
ENB (d=2) 460/115/193
House (d=2) 12967/3241/5405
SC2 =1 TSCP
Taxi (d=2) 18000/4500/7500
SC2 =0.8
SC2 =0.8
SC2 =1
SC2 =1
SC2 =0.8
CHR
Max
SC2 rect
SC2 rect
SC2 rect
CHR
TSCP
CHR
Max
Max
TSCP
1.4
1.6 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
6
8 10 log10 vol.
Wage (d=2) 6000/1500/2500
0.7
0.8 0.9 Joint cov.
1.0
1.640 1.635 1.630 log10 vol.
Jura (d=3) 215/53/91
Max SC2 =0.8
SC2 =0.8
TSCP
TSCP
SC2 rect
CHR
SC2 rect
CHR
CHR
SC2 =0.8
Max 0.53 0.54 0.55 0.56 log10 vol.
0.7
0.8 0.9 Joint cov.
Max
1.0
3.5
4.0 log10 vol.
SF1 (d=3) 193/48/82
4.5
0.7
0.8 0.9 Joint cov.
1.0
CHR
Max
SC2 rect
SC2 =0.8
SC2 =0.8
SC2 rect
SC2 =1
CHR 2.0
0.7
0.8 0.9 Joint cov.
1.0
TSCP SC2 rect 1
0 log10 vol.
1
0.7
0.8 0.9 Joint cov.
1.0
3.5 3.6 log10 vol.
Households (d=4) 4324/1081/1802
CHR
SC2 rect
TSCP
SC2 =0.8
TSCP
SC2 rect
Max
Max 0.8 0.9 Joint cov.
Max
1.0
18.0
Air (d=6) 6000/1500/2500
18.5 19.0 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
2.4
ATP1d (d=6) 202/50/85
CHR
SC2 =1
Max
TSCP
Max
TSCP
SC2 rect
CHR
SC2 rect
Max
SC2 rect 0.8 0.9 Joint cov.
1.0
0.7
0.8 0.9 Joint cov.
1.0
15.0
RF2 (d=8) 4607/1151/1921
SC2 rect
SC2 =0.8
SC2 =0.8
CHR
TSCP
TSCP
SC2 rect
5
0.7
0.8 0.9 Joint cov.
1.0
0.8 0.9 Joint cov.
1.0
SC2 =0.8
Max log10 vol.
0.7
TSCP
SC2 =1
4
16.0
SC2 =1
CHR
Max
15.5 log10 vol.
OSales (d=12) 333/83/140
CHR
SC2 rect
1.0
CHR 15.50 15.75 16.00 log10 vol.
RF1 (d=8) 5403/1350/2252 SC2 =1
0.8 0.9 Joint cov.
SC2 =0.8
SC2 =0.8
0.7
0.7
SC2 =1
SC2 =0.8
6 8 log10 vol.
2.5 log10 vol.
ATP7d (d=6) 177/44/75
TSCP
SC2 =1
1.0
SC2 =0.8
TSCP
0.7
0.8 0.9 Joint cov.
Stock (d=4) 570/142/238
SC2 =0.8
10
0.7
SC2 =1
SC2 =1
log10 vol.
1.0
SC2 =0.8
CHR
5
0.8 0.9 Joint cov.
CHR
CHR
SC2 =1
0.7
SC2 =1
Births-2 (d=4) 6000/1500/2500 SC2 rect
5
Student (d=3) 237/59/99
Max 1.5 log10 vol.
4 log10 vol.
Max
TSCP
SC2 =1
1.0
3
SF2 (d=3) 639/159/268
TSCP
1.0
SC2 =1
SC2 rect TSCP
0.8 0.9 Joint cov.
SCPF (d=3) 682/170/285
SC2 =1
SC2 =1
0.7
Max 4
5 log10 vol.
24
6 0.7
0.8 0.9 Joint cov.
1.0
52
54 56 log10 vol.
58
0.7
0.8 0.9 Joint cov.
1.0
WQ (d=14) 636/159/265
OES10 (d=16) 241/60/102
TSCP
OES97 (d=16) 200/50/84
SC2 =1
SC2 =1
SC2 rect
TSCP
SC2 =0.8
Max
SC2 =0.8
TSCP
SC2 =1
CHR
CHR
SC2 =0.8
SC2 rect
SC2 rect
CHR
Max 12.5
13.0 13.5 log10 vol.
0.7
0.8 0.9 Joint cov.
Max
1.0
54
SCM1d (d=16) 5881/1470/2452
56 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
56
58 log10 vol.
60
0.7
0.8 0.9 Joint cov.
1.0
SCM20d (d=16) 5379/1344/2243
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 =1
TSCP
TSCP
CHR
CHR
SC2 rect
SC2 rect
Max
Max 43.8
44.0 44.2 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
44.8
45.0 45.2 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
Random forest, α = 0.05 ANSUR-II (d=2) 1191/297/498
Bio (d=2) 27438/6859/11433
SC2 =0.8 TSCP
Births-1 (d=2) 6000/1500/2500
SC2 =0.8
SC2 =1
SC2 =1
SC2 =0.8
CHR
CHR
CHR
SC2 rect
SC2 rect
SC2 rect
TSCP
TSCP
Max
Max
SC2 =1
3.45
3.50 log10 vol.
3.550.7
0.8
0.9 Joint cov.
Max
1.0
5 6 log10 vol.
Blog (d=2) 12000/3000/5000
0.7
0.8
0.9 Joint cov.
1.0
5 6 log10 vol.
CalCOFI (d=2) 30000/7500/12500
SC2 =1 CHR
SC2 =0.8 CHR
Max
TSCP
TSCP
SC2 rect
SC2 =0.8
Max 0.8
0.9 Joint cov.
1.0
CHR 1.2
ENB (d=2) 460/115/193
1.4 1.6 log10 vol.
1.80.7
0.8
0.9 Joint cov.
1.0
0.5
House (d=2) 12967/3241/5405
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 =1
CHR
SC2 rect
TSCP
TSCP
SC2 rect
0.8
0.9 Joint cov.
1.0
8 10 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
1.420 1.415 1.410 0.7 log10 vol.
Jura (d=3) 215/53/91 SC2 =1
TSCP
SC2 =0.8
SC2 =0.8
SC2 rect
SC2 rect
CHR
SC2 =0.8
Max
Max 0.8
0.9 Joint cov.
1.0
1.0
SC2 =1
TSCP
0.7
0.9 Joint cov.
SCPF (d=3) 682/170/285
CHR
0.68 0.70 log10 vol.
0.8
CHR
SC2 rect
TSCP
1.0
CHR 6
Wage (d=2) 6000/1500/2500 SC2 =1
0.9 Joint cov.
SC2 =1
Max 2.00.7
0.8
SC2 =0.8
SC2 rect 1.8 log10 vol.
0.7 0.7
Taxi (d=2) 18000/4500/7500
CHR
1.6
0.6 log10 vol.
Max
TSCP
Max
1.0
SC2 =1
Max 3.80 3.85 3.90 3.95 0.7 log10 vol.
0.9 Joint cov.
SC2 rect
TSCP SC2 =0.8
0.8
EDM (d=2) 92/23/39
SC2 =1
SC2 rect
0.7
Max 4
5 log10 vol.
25
0.7
0.8
0.9 Joint cov.
1.0
4
5 log10 vol.
6 0.7
0.8
0.9 Joint cov.
1.0
SF1 (d=3) 193/48/82
SF2 (d=3) 639/159/268
Max
Student (d=3) 237/59/99
CHR
SC2 rect
SC2 =0.8 TSCP
SC2 rect
TSCP
SC2 =1
SC2 =1 SC2 =0.8 CHR
SC2 =1
TSCP
Max
SC2 =0.8
SC2 rect
Max 2.0 log10 vol.
2.5 0.7
0.8
0.9 Joint cov.
1.0
CHR 0
1 log10 vol.
Births-2 (d=4) 6000/1500/2500
0.7
0.8
0.9 Joint cov.
1.0
3.8 3.9 log10 vol.
Households (d=4) 4324/1081/1802
SC2 =1
SC2 =1
TSCP
TSCP
TSCP
CHR
SC2 =0.8
SC2 rect
Max
Max
Max
SC2 rect
0.7
0.8
0.9 Joint cov.
1.0
18.5 19.0 19.5 log10 vol.
Air (d=6) 6000/1500/2500
0.7
0.8
0.9 Joint cov.
1.0
2.8
ATP1d (d=6) 202/50/85
SC2 =1
SC2 =0.8
CHR
TSCP
SC2 =1 SC2 =0.8
SC2 rect
Max SC2 rect
Max
Max
CHR
CHR
0.8
0.9 Joint cov.
1.0
16.5
RF1 (d=8) 5403/1350/2252
17.0 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
16.2
CHR
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 rect
SC2 rect
SC2 rect
TSCP
TSCP
CHR
Max
Max 70.7
0.8
0.9 Joint cov.
6 log10 vol.
WQ (d=14) 636/159/265
7
0.7
0.8
0.9 Joint cov.
1.0
55.0
OES10 (d=16) 241/60/102 SC2 =0.8
SC2 =0.8
Max
SC2 =1
SC2 =1
TSCP
TSCP
SC2 rect
SC2 rect
CHR
Max
CHR
SC2 =0.8
CHR 13.0
13.5 14.0 0.7 log10 vol.
0.8
0.9 Joint cov.
1.0
SCM1d (d=16) 5881/1470/2452
60.0 62.5 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
SCM20d (d=16) 5379/1344/2243 SC2 =0.8
SC2 =1
SC2 =1
TSCP
TSCP
CHR
CHR
SC2 rect
SC2 rect
Max
Max 45.6
45.8 46.0 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
1.0
57.5 60.0 0.7 log10 vol.
0.8
0.9 Joint cov.
1.0
Max 57.5
SC2 =0.8
0.9 Joint cov.
OES97 (d=16) 200/50/84
SC2 rect
SC2 =1
0.8
Max
1.0
TSCP
1.0
OSales (d=12) 333/83/140
SC2 =1
6 log10 vol.
0.9 Joint cov.
TSCP
SC2 =0.8
5
16.4 16.6 16.8 0.7 log10 vol.
RF2 (d=8) 4607/1151/1921 CHR
SC2 =1
0.8
SC2 rect
SC2 =0.8 100.7
3.2 0.7
ATP7d (d=6) 177/44/75
TSCP
8 log10 vol.
3.0 log10 vol.
TSCP
SC2 =1
6
1.0
SC2 =0.8
SC2 =0.8
10 log10 vol.
0.9 Joint cov.
SC2 =1
CHR
5
0.8
Stock (d=4) 570/142/238
CHR
SC2 rect
0.7
46.4
46.6 46.8 0.7 log10 vol.
26
0.8
0.9 Joint cov.
1.0
60
62 log10 vol.
64
0.7
0.8
0.9 Joint cov.
1.0
Linear model, α = 0.10 ANSUR-II (d=2) 1191/297/498
Bio (d=2) 27438/6859/11433
CHR SC2 =0.8
Births-1 (d=2) 6000/1500/2500
SC2 =0.8
SC2 =0.8
SC2 =1
SC2 =1
SC2 rect
CHR
TSCP
CHR
SC2 rect
SC2 rect
TSCP
TSCP
Max
Max
SC2 =1
3.24 3.26 3.28 3.30 log10 vol.
0.7
0.8 0.9 Joint cov.
Max
1.0
5 6 log10 vol.
Blog (d=2) 12000/3000/5000
0.7
0.8 0.9 Joint cov.
1.0
4
CalCOFI (d=2) 30000/7500/12500
Max SC2 =1
SC2 =1
SC2 =1
CHR
SC2 =0.8
SC2 =0.8
SC2 rect
TSCP
TSCP
TSCP
Max 0.7
0.8 0.9 Joint cov.
1.0
1.4 1.6 log10 vol.
1.8 0.7
0.8 0.9 Joint cov.
1.0
0.6 log10 vol.
House (d=2) 12967/3241/5405
Max
CHR
CHR
SC2 rect
TSCP
SC2 rect
CHR
Max 0.8 0.9 Joint cov.
1.0
TSCP 6
Wage (d=2) 6000/1500/2500
8 10 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
1.67 1.66 log10 vol.
Jura (d=3) 215/53/91 SC2 =0.8
SC2 =1
SC2 =1
SC2 =1
SC2 =0.8
CHR
TSCP
TSCP
Max
CHR
SC2 rect
TSCP
SC2 rect
CHR
SC2 rect
Max 0.54
0.7
0.8 0.9 Joint cov.
3.5 4.0 log10 vol.
SF1 (d=3) 193/48/82
0.7
0.8 0.9 Joint cov.
1.0
3
SF2 (d=3) 639/159/268 CHR
CHR
Max
SC2 rect
TSCP
SC2 rect
TSCP
Max
CHR
Max
SC2 rect
1.5 0.7
0.8 0.9 Joint cov.
1.0
0 log10 vol.
Births-2 (d=4) 6000/1500/2500
1
0.7
0.8 0.9 Joint cov.
1.0
3.55 3.60 3.65 3.70 log10 vol.
Households (d=4) 4324/1081/1802
SC2 =1
TSCP
CHR
TSCP
SC2 rect
SC2 rect
Max log10 vol.
10
0.7
0.8 0.9 Joint cov.
1.0
0.7
0.8 0.9 Joint cov.
1.0
TSCP
SC2 =0.8
5
1.0
SC2 =0.8
SC2 =1
Max
0.8 0.9 Joint cov.
SC2 =1
SC2 =0.8
CHR
0.7
Stock (d=4) 570/142/238
CHR
SC2 rect
5
SC2 =0.8
SC2 =0.8
1.0 log10 vol.
4 log10 vol.
SC2 =1
TSCP
0.5
1.0
Student (d=3) 237/59/99
SC2 =1
SC2 =1
0.8 0.9 Joint cov.
Max
1.0
SC2 =0.8
0.7
SCPF (d=3) 682/170/285
SC2 =0.8
0.52 log10 vol.
1.0
Max
SC2 rect
0.7
0.8 0.9 Joint cov.
SC2 =0.8
SC2 =0.8
2.30 log10 vol.
0.7
SC2 =1
SC2 =1
2.25
0.7
Taxi (d=2) 18000/4500/7500
SC2 =0.8
TSCP
1.0
CHR 1.2
ENB (d=2) 460/115/193 SC2 =1
0.8 0.9 Joint cov.
SC2 rect
SC2 rect
4.0
0.7
Max
SC2 =0.8
3.9 log10 vol.
6
EDM (d=2) 92/23/39
CHR
3.8
5 log10 vol.
Max 18.0
18.5 19.0 log10 vol.
27
0.7
0.8 0.9 Joint cov.
1.0
4.0 4.2 log10 vol.
0.7
0.8 0.9 Joint cov.
1.0
Air (d=6) 6000/1500/2500
ATP1d (d=6) 202/50/85
SC2 =1
ATP7d (d=6) 177/44/75
SC2 =0.8
SC2 =0.8
TSCP
SC2 =0.8 CHR
TSCP
SC2 =1
SC2 =1
TSCP
Max
Max
SC2 rect
SC2 rect
CHR
Max
CHR 6
8 log10 vol.
0.7
0.8 0.9 Joint cov.
SC2 rect
1.0
19.5 20.0 log10 vol.
RF1 (d=8) 5403/1350/2252
20.5 0.7
0.8 0.9 Joint cov.
1.0
17.75
RF2 (d=8) 4607/1151/1921
SC2 =1 CHR
CHR SC2 =0.8
SC2 =0.8
SC2 rect
CHR
TSCP
TSCP
SC2 rect
TSCP
Max
1.0
10 log10 vol.
WQ (d=14) 636/159/265
11
0.7
0.8 0.9 Joint cov.
1.0
62.5 65.0 log10 vol.
OES10 (d=16) 241/60/102 SC2 =0.8
SC2 =0.8
Max
SC2 =1
SC2 =1
SC2 rect
TSCP
SC2 =0.8
CHR
CHR
SC2 =1
SC2 rect
SC2 rect
13.5 14.0 log10 vol.
0.7
0.8 0.9 Joint cov.
0.7
0.8 0.9 Joint cov.
1.0
TSCP
Max 13.0
67.5
OES97 (d=16) 200/50/84
TSCP
CHR
1.0
SC2 =1
Max 0.8 0.9 Joint cov.
0.8 0.9 Joint cov.
SC2 =0.8
SC2 rect
Max
0.7
OSales (d=12) 333/83/140
SC2 =1
11.0 11.5 12.0 12.5 0.7 log10 vol.
18.00 18.25 log10 vol.
Max
1.0
58
SCM1d (d=16) 5881/1470/2452
log10 vol.
60
0.7
0.8 0.9 Joint cov.
1.0
60
62 log10 vol.
64
0.7
0.8 0.9 Joint cov.
1.0
SCM20d (d=16) 5379/1344/2243 TSCP
SC2 =1
CHR
SC2 =0.8 TSCP
SC2 =1
CHR
SC2 =0.8
SC2 rect
SC2 rect
Max
Max 45.0
45.2 log10 vol.
45.4 0.7
0.8 0.9 Joint cov.
1.0
47.4
47.6 47.8 log10 vol.
48.0 0.7
0.8 0.9 Joint cov.
1.0
Linear model, α = 0.05 ANSUR-II (d=2) 1191/297/498
Bio (d=2) 27438/6859/11433
Births-1 (d=2) 6000/1500/2500
CHR
SC2 =1 TSCP
SC2 =1
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 =0.8
CHR
CHR
SC2 rect
SC2 rect
SC2 rect
TSCP
TSCP
Max
Max 3.38 3.40 3.42 3.44 0.7 log10 vol.
0.8
0.9 Joint cov.
Max
1.0
5
Blog (d=2) 12000/3000/5000
6 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
5 6 log10 vol.
CalCOFI (d=2) 30000/7500/12500
CHR
SC2 =1
SC2 =0.8
TSCP
SC2 =0.8
SC2 rect
SC2 rect
SC2 rect
TSCP
SC2 =0.8
Max 0.8
0.9 Joint cov.
1.0
1.0
SC2 =1
TSCP
4.2 4.3 0.7 log10 vol.
0.9 Joint cov.
Max
Max
4.1
0.8
EDM (d=2) 92/23/39
CHR
SC2 =1
0.7
CHR 1.4
1.6 1.8 log10 vol.
28
2.0 0.7
0.8
0.9 Joint cov.
1.0
0.65
0.70 0.75 log10 vol.
0.80 0.7
0.8
0.9 Joint cov.
1.0
ENB (d=2) 460/115/193
House (d=2) 12967/3241/5405
Taxi (d=2) 18000/4500/7500
CHR
SC2 =1
Max
Max
SC2 =1
SC2 =0.8
TSCP
SC2 =0.8
SC2 =1
SC2 =0.8
TSCP
SC2 rect
SC2 rect
CHR
CHR
Max
TSCP
2.40
2.45 2.50 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
SC2 rect
6
8 10 log10 vol.
Wage (d=2) 6000/1500/2500
120.7
0.8
0.9 Joint cov.
1.0
1.46
Jura (d=3) 215/53/91
SC2 =1
TSCP
CHR
SC2 =0.8
TSCP
SC2 =0.8
SC2 rect
SC2 rect
SC2 rect
CHR
SC2 =0.8
Max
Max 0.8
0.9 Joint cov.
1.0
4.0
4.5 5.0 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
4
SF2 (d=3) 639/159/268 CHR SC2 rect
Max
SC2 =0.8
TSCP
SC2 rect
0.8
0.9 Joint cov.
0.5
Births-2 (d=4) 6000/1500/2500
1.0 log10 vol.
1.5 0.7
0.8
0.9 Joint cov.
1.0
3.8
Households (d=4) 4324/1081/1802
CHR
CHR
CHR
TSCP
SC2 rect
SC2 rect
Max 10
0.7
0.8
0.9 Joint cov.
18.5
Air (d=6) 6000/1500/2500
19.0 19.5 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
4.2 4.4 log10 vol.
ATP1d (d=6) 202/50/85
SC2 =1
SC2 =0.8
SC2 =1
SC2 =1
Max
Max
SC2 rect
TSCP
TSCP
SC2 =0.8
SC2 rect
SC2 rect
Max
CHR 8 log10 vol.
10 0.7
0.8
0.9 Joint cov.
1.0
21.0 21.5 0.7 log10 vol.
0.8
0.9 Joint cov.
1.0
18.5
RF2 (d=8) 4607/1151/1921
CHR
CHR
SC2 =1
SC2 rect
SC2 =0.8
SC2 rect
SC2 =0.8
CHR
TSCP
TSCP
SC2 rect
Max 13 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
1.0
19.0 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
OSales (d=12) 333/83/140
SC2 =0.8
12
0.9 Joint cov.
TSCP
SC2 =1
Max
0.8
CHR 20.5
RF1 (d=8) 5403/1350/2252 SC2 =1
0.7
ATP7d (d=6) 177/44/75
SC2 =0.8
CHR
6
1.0
Max
1.0
TSCP
0.9 Joint cov.
SC2 =0.8
SC2 =0.8
log10 vol.
0.8
TSCP
TSCP
5
0.7
SC2 =1
SC2 =1
Max
3.9 4.0 log10 vol.
Stock (d=4) 570/142/238
SC2 =0.8
SC2 rect
1.0
CHR
1.0
SC2 =1
0.9 Joint cov.
SC2 =0.8
Max 0.7
0.8
SC2 =1
TSCP
2.0 log10 vol.
6 0.7
Student (d=3) 237/59/99
SC2 =0.8
1.5
5 log10 vol.
TSCP
SC2 rect
CHR
1.0
Max 3.5
SC2 =1
SC2 =1
0.9 Joint cov.
SC2 =1
SF1 (d=3) 193/48/82 Max
0.8
SCPF (d=3) 682/170/285
TSCP
0.7
1.44 0.7
CHR
SC2 =1
0.650 0.675 0.700 log10 vol.
1.45 log10 vol.
Max 10
11 log10 vol.
29
120.7
0.8
0.9 Joint cov.
1.0
62.5
65.0 67.5 log10 vol.
70.00.7
0.8
0.9 Joint cov.
1.0
WQ (d=14) 636/159/265
OES10 (d=16) 241/60/102
OES97 (d=16) 200/50/84
Max
SC2 =0.8
SC2 =0.8
SC2 rect
SC2 =1
SC2 =1
TSCP
TSCP
TSCP
SC2 rect
SC2 rect
CHR
Max
CHR
SC2 =0.8
CHR
SC2 =1
13.5
14.0 14.5 log10 vol.
0.7
0.8
0.9 Joint cov.
1.0
Max 62
SCM1d (d=16) 5881/1470/2452
64 log10 vol.
66 0.7
0.8
0.9 Joint cov.
1.0
SCM20d (d=16) 5379/1344/2243
SC2 =1
SC2 =0.8
SC2 =0.8
SC2 =1
TSCP
TSCP
CHR
CHR
SC2 rect
SC2 rect
Max
Max 46.6
46.8 log10 vol.
47.0 0.7
0.8
0.9 Joint cov.
1.0
48.8
49.0 49.2 0.7 log10 vol.
30
0.8
0.9 Joint cov.
1.0
64
66 log10 vol.
68 0.7
0.8
0.9 Joint cov.
1.0