Probabilistic Conformal Prediction with Approximate Conditional Validity

## Abstract

Abstract We develop a new method for generating prediction sets that combines the flexibility of conformal methods with an estimate of the conditional distribution P Y ∣ X subscript P conditional 𝑌 𝑋 \textup{P}_{Y\mid X} P start_POSTSUBSCRIPT italic_Y ∣ italic_X end_POSTSUBSCRIPT . Existing methods, such as conformalized quantile regression and probabilistic conformal prediction, usually provide only a marginal coverage guarantee. In contrast, our approach extends these frameworks to achieve approximately conditional coverage, which is crucial for many practical applications. Our prediction sets adapt to the behavior of the predictive distribution, making them effective even under high heteroscedasticity.
While exact conditional guarantees are infeasible without assumptions on the underlying data distribution, we derive non-asymptotic bounds that depend on the total variation distance of the conditional distribution and its estimate. Using extensive simulations, we show that our method consistently outperforms existing approaches in terms of conditional coverage, leading to more reliable statistical inference in a variety of applications.

## 1 Introduction

Conformal predictions are commonly used to construct prediction sets. Under minimal assumptions, they offer finite-sample validity (Vovk et al., ; Shafer & Vovk, ). However, significant challenges arise with high heteroskedasticity, often leading to incorrect inferences (Dewolf et al., ). The split-conformal approach uses a set of $n$ calibration data points $\{(X_{k},Y_{k})\}_{k\in[n]}$ with $X_{k}\in\mathbb{R}^{d}$ and $Y_{k}\in\mathcal{Y}$ to create a prediction set $\mathcal{C}_{\alpha}(x)$ where $\alpha\in(0,1)$. For each $x\in\mathbb{R}^{d}$, the prediction set based on a conformity score function $V\colon\mathbb{R}^{d}\times\mathcal{Y}\to\mathbb{R}$, is given by

$$ $\mathcal{C}_{\alpha}(x)=\left\{{y\in\mathcal{Y}\colon V(x,y)\leq Q_{1-\alpha} \left({\frac{1}{n+1}\sum\nolimits_{k=1}^{n}\delta_{V(X_{k},Y_{k})}+\frac{1}{n+ 1}\delta_{\infty}}\right)}\right\},$ (1) $$

where $\delta_{x}$ is the Dirac mass and $Q_{1-\alpha}$ is the $(1-\alpha)$-quantile of the adjusted empirical score distribution $\frac{1}{n+1}\sum_{k=1}^{n}\delta_{V(X_{k},Y_{k})}+\frac{1}{n+1}\delta_{\infty}$.
If the calibration data $\{(X_{k},Y_{k})\}_{k\in[n]}$ is drawn i.i.d. from a population distribution $\textup{P}_{X,Y}$, then for any new data point $(X_{n+1},Y_{n+1})\sim\textup{P}_{X,Y}$ sampled independently of the calibration data, the conformal theory ensures the marginal validity of $\mathcal{C}_{\alpha}(X_{n+1})$, meaning that
$\mathbb{P}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})}\right)\geq 1-\alpha$.
This marginal guarantee can hide significant discrepancies in the coverage of different regions of the input space $\mathbb{R}^{d}$; see e.g. Izbicki et al. (). Conditional validity is a more desirable guarantee than marginal validity:
for any $x\in\mathbb{R}^{d}$, the set $\mathcal{C}_{\alpha}(x)$ is conditionally valid if

$$ $\mathbb{P}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})\mid X_{n+1}=x}\right) \geq 1-\alpha.$ (2) $$

However, this property cannot be achieved without further assumptions about the data distribution; see Vovk (); Lei & Wasserman ().
For practical purposes, it is enough to construct sets $\mathcal{C}_{\alpha}$ which approximate (2), and ideally achieve it asymptotically under suitable conditions in the limit of large sample size $n$.

Related work. A great deal of research has been devoted to this problem, starting with the case where $\mathcal{Y}=\mathbb{R}$.
For example Romano et al. (); Kivaranovic et al. () proposed methods based on estimate of the lower and upper conditional quantile function $\hat{q}_{\alpha/2}$ and $\hat{q}_{1-\alpha/2}$, to define a quantile-based conformity score:
$V(x,y)=\max\left\{{\hat{q}_{\alpha/2}(x)-y,y-\hat{q}_{1-\alpha/2}(x)}\right\}$ which is then conformalized.
Improvements of this conformity score are investigated in Kivaranovic et al. (); Sesia & Candès (). Sesia & Candès () have shown that the constructed interval converges to the narrowest possible interval that achieve conditional coverage under mild assumptions. We stress that these methods are specific to the case $\mathcal{Y}=\mathbb{R}$; in addition, when the conditional distribution $\textup{P}_{Y\mid X}$ is multimodal, restricting prediction to intervals is suboptimal; see Wang et al. () for a discussion and examples.
It has been suggested in Guan (); Alaa et al. () to partition $\mathbb{R}^{d}$ and learn a specific quantile for each regions. Splitting the space $\mathbb{R}^{d}$ into multiple regions typically leads to an increase in the length of the prediction set; see Romano et al. (); Melki et al. () for comments.

A number of studies have focused on the construction of prediction sets using an estimator of the conditional distribution $\textup{P}_{Y\mid X}$. Again, in the case where $\mathcal{Y}=\mathbb{R}$, Cai et al. (); Lei & Wasserman () have constructed prediction intervals based on an estimator of the conditional density and have established the asymptotic validity under appropriate conditions. Han et al. () use kernel density estimation to construct asymmetric prediction bands. However, this method cannot handle bimodality as it generates a single interval.
On the other hand, Sesia & Romano () partition the domain of $Y$ into bins to create a histogram approximation of $\textup{P}_{Y\mid X}$. The authors showed that their method satisfies the marginal validity while achieving the asymptotic conditional coverage; see also Lei et al. (). Asymptotic conditional coverage is also obtained in Sesia & Candès (); Cauchois et al. () using quantile regression-based methods, or using cumulative distribution function estimators (Izbicki et al., ; Chernozhukov et al., ). Conditionally valid prediction sets have been shown to improve the robustness to perturbations (Gendler et al., ).
Guha et al. () proposed a novel approach that converts regression tasks into classification problems by binning the output space and discretizing the labels. Leveraging this discretization, they approximate the conditional density to construct prediction sets that correspond to regions of Highest Predictive Density (HPD). A key limitation of these methods lies in the discretization process, as the number of labels required for complex scenarios can become computationally prohibitive. Diamant et al. () introduced an approach that estimates conditional densities using neural networks parameterized by splines, offering a more flexible representation.

Very few studies have addressed the scenario where the prediction target is multi-dimensional, i.e., $\mathcal{Y}=\mathbb{R}^{q}$ with $q>1$.
Wang et al. () developed the PCP method, based on implicit conditional generative models (CGMs). These CGMs allow for the generation of samples from the conditional distribution without requiring an explicit closed-form expression. The PCP method constructs prediction sets as unions of balls, whose centers are generated from the CGM. However, in PCP, the radius of these balls is fixed across the space, which can introduce significant limitations. In regions with low variability, a fixed radius may result in over-coverage, while for highly dispersed conditional distributions, it may lead to under-coverage, failing to capture the full extent of the relevant space.
We observed that the performance of PCP deteriorates as the number of balls increases, exacerbating the heteroscedasticity problem. This underscores the need for a more adaptive methodology that can dynamically adjust the size of prediction sets in response to local variability in the data.

Our work addresses these challenges through the following main contributions:

- •
We propose a new CP^2 method for constructing conditional confidence sets that adapts to the local structure of the data distribution, capable of addressing both classical regression problems where $\mathcal{Y}=\mathbb{R}$ and more complex multi-dimensional prediction tasks where $\mathcal{Y}=\mathbb{R}^{q}$. Our approach is versatile in accommodating scenarios involving either an explicit conditional density estimator or an implicit generative model; see Section 2.
- •
We develop a theoretical framework to analyze the properties of the proposed CP^2 method, establishing both its marginal and approximate conditional validity. Furthermore, we demonstrate that asymptotic conditional coverage is attainable under a weak consistency assumption on the predictive distribution; see Section 3.
- •
We demonstrate the effectiveness of the proposed method through a series of experiments on synthetic and real-world datasets. The results indicate that our approach consistently outperforms existing methods in terms of conditional coverage. Specifically, it excels in handling classical regression problems, effectively addressing multimodality, and proves robust in the more challenging setting of multidimensional prediction tasks; see Section 4.

## 2 The CP 2 framework

We want to construct marginally valid predictive sets with approximate conditional validity. We follow the split-conformal approach to conformal inference due to its computational feasibility with large datasets; see along others Papadopoulos et al. (); Papadopoulos (); Papadopoulos et al. (); Romano et al. (); Kivaranovic et al. (). The first step involves splitting the data samples into two disjoint subsets, the training set $\mathcal{T}=\{(\tilde{X}_{k},\tilde{Y}_{k})\}_{k=1}^{m}$ and calibration set $\mathcal{C}=\{(X_{k},Y_{k})\}_{k=1}^{n}$. It is assumed in the sequel that the training and calibration data are mutually independent and i.i.d. with distribution $\textup{P}_{X,Y}$ over the feature vectors $X\in\mathbb{R}^{d}$ and response variables $Y\in\mathcal{Y}$.
The target set $\mathcal{Y}$ can be either finite or continuous. We must now predict the unknown value of $Y_{n+1}$ given an input $X_{n+1}$, independent from $\mathcal{T}\cup\mathcal{C}$. Our goal is to construct a prediction set $\mathcal{C}_{\alpha}(X_{n+1})$ that contains the unobserved output $Y_{n+1}$ with probability close to $1-\alpha$, where $\alpha\in(0,1)$ is the user-specified confidence level.

An estimator $\Pi_{Y\mid X}$ of the conditional probability $\textup{P}_{Y\mid X}$ is learnt using the training data. There is a rich body of research on nonparametric conditional density estimation, with the most common methods relying on smoothing techniques such as kernel smoothing and local polynomial fitting.
An alternative approach involves transforming the conditional density estimation task into a regression problem, allowing the application of nonparametric regression methods to approximate the conditional density. More recently, generative methods leveraging deep neural networks have been developed for nonparametric conditional density estimation, which enable sampling from the conditional distribution; see Abadi et al. (); Zhou et al. () for examples of these techniques.
In the sequel, the choice of this algorithm is treated as a black box. There are three main ingredients for our approach:

- 1.
In classical conformal prediction methods, the shape of the prediction set is specified by the score function $V(x,y)$; see (1). The first element of our construction is a family of explicitly given confidence sets $\mathcal{R}_{z}(x;t)$ parameterized by $t\in\mathsf{T}$ where $\mathsf{T}$ is a subset of $\mathbb{R}$ and $z\in\mathcal{Z}$ auxiliary variables. The index set $\mathsf{T}$ can, in most cases, be taken as either $\mathsf{T}=\mathbb{R}$ or $\mathsf{T}=\mathbb{R}_{+}$. The key assumptions for the confidence set are as follows:
(a) The size of $\mathcal{R}_{z}(x;t)$ increases with $t\in\mathsf{T}$ for any $z\in\mathcal{Z}$. In addition, by choosing a sufficiently large value of $t$, the entire output space $\mathcal{Y}$ can be covered.
(b) There exists a form of continuity for $t\mapsto\mathcal{R}_{z}(x;t)$. In mathematical terms, the following assumption should hold:
H 1.
For any $(x,z)\in\mathbb{R}^{d}\times\mathcal{Z}$, the confidence sets $\{\mathcal{R}_{z}(x;t)\}_{t\in\mathsf{T}}$ are non-decreasing, $\Pi_{Y|X=x}(\cap_{t\in\mathsf{T}}\mathcal{R}_{z}(x;t))=0$, $\cup_{t\in\mathcal{T}}\mathcal{R}_{z}(x;t)=\mathcal{Y}$; in addition, for any $t\in\mathsf{T}$, $\cap_{t^{\prime}>t}\mathcal{R}_{z}(x;t^{\prime})=\mathcal{R}_{z}(x;t)$.
Example 2.1.
For instance, $\mathcal{R}(x;t)$ can be chosen as a ball centered around an estimate of conditional mean $\textup{P}_{Y\mid X}$ with radius $t$ and there is no auxiliary variables then (we then remove subscript $z$). Examples of confidence intervals specialized to the case where $\mathcal{Y}=\mathbb{R}$ are given in Table 1.
Example 2.2.
If the predictive distribution is multimodal, a ball centered around the predictive mean often fails to provide an informative prediction set. Ideally, $\mathcal{R}_{z}(x;t)$ should correspond to the set with the highest predictive density (HPD) of $\textup{P}_{Y\mid X}$. However, HPD regions are difficult to determine in practice, even when the conditional predictive density is available. Following Wang et al. (), we may define the prediction sets as $\mathcal{R}_{z}(x;t):=\cup_{i=1}^{M}\mathrm{B}(y_{i},t)$, and they depends on an exogenous variables $z=(y_{1},\ldots,y_{M})\in\mathcal{Z}$, where each $y_{i}$ is sampled conditionally independently from the conditional generative model $\Pi_{Y\mid X=x}$.
In a general setting, the auxiliary variables $z\in\mathcal{Z}$ are sampled according to a kernel $\bar{\Pi}_{Z\mid X=x}$ for a given $x$.
- 2.
To localize conformal prediction methods, it is convenient to introduce a function $f_{\tau}(\lambda)$ parameterized by $\tau$. Such function aims to transform the conformity score $\lambda$ and was introduced (albeit in a slightly different form) in Deutschmann et al. (); Han et al. (). Examples of such a function are $f_{\tau}(\lambda)=\tau\lambda$ and $f_{\tau}(\lambda)=\tau+\lambda$. We assume that:
H 2.
There exists $\varphi\in\mathsf{T}$ such that $\tau\in\mathsf{T}\mapsto f_{\tau}(\varphi)$ is increasing and bijective.
In addition, $\lambda\in\mathsf{T}\mapsto f_{\tau}(\lambda)$ is increasing for any $\tau\in\mathsf{T}$.
We define $\tau_{x,z}$ using the estimated predictive density $\Pi_{Y\mid X=x}$ according to
$\tau_{x,z}=\inf\left\{{\tau\in\mathsf{T}\colon\Pi_{Y\mid X=x}(\mathcal{R}_{z}(
x;f_{\tau}(\varphi)))\geq 1-\alpha}\right\}.$
(3)
It is easily shown that for any $\alpha\in(0,1)$, $x\in\mathbb{R}^{d}$, $z\in\mathcal{Z}$, then $\tau_{x,z}\in\mathsf{T}$ and $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))\geq 1-\alpha$; see Lemma A.3.
- 3.
Finally, we need to introduce the conformity score for the considered setup. The natural choice is the minimal size of the set required to cover the observation $y$ at the input $x$ for the auxiliary variables $z$: $\lambda_{x,y,z}=\inf\left\{{t\in\mathsf{T}\colon y\in\mathcal{R}_{z}(x;t)}\right\}$.
- 4.
The resulting procedure works as follows. For $k\in\{1,\dots,n\}$, we set $\bar{\tau}_{k}:=\tau_{X_{k},Z_{k}}$ and $\bar{\lambda}_{k}:=\lambda_{X_{k},Y_{k},Z_{k}}$ where $\{Z_{k}\}_{k=1}^{n}$ are sampled conditionally independently from $\bar{\Pi}_{Z\mid X=X_{k}}$. Given $X_{n+1}\in\mathbb{R}^{d}$, we sample $Z_{n+1}\sim\bar{\Pi}_{Z\mid X=X_{n+1}}$ conditionally independently from $\{(X_{k},Y_{k},Z_{k})\}_{k=1}^{n}$,
and construct the resulting CP^2 prediction set as
$\textstyle\mathcal{C}_{\alpha}(X_{n+1})=\mathcal{R}_{Z_{n+1}}\left({X_{n+1};f_
{\bar{\tau}_{n+1}}\bigl{(}Q_{1-\alpha}(\mu_{n})\bigr{)}}\right),$
(4)
where $Q_{1-\alpha}(\mu_{n})$ is the $1-\alpha$ quantile of the distribution $\mu_{n}$, and is given by
$\mu_{n}=\frac{1}{n+1}\sum\nolimits_{k=1}^{n}\delta_{f_{\bar{\tau}_{k}}^{-1}(
\bar{\lambda}_{k})}+\frac{1}{n+1}\delta_{\infty}.$
(5)
The transformation $\{v\mapsto f_{\tau}(v)\}_{\tau\in\mathsf{T}}$ balances the following two factors: (a) The optimal parameter $\lambda_{x,y,z}$ ensuring that $y$ is included in the confidence set $\mathcal{R}_{z}(x;\lambda_{x,y,z})$; (b) The parameter $\tau_{x,z}$ obtained from the probabilistic model $\Pi_{Y\mid X=x}$.

Table: Table 1: Confidence sets $\mathcal{R}(x;t)$ found in the literature and also discussed in Gupta et al. ().

We stress that CP^2 is a general framework that can be adapted to many choices for conditional predictive density estimates, constructing the family of confidence sets, and selecting the calibration function $f_{\tau}(\lambda)$. However, we start with the simple example that shows that CP^2 is more general than the classical split-conformal CP approach.

Figure: Figure 1: Predictions sets obtained via the standard CP and CP^2 methods.
Refer to caption: x1.png

###### H 1 .

###### Example 2.1 .

###### Example 2.2 .

###### H 2 .

##### Simple example of CP 2 .

We begin with a simple application of CP^2 to highlight its differences from the basic conformal approach, with $\mathcal{Y}=\mathbb{R}$. The calibration sets are defined as: $\mathcal{R}(x;t)=\{y\in\mathcal{Y}\colon|y-\text{pred}(x)|\leq t\}$ for $t\in\mathbb{R}_{+}$. There are no auxiliary variables $z$ in this case, so we omit $z$ from the notation. Assumption 1 is easily satisfied with $\mathsf{T}=\mathbb{R}_{+}$. We then take $f_{\tau}(\lambda)=\tau\lambda$, $\tau\in\mathbb{R}_{+}$ and $\varphi=1$: 2 is also satisfied. Note that $f^{-1}_{\tau}(\lambda)=\lambda/\tau$ for $\tau\in\mathbb{R}_{+}^{*}$.
Using the CP^2 approach, we find that $\lambda_{x,y}=|y-\mathrm{pred}(x)|$, which corresponds to a standard conformity score. The classical conformal prediction method defines the $1-\alpha$ quantile based on the associated empirical measure $\nu_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{\bar{\lambda}_{k}}+\frac{1}{n+1}
\delta_{\infty}$, where $\bar{\lambda}_{k}=|Y_{k}-\mathrm{pred}(X_{k})|$. CP^2 differs from the basic conformal approach by introducing $\tau_{x}=\operatornamewithlimits{\arg\min}\{\tau\in\mathsf{T}\colon\Pi_{Y\mid X
=x}([\mathrm{pred}(x)\pm\tau])\geq 1-\alpha\}$ as in (3), where $[\mathrm{pred}(x)\pm\tau]=[\mathrm{pred}(x)-\tau,\mathrm{pred}(x)+\tau]$.
The prediction set becomes $[\mathrm{pred}(x)\pm f_{\tau_{x}}(Q_{1-\alpha}(\mu_{n}))]$, where $\mu_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{\bar{\lambda}_{k}/\tau_{k}}+\frac{1
}{n+1}\delta_{\infty}$.
In this setting, the choice of $\varphi$ is irrelevant. Take $\varphi>0$ and denote $\tau^{\varphi}_{x}=\operatornamewithlimits{\arg\min}\{\tau\in\mathsf{T}\colon
\Pi_{Y\mid X=x}([\mathrm{pred}(x)\pm f_{\tau}(\varphi)])\geq 1-\alpha\}$. It is easily seen that $\tau_{x}^{\varphi}=\tau_{x}/\varphi$. The prediction set becomes
$[\mathrm{pred}(x)\pm f_{\tau^{\varphi}_{x}}(Q_{1-\alpha}(\mu^{\varphi}_{n}))]$, where $\mu^{\varphi}_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{\bar{\lambda}_{k}/\bar{
\tau}_{k}^{\varphi}}+\frac{1}{n+1}\delta_{\infty}$ with $\bar{\tau}_{k}^{\varphi}=\tau^{\varphi}_{X_{k}}$. Note that $Q_{1-\alpha}(\mu^{\varphi}_{n}))=\varphi Q_{1-\alpha}(\mu_{n}))$ and thus
$f_{\tau^{\varphi}_{x}}(Q_{1-\alpha}(\mu^{\varphi}_{n}))=f_{\tau_{x}}(Q_{1-
\alpha}(\mu_{n}))$ showing that, the prediction set does not depend on the choice of $\varphi$. To illustrate the advantage of our method, in Figure 1 we present the prediction sets obtained with the classical CP method and CP^2 in the case of a Neal’s funnel-shaped distribution in 2 dimensions; see Neal (, Section 9).

The natural approach for the general case is to use the conditional distribution $\Pi_{Y\mid X=x}$ or its estimate to find Highest Predictive Density (HPD) regions and calibrate their size with the help of CP^2. We develop the respective general algorithm CP^2-HPD in Appendix B.1. However, the procedures to find HPDs are usually highly non-trivial and we only investigate this approach experimentally for synthetic data; see Section 4.1. Next, we provide a specific implementation of our general CP^2 framework that is universally applicable.

##### CP 2 with Implicit Conditional Generative Model: CP 2 -PCP .

We also develop a second instance of the CP^2 algorithm, inspired by Wang et al. ().
Unlike CP^2-HPD, this approach does not require the conditional density. Instead, it is designed for cases where the conditional generative model (CGM) $\Pi_{Y\mid X}$ is implicit, meaning we cannot evaluate it pointwise while being able to sample from it.
For each calibration point $X_{k}$, we draw $M$ random variables $\{\hat{Y}_{k,i}\}_{i=1}^{M}$ from $\Pi_{Y\mid X=X_{k}}$.
We denote $Z_{k}=(\hat{Y}_{k,1},\ldots,\hat{Y}_{k,M})$ and consider the confidence sets as the union of spheres centered around the sample points $\mathcal{R}_{Z_{k}}(X_{k};t)=\cup_{i=1}^{M}\mathrm{B}(\hat{Y}_{k,i},t)$. With such choice, we get $\bar{\lambda}_{k}=\min_{i=1}^{M}\|Y_{k}-\hat{Y}_{k,i}\|$.
We then draw a second sample $\{\tilde{Y}_{k,j}\}_{j=1}^{\tilde{M}}$, and compute $\bar{\tau}_{k}=\{t\in\mathbb{R}_{+}\colon\tilde{M}^{-1}\sum_{j=1}^{\tilde{M}}
\mathds{1}_{\mathcal{R}_{Z_{k}}(X_{k};f_{t}(\varphi))}(\tilde{Y}_{k,j})\geq 1-\alpha\}$.
It is easily seen that

$$ $\textstyle\bar{\tau}_{k}=\left({t\mapsto f_{t}(\varphi)}\right)^{-1}\left\{{Q_ {1-\alpha}\left({\frac{1}{\tilde{M}}\sum_{j=1}^{\tilde{M}}\delta_{\min_{i=1}^{ M}\|\tilde{Y}_{k,j}-\hat{Y}_{k,i}\|}}\right)}\right\}.$ $$

Given a new input $X_{n+1}\in\mathbb{R}^{d}$, we sample $Z_{n+1}=(\hat{Y}_{n+1,1},\ldots,\hat{Y}_{n+1,M})$ and obtain prediction set as follows

$$ $\textstyle\mathcal{C}_{\alpha}(X_{n+1})=\left\{{y\in\mathcal{Y}\colon\min_{i=1 }^{M}\|y-\hat{Y}_{n+1,i}\|\leq f_{\bar{\tau}_{n+1}}\left({Q_{1-\alpha}\left({ \mu_{n}}\right)}\right)}\right\},$ $$

where $\mu_{n}$ is given in (5).
The CP^2-PCP method employs the same confidence set $\mathcal{R}_{z}(x;t)$ as the one used by PCP. This method effectively captures multimodalities using balls centered at likely outputs $\hat{Y}_{n+1,i}$.
Furthermore, the conformity scores used by PCP correspond to our $\lambda_{x,y,z}$. However, the key distinction between the two algorithms lies in the additional parameter $\tau_{x,z}$ for CP^2-PCP, which requires the generation of a second random sample from $\Pi_{Y\mid X=x}$. This method is especially useful when solving equation (3) is intractable. We summarize CP^2-PCP in Algorithm 1.

Figure: Algorithm 1 CP^2-PCP

## 3 Theoretical Guarantees

In this section, we provide both marginal and conditional guarantees for the prediction set $\mathcal{C}_{\alpha}(x)$ given in (4).
The validity of these guarantees is ensured by the exchangeability of the calibration data, with the exception of Theorem 3.3 which relies on a concentration inequality and thus requires i.i.d. calibration data.
The following theorem establishes marginal validity of the predictive set defined by CP^2.

###### Theorem 3.1 .

Assume 1-2.
Then, for any $\alpha\in(0,1)$, it holds $1-\alpha\leq\mathbb{P}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})}\right)$.
Moreover, if the conformity scores $\{f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})\}_{k=1}^{n+1}$ are almost surely distinct, then it also holds that
$\mathbb{P}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})}\right)<1-\alpha+(n+1
)^{-1}.$

The proof is postponed to Section A.1. Moreover, the upper bound on the coverage always holds when the distribution of $f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})$ is continuous.
Now, we will investigate the conditional validity.
Denote by $\mathrm{d}_{\mathrm{TV}}$ the total variation distance and by $\mathbb{P}^{\mathcal{T}}$ the conditional probability given the training data.

###### Theorem 3.2 .

Assume 1-2, and let $\alpha\in(0,1)$.
For any $x\in\mathbb{R}^{d}$ and $z\in\mathcal{Z}$, it holds

$$ $\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(x)\mid(X_{n+1},Z _{n+1})=(x,z)}\right)\geq 1-\alpha-\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X =x};\Pi_{Y\mid X=x})-p_{n+1}^{(x,z)},$ $$

where
$p_{n+1}^{(x,z)}=\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\mu_{n})<f_{\bar{
\tau}_{n+1}}^{-1}(\bar{\lambda}_{n+1})\leq\varphi\,|\,(X_{n+1},Z_{n+1})=(x,z)}\right)$.

The proof is postponed to Section A.1.
The more accurately the estimator $\Pi_{Y\mid X=x}$ approximates the true conditional distribution, the closer the result will be to $1-\alpha$. The second term in the lower bound is $p_{n+1}^{(x,z)}$. Its expected value is upper bounded by $\mathbb{E}[p_{n+1}^{(X,Z)}]\leq\alpha$, and non-asymptotic bounds for this error term are developed in Section A.2.

We will now briefly discuss the asymptotic conditional coverage guarantee; details are provided in the supplementary paper. Assuming the availability of an oracle for the predictive distribution, i.e., $\textup{P}_{Y\mid X=x}=\Pi_{Y\mid X=x}$,we get under 1 and 2, that for any $t\in\mathbb{R}$,

$$ $\displaystyle\mathbb{P}\left({\lambda_{X,Y,Z}\leq f_{\tau_{X,Z}}(t)\,|\,X=x,Z= z}\right)$ $\displaystyle=\mathbb{P}\left({Y\in\mathcal{R}_{z}(x;f_{\tau_{x,z}}(t))\,|\,X= x,Z=z}\right)$ $\displaystyle=\Pi_{Y\mid X=x}\left({\mathcal{R}_{z}(x;f_{\tau_{x,z}}(t))} \right),$ $$

where $(X,Y,Z)$ follows the same distribution than $(X_{k},Y_{k},Z_{k})$, $k\in\{1,\dots,n\}$.
Note that $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(t)))\geq 1-\alpha$ if and only if $t\geq\varphi$, which implies that

$$ $\mathbb{P}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})\leq t\mid(X,Z)=(x,z))\geq 1-\alpha$ if and only if $t\geq\varphi$ . (6) $$

From (6) it is easily seen that the $(1-\alpha)$-quantile of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ is $\varphi$.
The Glivenko–Cantelli Theorem (Van der Vaart, , Theorem 19.1) demonstrates that $\sup_{t\in\mathbb{R}}|\mu_{n}(-\infty,t]-\mathbb{P}(f_{\tau_{X,Z}}^{-1}(
\lambda_{X,Y,Z})\leq t)|\to 0$ almost surely as $n\to\infty$, where $\mu_{n}$ is defined in (5).
Since the convergence of the c.d.f. implies the convergence of the quantile function (Van der Vaart, , Lemma 21.2), we deduce that $Q_{1-\alpha}(\mu_{n})\to\varphi$ almost-surely as $n\to\infty$.
Under weak additional conditions this implies that $\lim_{n\to\infty}p_{n+1}^{(x,z)}=0$, $\bar{\Pi}_{Z\mid X}\times\textup{P}_{X}$-almost everywhere, where $\bar{\Pi}_{Z\mid X}$ is the Markov kernel used to draw the auxiliary variables $z$; see Section A.4. In this case, Theorem 3.1 implies the asymptotic validity of CP^2.

In practice, the oracle is unavailable. In the following theorem, we examine the asymptotic conditional conformal validity as the size of the training dataset,
$m_{n}$, goes to infinity with $n$. In most cases, $\lim_{n\to\infty}m_{n}/n=\gamma>0$, but this is not required here.
To make the dependency of the estimator on the size of the training set explicit, we will denote the conditional distribution a $\Pi_{Y\mid X}^{(m_{n})}$.
Consider the following assumption.

###### H 3 .

There exists sequence $(r_{n})$ such that
$\lim\limits_{n\to\infty}\mathbb{P}(\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y};
\textup{P}_{X}\times\Pi_{Y\mid X}^{(m_{n})})\leq r_{n})=1.$

In most interesting case, we have $\lim_{n\to\infty}r_{n}=0$. Such types of bounds can be deduced from Devroye & Lugosi (, Chapter 9).
Let $(X,Y,Z)$ and $(X,\hat{Y},Z)$ be random variables distributed according to $\textup{P}_{X,Y}\times\bar{\Pi}_{Z\mid X}$ and $\textup{P}_{X}\times\Pi_{Y\mid X}^{(m_{n})}\times\bar{\Pi}_{Z\mid X}$, respectively.

###### Theorem 3.3 .

Assume 1-2-3 hold.
If the distributions of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ and $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$ are continuous, then, it holds

$$ $\left\lvert{\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n +1})\,|\,X_{n+1},Z_{n+1}}\right)-1+\alpha}\right\rvert=\operatorname{\mathrm{O }}_{\mathbb{P}}\left({\sqrt{n^{-1}\log n}+r_{n}}\right).$ $$

In Lei et al. (); Izbicki et al. (); Sesia & Candès (), the asymptotic conditional validity is demonstrated by assuming the consistency of their methods’ estimators. For instance, Romano et al. () assume that the conditional quantile regressor converges in $L^{2}$ towards the true quantile with high probability.

###### Theorem 3.1 .

###### Theorem 3.2 .

###### H 3 .

###### Theorem 3.3 .

## 4 Numerical Experiments

Figure: (a) Data and MDN estimate
Refer to caption: x2.png

In this section, we conduct a comprehensive analysis demonstrating the advantage of CP^2 compared to standard and adaptive split conformal algorithms. Specifically, we benchmark our algorithm against several state-of-the-art methods: Conformalized Quantile Regression (Romano et al., ), Conformalized Histogram Regression (Sesia & Romano, ) and Probabilistic Conformal Prediction (Wang et al., ). All these method share some key aspects: they are built on top of the pre-trained models and do not require access to training data or the model’s internals on both calibration and prediction steps. We aim to answer these specific questions: how does CP^2 performs in terms of coverage, conditional coverage and predictive set volume when compared to state-of-the-art methods on synthetic and real data.

### 4.1 Synthetic Data Experiment

In this example, $(X_{k},Y_{k})$ is sampled from a mixture of $P=4$ Gaussians; see Figure 2(a). The number of training and calibration samples is $m=10^{4}$ and $n=10^{3}$, respectively. We fit a Mixture Density Network (MDN) as an explicit generative model, $\gamma_{Y\mid X=x}(y)=\sum_{\ell=1}^{P}\pi_{\ell}(x)\mathcal{N}(y;\mu_{\ell}(x
),\sigma^{2}_{\ell}(x))$, where $\mu_{\ell}(\cdot)$, $\sigma_{\ell}(\cdot)$ and $\pi_{\ell}(\cdot)$ are all modeled by fully connected 2-layers neural networks (the condition $\sum_{\ell=1}^{P}\pi_{\ell}(x)=1$ is ensured by using softmax activation functions).
We use CP^2-HPD (the calculation of the HPD rates as well as $\tau_{x}$ and $\lambda_{x,y}$ is explicit in this case) with $f_{t}(v)=tv$.
The parameters of the MDN are trained by maximizing the likelihood on the training set.

We compare the plain CP^2-HPD, PCP (with the same MDN as CP^2-HPD and $M=50$ draws) and CQR. All methods achieve the desired marginal coverage $1-\alpha=0.9$. We illustrate the conditional coverage in Figure 2(b) and the lengths of the predictive sets in Figure 2(c).
CP with a fixed-width predictive set performs poorly in this multimodal example, both in terms of the size of the confidence set and the conditional coverage.
CP^2-HPD and CQR perform similarly in terms of conditional coverage (which remains close to $1-\alpha=0.9$).
The conditional coverage of PCP varies between 0.85 and 0.95.
CP^2-HPD produces shorter prediction sets compared to CQR and PCP. This is because CP^2-HPD uses an HPD confidence set that is more suitable for multimodal applications than the interval produced by CQR.

Figure: Figure 3: Worst-slab coverage on real data. Results averaged over 50 random splits of each dataset. Calibration and test set sizes set to 2000, 50 conditional samples for PCP, CP^2 and $\Pi_{Y\mid X}$. Worst-slab coverage parameter $(1-\delta)=0.1$. Nominal coverage level is $(1-\alpha)=0.9$ and is shown in dashed black. Methods with conditional coverage below $0.75$ shown as cross-hatched on horizontal axis.
Refer to caption: x5.png

Figure: Figure 4: Sizes of the prediction sets on real data. We divide the size of the set by the standard deviation of response to present the results on the same scale.
Refer to caption: x6.png

### 4.2 Real-world Regression Data Experiments

In this section, we study the performance of CP^2-PCP on several real world regression datasets. In particular, we used the same datasets as in Wang et al. () to facilitate comparisons.

Datasets.
We use publicly available regression datasets, which are also considered in Romano et al. (); Wang et al. (). Some of them come from the UCI repository: bike sharing (bike), protein structure (bio), blog feedback (blog), Facebook comments (fb1 and fb2). Other datasets come from US Department of Health surveys (meps19, meps20 and meps21), and from weather forecasts (temp) (Cho et al., ).

Methods.
We compare the proposed CP^2-PCP method with Probabilistic Conformal Prediction (PCP; Wang et al. ()), Conformalized Quantile Regression (CQR; Romano et al. ()) and Conformalized Histogram Regression (CHR; Sesia & Romano ()). We also consider CQR$2$ which is a modification of CQR that uses inverse quantile nonconformity score. For our method and PCP we use a Mixture Density Network Bishop () to estimate the conditional distribution $\textup{P}_{Y\mid X}$, since it was chosen in Wang et al. () as best-performing. We also consider different choices of $f_{t}$ for our method: CP^2-PCP-L stands for CP^2-PCP with $f_{t}(v)=tv$ and CP^2-PCP-D stands for CP^2-PCP with $f_{t}(v)=t+v$. Our implementation of CP^2-PCP is summarized in Algorithm 1.
Additionally, we consider $\Pi_{Y\mid X}$ which is a special case of CP^2-PCP with $f_{t}(v)=t$.

Metrics.
Empirical coverage (marginal and conditional) is the main quantity of interest for prediction sets. We evaluate worst-slab conditional coverage (Cauchois et al., ; Romano et al., ) in our experiments, see details in Appendix B.3. We also measure the total size of the predicted sets, scaled by the standard deviation of the response $Y$.

Experimental setup.
Our experimental setup largely follows the approach outlined in Wang et al. (). Specifically, we split each dataset into training, calibration, and testing sets. A Mixture Density Network (MDN) with 10 components is then trained to approximate the conditional distribution $\textup{P}_{Y\mid X}$. For each calibration and test point, we first compute the Gaussian Mixture parameters, forming $\Pi_{Y\mid X}$, and subsequently draw $M=5,20,50$ samples from these distributions, which yield $\mathcal{R}_{z}(x,t)$. This process is repeated across 50 different random splits of each dataset.
Results of the experiments for $M=50$ samples are presented in Figures 3 and 4, additional results are available in Appendix B. In terms of marginal coverage, all methods achieve the target $1-\alpha$ value, except for $\Pi_{Y\mid X}$.

Standard conformal prediction fails to maintain the conditional coverage as expected. We can also observe that PCP consistently struggles with conditional coverage. On all the datasets CP^2-PCP provides valid conditional coverage, while CQR fails on blog and temp. CHR method shows unstable performance not achieving conditional coverage more often than other methods but sometimes providing narrower predictions sets. Additionally, CP^2-PCP significantly outperforms quantile regression-based methods in terms of size of the prediction sets on bike, bio and temp datasets.

Finally, we assess conditional coverage with the help of clustering. We apply HDBSCAN (Campello et al., ; McInnes & Healy, ) method to cluster the test set and then compute coverage within clusters. Results for fb1 dataset are presented in Figure 5. We again observe that CP and PCP do not achieve conditional coverage and CHR and CQR performance is unstable. CP^2-PCP on the other hand maintains valid conditional coverage on all clusters and even on outliers (cluster label -1). Note that these are all outliers combined and they may not lie in the same region of the input space.

Figure: Figure 5: Conditional coverage for different clusters, fb1 dataset. We have used HDBSCAN algorithm with minimum cluster size of 100, min_samples hyper-parameter of 20 and $l_{2}$ metric. Cluster label -1 corresponds to the outliers. Sample size for sampling-based methods was set to 50. Nominal coverage equals $(1-\alpha)=0.9$ and is shown in dashed blacks.
Refer to caption: x7.png

Figure: Figure 6: Marginal coverage for multi-target datasets, 50 replications. Sample size was set to 1000. Nominal coverage equals $(1-\alpha)=0.9$ and is shown by dashed black line.
Refer to caption: x8.png

### 4.3 Real-world Regression Data with Multi-dimensional Targets

We also study CP2 family of algorithms on the multi-target regression problems. Since selecting the threshold $\tau$ for our methods is not dependent on the number of dimensions in $Y$ their application is straightforward. On the other hand, most other methods are inherently one-dimensional thus require the use of the Bonferroni correction (Dunn, ). Each coordinate is treated independently with miscoverage level adjusted to $\alpha/d$, where $d$ is the number of targets. As a result, for quantile regression-based methods prediction sets are rectangular cuboids, formed as a product of the corresponding intervals.

Datasets. We consider open-source multidimensional regression datasets: river flow data rf1 and rf2 (Xioufis et al., ), supply chain management scm1d and scm20d (Xioufis et al., ), indoor localisation indoor (Torres-Sospedra et al., ), GPU computation time sgemm_small(^1^11The full dataset contains 241600 examples. Due to computational constraints we randomly subsample 10000 examples for each replication of our experiment.) (Ballester-Ripoll et al., ).

We use the same metrics as before: marginal coverage and worst-slab coverage. Evaluating the difference in prediction set size is more complex in case of multiple dimensions. Due to computational constraints we perform pairwise comparisons between our methods and selected baselines, measuring approximate areas of 2D projections of the prediction sets  (Wang et al., ). These results can be found on Figure 8. We approximate areas using a grid and fewer samples.

Since our methods naturally extend beyond one dimension, the experimental setup is almost identical. We use the same underlying model for $\textup{P}_{Y\mid X}$, the prediction set is now a union of $d$-dimensional balls of the same radius around the sampled centers. The number of samples is increased to 1000.

Figure: Figure 7: Conditional coverage for multi-target datasets targets, 50 replications. Sample size was set to 1000. Nominal coverage equals $(1-\alpha)=0.9$ and is shown in dashed black. Worst-slab coverage parameter $(1-\delta)=0.1$.
Refer to caption: x9.png

Figure: Figure 8: Average rank of the projected set size. For each pair of targets area of the corresponding 2D projection of the prediction set is calculated. For each test point and each pair of targets methods are ranked. Lower rank is smaller area. This graph shows averaged results of 10 replications.
Refer to caption: x10.png

Results. In Figure 6 we show marginal coverage attained by different algorithms. As expected, naive application of 1D techniques CQR, CQR2 and CHR to multiple outputs produces significant overcover. PCP and CP^2 methods naturally extend to multidimensional targets and provide correct marginal coverage.

In Figure 7 we present the conditional coverage estimates for multi-target datasets. PCP significantly undercovers on rf1, rf2 and scm1d datasets, while CP^2 comes very close to the nominal coverage of $0.9$. In case of CQR, CQR2 and CHR, they still overcover (scm20d, sgemm) or perform comparably to our approach.

Figure 8 shows the aggregated results of the set size comparisons in multidimensional target setting. For each test point and each pair of axes we rank the methods by the area of the projection of the corresponding prediction set. The plot show average rank for each method, aggregated across all axes pairs and replications. Lower rank corresponds to smaller area, which is our goal. For datasets indoor, scm1d and sgemm_small our approach performs better, while also providing sharper conditional covergae, as was shown earlier. On the remaining datasets CP^2performs similarly to the competitors.

## 5 Conclusion

We address the challenge of conditional coverage in CP, and overcome previous negative results by assuming the knowledge of a good estimator of $\textup{P}_{Y\mid X}$.
Our proposed mechanism conformalized the conditional estimator $\Pi_{Y\mid X}$ to ensure marginal validity while maintaining similar conditional coverage guarantees.
Specifically, if experts can provide an accurate conditional estimator, our algorithm CP^2 generates nearly conditionally valid multidimensional prediction sets. This approach offers a practical solution for tackling heteroscedasticity in various machine learning applications.

#### Acknowledgments

Part of this work has been carried out under the auspice of the Lagrange Mathematics and Computing Research Center.
E.M. is Funded by the European Union (ERC, Ocean, 101071601). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## Appendix A Additional Results and Calculations

In this section, we analyze the theoretical results of Section 3.
First, let’s recall the definition of the quantile function for any distribution $\mu_{n}$ living in $\mathbb{R}$.
For any $\alpha\in(0,1)$, the quantile $Q_{1-\alpha}(\mu_{n})$ is defined by

$$ $Q_{1-\alpha}(\mu_{n})=\inf\left\{{t\in\mathbb{R}\colon\mu_{n}((-\infty,t])\geq 1 -\alpha}\right\}.$ $$

Given a measure $\Pi_{Y\mid X=x}$ defined on $\sigma(\mathcal{Y})$, we consider for all $x\in\mathbb{R}^{d}$, $z\in\mathcal{Z}$, the parameters $\tau_{x,z}$ and $\lambda_{x,y,z}$ given by

$$ $\displaystyle\tau_{x,z}=\inf\left\{{\tau\in\mathsf{T}\colon\Pi_{Y\mid X=x}( \mathcal{R}_{z}(x;f_{\tau}(\varphi)))\geq 1-\alpha}\right\},$ (7) $\displaystyle\lambda_{x,y,z}=\inf\left\{{\lambda\in\mathsf{T}\colon y\in \mathcal{R}_{z}(x;\lambda)}\right\},$ $$

where $\varphi$ is chosen as in 2, and by convention we set $\inf\emptyset=\infty$.
We denote by $\delta_{v}$ the Dirac measure at $v\in\mathbb{R}$, and write $\bar{\tau}_{k}=\tau_{X_{k},Z_{k}}$ and $\bar{\lambda}_{k}=\lambda_{X_{k},Y_{k},Z_{k}}$.
In this Appendix, we study the coverage of the prediction set given $\forall(x,z)\in\mathbb{R}\times\mathcal{Z}$ by

$$ $\mathcal{C}_{\alpha}(x)=\mathcal{R}_{z}\left({x;f_{\tau_{x,z}}\bigl{(}Q_{1- \alpha}(\mu_{n})\bigr{)}}\right),$ $$

where the distribution $\mu_{n}$ is defined as

$$ $\mu_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{f_{\bar{\tau}_{k}}^{-1}(\bar{ \lambda}_{k})}+\frac{1}{n+1}\delta_{\infty}.$ $$

The key idea behind the choice of $\bar{\tau}_{k}$ is to ensure that the conditional coverage of the prediction set $\mathcal{C}_{\alpha}(X_{k})$ is approximately $1-\alpha$ when the empirical distribution $\Pi_{Y\mid X=X_{k}}$ is close to $\textup{P}_{Y\mid X=X_{k}}$. In other words, $\bar{\tau}_{k}$ is chosen such that the probability of the observed value $Y_{k}$ given $X_{k}$ falling inside the prediction set $\mathcal{C}_{\alpha}(X_{k})$ is close to $1-\alpha$.
On the other hand, the parameter $\bar{\lambda}_{k}$ is used to ensure that the prediction set $\mathcal{R}_{Z_{k}}(X_{k};\bar{\lambda}_{k})$ contains the observed value $Y_{k}$.
Moreover, note that $\bar{\tau}_{k}$ only depends on the input data $(X_{k},Z_{k})$, while $\bar{\lambda}_{k}$ depends on $(X_{k},Y_{k},Z_{k})$. Thus, the i.i.d. property of $\{(X_{k},Y_{k},Z_{k})\colon k\in[n+1]\}$ ensures that the $\{(\bar{\tau}_{k},\bar{\lambda}_{k})\}_{k=1}^{n+1}$ are also i.i.d.

### A.1 Proof of Theorems 3.1 and 3.2

###### Lemma A.1 .

Assume 1 hold. For any $(x,y,z)\in\mathbb{R}^{d}\times\mathcal{Y}\times\mathcal{Z}$, $\lambda_{x,y,z}$ exists in $\mathsf{T}$, and we have $y\in\mathcal{R}_{z}(x;\lambda_{x,y,z})$.

###### Proof.

Let $(x,y,z)\in\mathbb{R}^{d}\times\mathcal{Y}\times\mathcal{Z}$ be fixed. Since $\cap_{t\in\mathsf{T}}\mathcal{R}_{z}(x;t)=\emptyset$ and $\cup_{t\in\mathsf{T}}\mathcal{R}_{z}(x;t)=\mathcal{Y}$, we deduce the existence of $t_{0}$ and $t_{1}$ such that $y\notin\mathcal{R}_{z}(x;t_{0})$ and $y\in\mathcal{R}_{z}(x;t_{1})$. Therefore, $\{t\in\mathsf{T}\colon y\in\mathcal{R}_{z}(x;t)\}$ is non-empty and lower-bounded by $t_{0}$. Thus, the infimum $\lambda_{x,y,z}$ exists.
Now, let’s prove that $y\in\mathcal{R}_{z}(x;\lambda_{x,y,z})$.
Since $\lambda_{x,y,z}=\inf\{t\in\mathsf{T}\colon y\in\mathcal{R}_{z}(x;t)\}$, we deduce the existence of a decreasing sequence $\{\lambda_{n}\}_{n\in\mathbb{N}}$ such that $y\in\mathcal{R}_{z}(x;\lambda_{n})$ and $\lim_{n\to\infty}\lambda_{n}=\lambda_{x,y,z}$. By definition of $\{\lambda_{n}\}_{n\in\mathbb{N}}$, we have $y\in\cap_{n\in\mathbb{N}}\mathcal{R}_{z}(x;\lambda_{n})$. However, using 1, remark that

$$ $\displaystyle\cap_{n\in\mathbb{N}}\mathcal{R}_{z}(x;\lambda_{n})$ $\displaystyle=\cap_{n\in\mathbb{N}}\cap_{t>\lambda_{n}}\mathcal{R}_{z}(x;t)$ $\displaystyle=\cap_{t>\lim\limits_{n\to\infty}\lambda_{n}}\mathcal{R}_{z}(x;t)$ $\displaystyle=\cap_{t>\lambda_{x,y,z}}\mathcal{R}_{z}(x;t)=\mathcal{R}_{z}(x; \lambda_{x,y,z}).$ $$

Since $y\in\cap_{n\in\mathbb{N}}\mathcal{R}_{z}(x;\lambda_{n})$, it implies that $y\in\mathcal{R}_{z}(x;\lambda_{x,y,z})$.
∎

We will now present the proof for Theorem 3.1, which establishes the marginal validity of our proposed method.

###### Theorem A.2 .

Assume 1-2 hold, if $\{f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})\}_{k=1}^{n+1}$ are almost surely distinct, then it follows

$$ $1-\alpha\leq\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n +1})}\right)<1-\alpha+\frac{1}{n+1}.$ (8) $$

###### Proof.

Using Lemma A.1, we have

$$ $\displaystyle\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{ n+1})}\right)$ $\displaystyle=\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{R}_{Z_{n+1}} \left({X_{n+1},f_{\bar{\tau}_{n+1}}(Q_{1-\alpha}(\mu_{n}))}\right)}\right)$ $\displaystyle=\mathbb{P}^{\mathcal{T}}\left({\lambda_{n+1}\leq f_{\bar{\tau}_{ n+1}}(Q_{1-\alpha}(\mu_{n}))}\right).$ $$

Since $\lambda\mapsto f_{\bar{\tau}_{n+1}}(\lambda)$ is increasing by 2, we deduce that

$$ $\mathbb{P}^{\mathcal{T}}\left({\lambda_{n+1}\leq f_{\bar{\tau}_{n+1}}(Q_{1- \alpha}(\mu_{n}))}\right)=\mathbb{P}^{\mathcal{T}}\left({f_{\bar{\tau}_{n+1}}^ {-1}(\lambda_{n+1})\leq Q_{1-\alpha}(\mu_{n})}\right).$ $$

Denote by $V_{k}=f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})$, the exchangeability of the data $\{(X_{k},Y_{k},Z_{k})\colon k\in[n+1]\}$ implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({V_{n+1}\leq Q_{1-\alpha}\left({\sum_{k=1}^{n} \frac{\delta_{V_{k}}}{n+1}+\frac{\delta_{\infty}}{n+1}}\right)}\right)=\mathbb {P}^{\mathcal{T}}\left({V_{n+1}\leq Q_{1-\alpha}\left({\sum_{k=1}^{n+1}\frac{ \delta_{V_{k}}}{n+1}}\right)}\right)\\ =\frac{1}{n+1}\sum_{k=1}^{n+1}\mathbb{E}^{\mathcal{T}}\left[{\mathds{1}_{V_{k} }\leq Q_{1-\alpha}\left({\frac{1}{n+1}\sum_{k=1}^{n+1}\delta_{V_{k}}}\right)} \right]\\ =\mathbb{E}^{\mathcal{T}}\left[{\mathbb{E}^{\mathcal{T}}\left[{\mathds{1}_{V_{ I}}\leq Q_{1-\alpha}\left({\frac{1}{n+1}\sum_{k=1}^{n+1}\delta_{V_{k}}}\right) \,\bigg{|}\,V_{1},\ldots,V_{n+1}}\right]}\right],$ $$

where $I\sim\mathcal{U}nif(1,\ldots,n+1)$. Therefore, the definition of the quantile function implies the lower bound in (8). Moreover, if there are no ties between the $\{V_{k}\}_{k=1}^{n+1}$, then

$$ $\mathbb{P}^{\mathcal{T}}\left({f_{\bar{\tau}_{n+1}}^{-1}(\lambda_{n+1})\leq Q_ {1-\alpha}(\mu_{n})}\right)<1-\alpha+\frac{1}{n+1}.$ $$

∎

The following lemma provides conditions under which $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))\geq 1-\alpha$.

###### Lemma A.3 .

Assume 1-2 hold, and let $\alpha\in(0,1$), $x\in\mathbb{R}^{d}$, $z\in\mathcal{Z}$.
If $\Pi_{Y\mid X=x}$ is a probability measure, then $\tau_{x,z}$ is defined in $\mathsf{T}$ and $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))\geq 1-\alpha$.

###### Proof.

Let $x\in\mathbb{R}^{d}$ be such that $\Pi_{Y\mid X=x}$ is a probability measure, and fix $z\in\mathcal{Z}$. Since $\tau\mapsto f_{\tau}(\varphi)$ is increasing and bijective by 2, we have

$$ $\displaystyle\sup_{\tau\in\mathsf{T}}\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau }(\varphi)))$ $\displaystyle=\Pi_{Y\mid X=x}\left({\cup_{\tau\in\mathsf{T}}\mathcal{R}_{z}(x; f_{\tau}(\varphi))}\right)$ $\displaystyle=\Pi_{Y\mid X=x}\left({\cup_{t\in\mathsf{T}}\mathcal{R}_{z}(x;t)} \right)=1.$ $$

The previous equality shows the existence of $\tau\in\mathsf{T}$ such that $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau}(\varphi)))\geq 1-\alpha$. Therefore $\{\tau\in\mathsf{T}\colon\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau}(\varphi)))
\geq 1-\alpha\}$ is non-empty.
This proves the existence of $\tau_{x,z}=\inf\{\tau\in\mathsf{T}\colon\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{
\tau}(\varphi)))\geq 1-\alpha\}$ in $\mathsf{T}\cup\{-\infty\}$.
Moreover, $\tau_{x,z}>-\infty$, otherwise we would have

$$ $1-\alpha\leq\inf_{\tau\in\mathsf{T}}\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau} (\varphi)))=\Pi_{Y\mid X=x}\left({\cap_{t\in\mathsf{T}}\mathcal{R}_{z}(x;t)} \right)=0.$ $$

Therefore, we deduce that $\tau_{x,z}\in\mathsf{T}$.
Lastly, remark that

$$ $\displaystyle\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))$ $\displaystyle=\Pi_{Y\mid X=x}(\cap_{\tau>\tau_{x,z}}\mathcal{R}_{z}(x;f_{\tau} (\varphi)))$ $\displaystyle=\inf_{\tau>\tau_{x,z}}\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau} (\varphi)))\geq 1-\alpha.$ $$

∎

Now, we prove Theorem 3.2. This result guarantees that the conditional confidence intervals constructed by our method approximately satisfy the desired coverage of $1-\alpha$.
Given $(x,y)\in\mathbb{R}^{d}\times\mathcal{Z}$, let’s introduce

$$ $\displaystyle p_{n+1}^{(x,z)}=\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\mu_ {n})<f_{\tau_{x,z}}^{-1}\left({\lambda_{x,Y_{n+1},z}}\right)\leq\varphi\,|\,X_ {n+1}=x,Z_{n+1}=z}\right),$ $\displaystyle q_{n+1}^{(x,z)}=\mathbb{P}^{\mathcal{T}}\left({\varphi<f_{\tau_{ x,z}}^{-1}\left({\lambda_{x,Y_{n+1},z}}\right)\leq Q_{1-\alpha}(\mu_{n})\,|\,X _{n+1}=x,Z_{n+1}=z}\right).$ $$

###### Theorem A.4 .

Assume 1-2 hold, let $x\in\mathbb{R}^{d}$ be such that $\Pi_{Y\mid X=x}$ is a probability measure. For any $z\in\mathcal{Z}$, it follows that

$$ $1-\alpha-\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})-p_{n +1}^{(x,z)}\leq\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X _{n+1})\mid X_{n+1}=x,Z_{n+1}=z}\right)\\ \leq\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))+\mathrm{d}_{ \mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})+q_{n+1}^{(x,z)}.$ $$

###### Proof.

First, recall that $\mathcal{C}_{\alpha}(x)$ is given in (4), and $\lambda_{x,Y_{n+1},z}$ is defined in (7). Applying Lemma A.1, we know that $\lambda_{x,Y_{n+1},z}$ is defined in $\mathsf{T}$, and also that $Y_{n+1}\in\mathcal{R}_{z}(x;\lambda_{x,Y_{n+1},z})$. Hence, it holds

$$ $\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{ n+1}=x,Z_{n+1}=z}\right)\\ =\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{R}_{z}\left({x;f_{\tau_{x,z }}(Q_{1-\alpha}(\mu_{n}))}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ =\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left( {Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$ (9) $$

Let’s introduce the term $\mathbb{P}^{\mathcal{T}}(\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(\varphi)\,|
\,X_{n+1}=x,Z_{n+1}=z)$ as follows:

$$ $\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left({ Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ =\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left( {Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ \pm\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}( \varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$ (10) $$

Now, we will control the difference between the two terms of the previous equation.
Let $A$ and $B$ be defined as

$$ $\displaystyle A=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{x,z}}^{-1}(\lambda_{x, Y_{n+1},z})\leq Q_{1-\alpha}(\mu_{n})<\varphi\,|\,X_{n+1}=x,Z_{n+1}=z}\right),$ $\displaystyle B=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{x,z}}^{-1}(\lambda_{x, Y_{n+1},z})\leq\varphi\leq Q_{1-\alpha}(\mu_{n})\,|\,X_{n+1}=x,Z_{n+1}=z} \right).$ $$

We have

$$ $\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left({ Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ =A+B+\mathbb{P}^{\mathcal{T}}\left({\varphi<f_{\tau_{x,z}}^{-1}\left({\lambda_ {x,Y_{n+1},z}}\right)\leq Q_{1-\alpha}(\mu_{n})\,|\,X_{n+1}=x,Z_{n+1}=z}\right),$ $$

and also

$$ $\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}( \varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ =A+B+\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\mu_{n})<f_{\tau_{x,z}}^{-1} \left({\lambda_{x,Y_{n+1},z}}\right)\leq\varphi\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$ $$

Therefore, the difference between the terms introduced in (10) can be rewritten as

$$ $\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left({ Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ -\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}( \varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ =\mathbb{P}^{\mathcal{T}}\left({\varphi<f_{\tau_{x,z}}^{-1}\left({\lambda_{x,Y _{n+1},z}}\right)\leq Q_{1-\alpha}(\mu_{n})\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\ -\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\mu_{n})<f_{\tau_{x,z}}^{-1}\left ({\lambda_{x,Y_{n+1},z}}\right)\leq\varphi\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$ (11) $$

- 1.
By definition of the total variation distance, we have
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
\geq\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,\hat{Y}_{n+1},z}\leq f_{\tau_{x,
z}}(\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)-\mathrm{d}_{\mathrm{TV}}(\textup{
P}_{Y\mid X=x};\Pi_{Y\mid X=x}).$
Moreover, Lemma A.3 implies that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,\hat{Y}_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
=\mathbb{P}^{\mathcal{T}}\left({\hat{Y}_{n+1}\in\left\{{y\in\mathcal{Y}\colon
\lambda_{x,y,z}\leq f_{\tau_{x,z}}(\varphi)}\right\}\,|\,X_{n+1}=x,Z_{n+1}=z}
\right)\\
=\mathbb{P}^{\mathcal{T}}\left({\hat{Y}_{n+1}\in\mathcal{R}_{z}(x;f_{\tau_{x,z
}}(\varphi))\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
=\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))\geq 1-\alpha.$
Therefore, we deduce that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\geq 1-\alpha-\mathrm{d}_{\mathrm{TV}}
(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x}).$
Combining the previous result with (10) and (11) shows that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left({
Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\geq 1-\alpha-
\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})\\
-\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\mu_{n})<f_{\tau_{x,z}}^{-1}\left
({\lambda_{x,Y_{n+1},z}}\right)\leq\varphi\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$
Finally, using (9) gives a lower bound on $\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})\mid X_{
n+1}=x,Z_{n+1}=z}\right)$.
- 2.
By definition of the total variation distance, we have
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
\leq\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,\hat{Y}_{n+1},z}\leq f_{\tau_{x,
z}}(\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)+\mathrm{d}_{\mathrm{TV}}(\textup{
P}_{Y\mid X=x};\Pi_{Y\mid X=x}).$
Moreover, Lemma A.3 implies that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,\hat{Y}_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)=\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{
\tau_{x,z}}(\varphi))).$
Therefore, we deduce that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(
\varphi)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
\leq\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))+\mathrm{d}_{
\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x}).$
Finally, combining the previous result with (10) and (11) shows that
$\mathbb{P}^{\mathcal{T}}\left({\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}\left({
Q_{1-\alpha}(\mu_{n})}\right)\,|\,X_{n+1}=x,Z_{n+1}=z}\right)\\
\leq\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))+\mathrm{d}_{
\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})+q_{n+1}^{(x,z)}.$

∎

###### Lemma A.1 .

###### Proof.

###### Theorem A.2 .

###### Proof.

###### Lemma A.3 .

###### Proof.

###### Theorem A.4 .

###### Proof.

### A.2 Bound on p n + 1 ( x , z ) superscript subscript 𝑝 𝑛 1 𝑥 𝑧 p_{n+1}^{(x,z)} italic_p start_POSTSUBSCRIPT italic_n + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_x , italic_z ) end_POSTSUPERSCRIPT and q n + 1 ( x , z ) superscript subscript 𝑞 𝑛 1 𝑥 𝑧 q_{n+1}^{(x,z)} italic_q start_POSTSUBSCRIPT italic_n + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_x , italic_z ) end_POSTSUPERSCRIPT

The objective of this section is to study the conditional guarantee obtained in Theorem A.4. Under some assumptions, we have demonstrated that the conditional coverage is controlled as follows:

$$ $1-\alpha-\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})-p_{n +1}^{(x,z)}\leq\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X _{n+1})\mid X_{n+1}=x,Z_{n+1}=z}\right)\\ \leq\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))+\mathrm{d}_{ \mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})+q_{n+1}^{(x,z)},$ $$

In the following, we consider the cumulative density functions $F\colon t\mapsto\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})
\leq t)$ and $\hat{F}\colon t\mapsto\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,
\hat{Y},Z})\leq t)$, where $(X,Y,Z)\sim\textup{P}_{X}\times\textup{P}_{Y\mid X}\times\Pi_{Z\mid X}$ and $(X,\hat{Y},Z)\sim\textup{P}_{X}\times\Pi_{Y\mid X}\times\Pi_{Z\mid X}$.
We denote by $\mu$ and $\hat{\mu}$ the law of the random variables $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ and $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$.
Moreover, recall that $\mu_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{f_{\bar{\tau}_{X_{k}}}^{-1}(\bar{
\lambda}_{k})}+\frac{1}{n+1}\delta_{\infty}$.
Note, the quantile $Q_{1-\alpha}(\mu_{n})$ is an order statistic with a known distribution that converges to the true quantile $Q_{1-\alpha}(\mu)$. The quantile is defined for any $t\in(0,1)$ by

$$ $Q_{t}(\nu)=\inf\{u\in\mathbb{R}\colon\nu((-\infty,u])\geq t\},\qquad\text{ where $\nu\in\{\mu,\mu_{n},\hat{\mu}$\}}.$ (12) $$

###### Theorem A.5 .

Assume 1-2 hold, and let $x\in\mathbb{R}^{d}$ be such that $\Pi_{Y\mid X=x}$ is a probability measure.
For any $\epsilon\in[0,1-\alpha)$, if $p_{\epsilon}=\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})<Q_{
1-\alpha-\epsilon}(\mu))\leq 1-\alpha$, then it follows that

$$ $p_{n+1}^{(x,z)}\leq\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha-\epsilon}(\mu)<f _{\tau_{x,y}}^{-1}(\lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha}(\hat{\mu})\,|\,X_{n +1}=x,Z_{n+1}=z}\right)\\ +\exp\left({-np_{\epsilon}(1-p_{\epsilon})h\left({\frac{1-\alpha-p_{\epsilon}} {p_{\epsilon}(1-p_{\epsilon})}}\right)}\right),$ $$

where $h:u\mapsto(1+u)\log(1+u)-u$.

###### Proof.

Let $\epsilon\in[0,1-\alpha)$, $x\in\mathbb{R}^{d}$, and consider

$$ $\displaystyle A=\left\{{Q_{1-\alpha}(\mu_{n})<Q_{1-\alpha-\epsilon}(\mu)} \right\},$ $\displaystyle B_{x,z}=\left\{{y\in\mathcal{Y}\colon f_{\tau_{x,z}}(Q_{1-\alpha -\epsilon}(\mu))<\lambda_{x,y,z}\leq f_{\tau_{x,z}}(\varphi)}\right\}.$ $$

We have

$$ $\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{x,z}}\left({Q_{1-\alpha}(\mu_{n})} \right)<\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(\varphi)\,|\,X_{n+1}=x,Z_{n+1 }=z}\right)\\ \leq\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)+\mathbb{P }^{\mathcal{T}}\left({Y_{n+1}\in B_{x,z}\mid X_{n+1}=x,Z_{n+1}=z}\right).$ $$

Now, let’s upper bound the first term of the right-hand side equation.
First, remark that

$$ $\left\{{Q_{1-\alpha}(\mu_{n})<Q_{1-\alpha-\epsilon}(\mu)}\right\} \Leftrightarrow\left\{{\frac{1}{n+1}\sum_{k=1}^{n}\mathds{1}_{f_{\bar{\tau}_{k }}^{-1}(\bar{\lambda}_{k})<Q_{1-\alpha-\epsilon}(\mu)}\geq 1-\alpha}\right\}.$ $$

Thus, we deduce that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\mathbb{P} ^{\mathcal{T}}\left({\sum_{k=1}^{n}\mathds{1}_{f_{\bar{\tau}_{k}}^{-1}(\bar{ \lambda}_{k})<Q_{1-\alpha-\epsilon}(\mu)}\geq(n+1)(1-\alpha)}\right).$ $$

Recall that $p_{\epsilon}=\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})<Q_{
1-\alpha-\epsilon}(\mu))$, and also that we assume $p_{\epsilon}\leq 1-\alpha$.
Therefore, the Bennett’s inequality (Boucheron et al., , Theorem 2) implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\exp\left( {-np_{\epsilon}(1-p_{\epsilon})\,h\left({\frac{(n+1)(1-\alpha)-np_{\epsilon}}{ np_{\epsilon}(1-p_{\epsilon})}}\right)}\right),$ (13) $$

where $h:u\mapsto(1+u)\log(1+u)-u$.
Moreover, define

$$ $\displaystyle u_{\epsilon}=\frac{1-\alpha-p_{\epsilon}}{p_{\epsilon}(1-p_{ \epsilon})},$ $\displaystyle\tilde{u}_{\epsilon}=\frac{(n+1)(1-\alpha)-np_{\epsilon}}{np_{ \epsilon}(1-p_{\epsilon})}.$ $$

We have $\tilde{u}_{\epsilon}\leq u_{\epsilon}$, from the increasing property of $h$ it follows that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\exp\left( {-np_{\epsilon}(1-p_{\epsilon})h(u_{\epsilon})}\right).$ $$

Furthermore, the definition of $B_{x,z}$ gives

$$ $\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in B_{x,z}\mid X_{n+1}=x,Z_{n+1}=z} \right)\\ =\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{x,z}}(Q_{1-\alpha-\epsilon}(\mu))< \lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(\varphi)\mid X_{n+1}=x,Z_{n+1}=z} \right).$ $$

Moreover, for any $t\in(-\infty,\varphi)$, we have

$$ $\displaystyle\hat{F}(t)$ $\displaystyle=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X, \hat{Y},Z})\leq t}\right)$ $\displaystyle=\int\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{ X,\hat{Y},Z})\leq t\,\Big{|}\,X=x,Z=z}\right)\bar{\Pi}_{Z\mid X=x}(\mathrm{d}z )\,\textup{P}_{X}(\mathrm{d}x)$ $\displaystyle=\int\mathbb{P}^{\mathcal{T}}\left({\hat{Y}\in\mathcal{R}\left({x ,f_{\tau_{z,z}}(t)}\right)\,\Big{|}\,X=x,Z=z}\right)\bar{\Pi}_{Z\mid X=x}( \mathrm{d}z)\,\textup{P}_{X}(\mathrm{d}x).$ $$

Using 2, the bijective property of $\tau\mapsto f_{\tau}(\varphi)$ implies the existence of $\nu\in\mathsf{T}$, such that $f_{\nu}(\varphi)=f_{\tau_{z,z}}(t)$. Note that, $\nu<\tau_{x,z}$ otherwise it would lead to $f_{\nu}(\varphi)\geq f_{\tau_{x,z}}(\varphi)>f_{\tau_{x,z}}(t)$. The definition of $\tau_{x,z}$ shows that

$$ $\mathbb{P}^{\mathcal{T}}\left({\hat{Y}\in\mathcal{R}\left({x,f_{\nu}(\varphi)} \right)\,\Big{|}\,X=x,Z=z}\right)<1-\alpha.$ $$

Therefore, we deduce that $Q_{1-\alpha}(\hat{\mu})\geq\varphi$, and we can conclude that

$$ $\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in B_{x,z}\mid X_{n+1}=x,Z_{n+1}=z} \right)\\ \leq\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha-\epsilon}(\mu)<f_{\tau_{x,y}}^{ -1}(\lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha}(\hat{\mu})\mid X_{n+1}=x,Z_{n+1}=z }\right).$ (14) $$

Finally, combining (13) and (14) concludes the proof.
∎

Given $\alpha\in(0,1)$, define the threshold

$$ $\epsilon_{n}=\sqrt{\frac{8\alpha(1-\alpha)\log n}{n}}.$ (15) $$

###### Lemma A.6 .

If the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ is continuous, then for all $\epsilon\in[0,1-\alpha)$, we have $p_{\epsilon}=\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})<Q_{
1-\alpha-\epsilon}(\mu))=1-\alpha-\epsilon$.
Moreover, if $\epsilon_{n}\leq\frac{\alpha(1-\alpha)}{8}$, then it follows

$$ $\exp\left({-np_{\epsilon_{n}}(1-p_{\epsilon_{n}})h\left({\frac{1-\alpha-p_{ \epsilon_{n}}}{p_{\epsilon_{n}}(1-p_{\epsilon_{n}})}}\right)}\right)\leq\frac{ 1}{n},$ $$

where $h:u\mapsto(1+u)\log(1+u)-u$.

###### Proof.

First, recall that $Q_{1-\alpha-\epsilon}(\mu)$ is defined in (12).
If the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ is continuous, then we have

$$ $1-\alpha-\epsilon\leq F(Q_{1-\alpha-\epsilon}(\mu))=\sup_{\delta>0}F(Q_{1- \alpha-\epsilon}(\mu)-\delta)\\ \leq\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})<Q_{1- \alpha-\epsilon}(\mu)}\right)=p_{\epsilon}\leq 1-\alpha-\epsilon.$ $$

Therefore, we deduce that $p_{\epsilon}=1-\alpha-\epsilon$.
Let’s denote

$$ $\displaystyle\delta_{n}$ $\displaystyle=(n+1)(1-\alpha)-np_{\epsilon_{n}},$ $\displaystyle u_{n}$ $\displaystyle=\frac{(n+1)(1-\alpha)-np_{\epsilon_{n}}}{np_{\epsilon_{n}}(1-p_{ \epsilon_{n}})}.$ $$

For any $u\geq 0$, remark that $\log(1+u)\geq u-u^{2}/2$. Thus, we deduce

$$ $\displaystyle np_{\epsilon_{n}}(1-p_{\epsilon_{n}})\,h\left({u_{n}}\right)$ $\displaystyle\geq\delta_{n}\frac{(1+u_{n})\log(1+u_{n})-u_{n}}{u_{n}}$ $\displaystyle\geq\delta_{n}\frac{u_{n}(1-u_{n})}{2}.$ (16) $$

Now, let’s show that $u_{n}\leq 1/4$. We have

$$ $\displaystyle u_{n}$ $\displaystyle=\frac{(n+1)(1-\alpha)-np_{\epsilon_{n}}}{np_{\epsilon_{n}}(1-p_{ \epsilon_{n}})}$ $\displaystyle=\frac{1-\alpha}{np_{\epsilon_{n}}(1-p_{\epsilon_{n}})}+\frac{1- \alpha-p_{\epsilon_{n}}}{p_{\epsilon_{n}}(1-p_{\epsilon_{n}})}$ $\displaystyle=\frac{1-\alpha}{n(\alpha+\epsilon_{n})(1-\alpha-\epsilon_{n})}+ \frac{\epsilon_{n}}{(\alpha+\epsilon_{n})(1-\alpha-\epsilon_{n})}.$ $$

Therefore, $u_{n}\leq 1/4$ if and only if

$$ $\frac{1-\alpha}{n}+\epsilon_{n}\leq\frac{(\alpha+\epsilon_{n})(1-\alpha- \epsilon_{n})}{4}.$ $$

The function $\epsilon\in[0,1/2-\alpha]\mapsto(\alpha+\epsilon)(1-\alpha-\epsilon)$ is increasing. Since $\epsilon_{n}\leq\alpha(1-\alpha)/8\leq 1/2-\alpha$, it is sufficient to prove that

$$ $\frac{1-\alpha}{n}+\epsilon_{n}\leq\frac{\alpha(1-\alpha)}{4}.$ $$

Since $\epsilon_{n}\leq\alpha(1-\alpha)/8$, we just need to show that

$$ $\frac{1-\alpha}{n}\leq\frac{\alpha(1-\alpha)}{8},\qquad\text{i.e.,}\qquad\frac {8\alpha(1-\alpha)}{n}\leq\alpha^{2}(1-\alpha).$ (17) $$

Again, using the fact that $\epsilon_{n}\leq\alpha(1-\alpha)/8$, we deduce that

$$ $\frac{8\alpha(1-\alpha)}{n}=\frac{\epsilon_{n}^{2}}{\log n}\leq\frac{\alpha^{2 }(1-\alpha)^{2}}{8\log n}=\alpha^{2}(1-\alpha)\times\frac{(1-\alpha)}{8\log n}.$ $$

Since $\frac{(1-\alpha)}{8\log n}\leq 1$, we deduce that (17) holds. This concludes that $u_{n}\leq 1/4$.
Moreover, for any $u\in[0,0.25]$, we have

$$ $\delta_{n}\frac{u(1-u)}{2}\geq\frac{u\delta_{n}}{4}.$ $$

Plugging the previous line in (16) implies that

$$ $\displaystyle\exp\left({-np_{\epsilon_{n}}(1-p_{\epsilon_{n}})\,h\left({u_{n}} \right)}\right)$ $\displaystyle\leq\exp\left({-\frac{\left[{(n+1)(1-\alpha)-np_{\epsilon_{n}}} \right]^{2}}{4np_{\epsilon_{n}}(1-p_{\epsilon_{n}})}}\right)$ $\displaystyle\leq\exp\left({-\frac{\left({1-\alpha+n\epsilon_{n}}\right)^{2}}{ 4n(\alpha+\epsilon_{n})(1-\alpha-\epsilon_{n})}}\right)$ $\displaystyle\leq\exp\left({-\frac{n\epsilon_{n}^{2}}{4(\alpha+\epsilon_{n})(1 -\alpha-\epsilon_{n})}}\right).$ (18) $$

Lastly, since $\epsilon_{n}\leq\alpha$, it follows that

$$ $\displaystyle\frac{n\epsilon_{n}^{2}}{4(\alpha+\epsilon_{n})(1-\alpha-\epsilon _{n})}=\frac{2\alpha(1-\alpha)\log n}{(\alpha+\epsilon_{n})(1-\alpha-\epsilon_ {n})}\geq\log n.$ $$

Combining the previous line with (18) completes the proof.
∎

For any $\epsilon\in[0,\alpha)$, define

$$ $q_{\epsilon}=\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})<Q_{ 1-\alpha+\epsilon}(\mu)).$ $$

###### Theorem A.7 .

Assume 1-2 hold, and let $x\in\mathbb{R}^{d}$ be such that $\Pi_{Y\mid X=x}$ is a probability measure.
If the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ is continuous and $n^{-1}\log n\leq 8^{-3}\alpha(1-\alpha)$, then, it holds

$$ $q_{n+1}^{(x,z)}\leq\frac{1}{n}+\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}( \hat{\mu})<f_{\tau_{x,y}}^{-1}(\lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha+\epsilon _{n}}(\mu)\,|\,X_{n+1}=x,Z_{n+1}=z}\right),$ (19) $$

where $\epsilon_{n}$ is defined in (15).

###### Proof.

Let’s consider

$$ $\displaystyle A=\left\{{Q_{1-\alpha+\epsilon_{n}}(\mu)<Q_{1-\alpha}(\mu_{n})} \right\},$ $\displaystyle B_{x,z}=\left\{{y\in\mathcal{Y}\colon f_{\tau_{x,z}}(Q_{1-\alpha }(\hat{\mu}))<\lambda_{x,y,z}\leq f_{\tau_{x,z}}(Q_{1-\alpha+\epsilon_{n}}(\mu ))}\right\}.$ $$

We have

$$ $\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{x,z}}\left({Q_{1-\alpha}(\hat{\mu})} \right)<\lambda_{x,Y_{n+1},z}\leq f_{\tau_{x,z}}(Q_{1-\alpha}(\mu_{n}))\,|\,X_ {n+1}=x,Z_{n+1}=z}\right)\\ \leq\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)+\mathbb{P }^{\mathcal{T}}\left({Y_{n+1}\in B_{x,z}\mid X_{n+1}=x,Z_{n+1}=z}\right).$ (20) $$

Now, let’s upper bound the first term of the right-hand side equation.
First, remark that

$$ $\left\{{Q_{1-\alpha+\epsilon_{n}}(\mu)<Q_{1-\alpha}(\mu_{n})}\right\} \Leftrightarrow\left\{{\frac{1}{n+1}\sum_{k=1}^{n}\mathds{1}_{f_{\bar{\tau}_{k }}^{-1}(\bar{\lambda}_{k})<Q_{1-\alpha+\epsilon_{n}}(\mu)}<1-\alpha}\right\}.$ $$

Thus, we deduce that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\mathbb{P} ^{\mathcal{T}}\left({\sum_{k=1}^{n}\mathds{1}_{f_{\bar{\tau}_{k}}^{-1}(\bar{ \lambda}_{k})<Q_{1-\alpha+\epsilon_{n}}(\mu_{n})}<(n+1)(1-\alpha)}\right).$ $$

Recall that $q_{\epsilon_{n}}=\mathbb{P}^{\mathcal{T}}(f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})
<Q_{1-\alpha+\epsilon_{n}}(\mu))$, and also that $q_{\epsilon_{n}}<1$ since the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ is continuous with $1-\alpha+\epsilon_{n}<1$.
Therefore, the Bennett’s inequality (Boucheron et al., , Theorem 2) implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\exp\left( {-nq_{\epsilon_{n}}(1-q_{\epsilon_{n}})\,h\left({\frac{(n+1)(1-\alpha)-nq_{ \epsilon_{n}}}{nq_{\epsilon_{n}}(1-q_{\epsilon_{n}})}}\right)}\right),$ $$

where $h:u\mapsto(1+u)\log(1+u)-u$.
Moreover, define

$$ $\displaystyle u_{\epsilon_{n}}=\frac{1-\alpha-q_{\epsilon_{n}}}{q_{\epsilon_{n }}(1-q_{\epsilon_{n}})},$ $\displaystyle\tilde{u}_{\epsilon_{n}}=\frac{(n+1)(1-\alpha)-nq_{\epsilon_{n}}} {nq_{\epsilon_{n}}(1-q_{\epsilon_{n}})}.$ $$

We have $\tilde{u}_{\epsilon_{n}}\leq u_{\epsilon_{n}}$, from the increasing property of $h$ combined with Lemma A.6, it follows that

$$ $\mathbb{P}^{\mathcal{T}}\left({A\mid X_{n+1}=x,Z_{n+1}=z}\right)\leq\exp\left( {-nq_{\epsilon_{n}}(1-q_{\epsilon_{n}})h(u_{\epsilon_{n}})}\right)\leq n^{-1}.$ $$

The previous inequality combined with (20) concludes the proof.
∎

###### Theorem A.5 .

###### Proof.

###### Lemma A.6 .

###### Proof.

###### Theorem A.7 .

###### Proof.

### A.3 Proof of Theorem 3.3

###### Theorem A.8 .

Assume 1-2 and suppose the distributions of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ and $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$ are continuous.
For any $\alpha\in(0,1)$ and $\rho>0$, it holds

$$ $\mathbb{P}^{\mathcal{T}}\left({\left\lvert{\mathbb{P}^{\mathcal{T}}\left({Y_{n +1}\in\mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)-1+\alpha} \right\rvert>\rho}\right)\\ \leq\frac{2n^{-1}+\sqrt{128\alpha(1-\alpha)n^{-1}\log n}+4\mathrm{d}_{\mathrm{ TV}}(\textup{P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X})}{\rho}.$ $$

###### Proof.

Let $\rho>0$ be fixed.
Applying Theorem 3.2, we obtain that

$$ $1-\alpha-\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})-p_{n +1}^{(x,z)}\leq\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X _{n+1})\mid X_{n+1}=x,Z_{n+1}=z}\right)\\ \leq\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))+\mathrm{d}_{ \mathrm{TV}}(\textup{P}_{Y\mid X=x};\Pi_{Y\mid X=x})+q_{n+1}^{(x,z)}.$ (21) $$

###### Theorem A.8 .

###### Proof.

##### Step 1: Lower bound.

Using the Markov’s inequality implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in \mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)<1-\alpha-\rho}\right )\\ \leq\mathbb{P}^{\mathcal{T}}\left({\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X };\Pi_{Y\mid X})+p_{n+1}^{(X,Z)}<\rho}\right)\leq\frac{\mathbb{E}^{\mathcal{T} }\big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X};\Pi_{Y\mid X})}\big{]}+ \mathbb{E}^{\mathcal{T}}\big{[}{p_{n+1}^{(X,Z)}}\big{]}}{\rho}.$ (22) $$

Moreover, using Theorem A.5 with $\Phi(\epsilon)=\epsilon[(u_{\epsilon}^{-1}-1)\log(1+u_{\epsilon})-1]$ and $u_{\epsilon}=\epsilon(\alpha+\epsilon)^{-1}(1-\alpha-\epsilon)$, it holds

$$ $\mathbb{E}^{\mathcal{T}}\big{[}{p_{n+1}^{(X,Z)}}\big{]}=\mathbb{P}^{\mathcal{T }}\left({Q_{1-\alpha}(\mu_{n})<f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})\leq\varphi }\right)\\ \leq\exp\left({-n\Phi(\epsilon)}\right)+\mathbb{P}^{\mathcal{T}}\left({Q_{1- \alpha-\epsilon}(\mu)<f_{\bar{\tau}_{n+1}}^{-1}(\bar{\lambda}_{n+1})\leq Q_{1- \alpha}(\hat{\mu})\,\big{|}\,X_{n+1}=x,Z_{n+1}=z}\right).$ $$

By Lemma A.6, if $n^{-1}\log n\leq 8^{-3}\alpha(1-\alpha)$, then, setting $\epsilon_{n}=\sqrt{8\alpha(1-\alpha)n^{-1}\log n}$ ensures that $\exp(-n\Phi(\epsilon_{n}))\leq n^{-1}$. We assume in the following that $n^{-1}\log n\leq 8^{-3}\alpha(1-\alpha)$, because, if it not the case, the final upper bound obtained at the end of the proof is still valid. Thus, we get

$$ $\mathbb{E}^{\mathcal{T}}\big{[}{p_{n+1}^{(X,Z)}}\big{]}\leq n^{-1}+\mathbb{P}^ {\mathcal{T}}\left({Q_{1-\alpha-\epsilon_{n}}(\mu)<f_{\bar{\tau}_{n+1}}^{-1}( \bar{\lambda}_{n+1})\leq Q_{1-\alpha}(\hat{\mu})}\right).$ (23) $$

Let’s define $\bar{\gamma}$ by

$$ $\bar{\gamma}=\min(1,1-\alpha+\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y};\textup {P}_{X}\times\Pi_{Y\mid X})).$ $$

We now show that $Q_{1-\alpha}(\hat{\mu})\leq Q_{\bar{\gamma}}(\mu)$.
By continuity of the cumulative density function of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$, we have

$$ $\displaystyle 1-\alpha$ $\displaystyle=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X, \hat{Y},Z})\leq Q_{1-\alpha}(\hat{\mu})}\right)$ $\displaystyle\geq\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X ,Y,Z})\leq Q_{1-\alpha}(\hat{\mu})}\right)-\mathrm{d}_{\mathrm{TV}}(\textup{P} _{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Hence, it follows that

$$ $\displaystyle\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z })\leq Q_{1-\alpha}(\hat{\mu})}\right)\leq\bar{\gamma}\leq\mathbb{P}^{\mathcal {T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})\leq Q_{\bar{\gamma}}(\mu)} \right).$ $$

Thus, the previous line implies that $Q_{1-\alpha}(\hat{\mu})\leq Q_{\bar{\gamma}}(\mu)$.
Once again, using the continuity of the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$, we can write

$$ $\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\hat{\mu})<f_{\tau_{x,y}}^{-1}( \lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha+\epsilon_{n}}(\mu)}\right)\leq\mathbb{P }^{\mathcal{T}}\left({Q_{1-\alpha-\epsilon_{n}}(\mu)<f_{\bar{\tau}_{n+1}}^{-1} (\bar{\lambda}_{n+1})\leq Q_{\bar{\gamma}}(\mu)}\right)\\ =F\left({Q_{\bar{\gamma}}(\mu)}\right)-F\left({Q_{1-\alpha-\epsilon_{n}}(\mu)} \right)=\bar{\gamma}+\epsilon_{n}-1+\alpha\\ =n^{-1}+\sqrt{8\alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}}(\textup{ P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Plugging the previous inequality inside (23) yields

$$ $\mathbb{E}^{\mathcal{T}}\big{[}{p_{n+1}^{(X,Z)}}\big{]}\leq n^{-1}+\sqrt{8 \alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y}; \textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Therefore, (22) implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in \mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)<1-\alpha-\rho}\right )\\ \leq\frac{n^{-1}+\sqrt{8\alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}} (\textup{P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X})+\mathbb{E}^{\mathcal{T}} \big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X};\Pi_{Y\mid X})}\big{]}}{ \rho}.$ (24) $$

##### Step 2: Upper bound.

Using (21), we obtain

$$ $\mathbb{P}^{\mathcal{T}}\left({\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in \mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)>1-\alpha+\rho}\right )\\ \leq\mathbb{P}^{\mathcal{T}}\left({\Pi_{Y\mid X=X}(\mathcal{R}_{Z}(X;f_{\tau_{ X,Z}}(\varphi)))+\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X=X};\Pi_{Y\mid X= X})+q_{n+1}^{(X,Z)}>1-\alpha+\rho}\right).$ $$

The continuity of the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$ implies

$$ $1-\alpha=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y} ,Z})\leq\varphi}\right)=\int\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}( \varphi)))\bar{\Pi}_{Z\mid X=x}(\mathrm{d}z)\textup{P}_{X}(\mathrm{d}x).$ $$

Since $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))\geq 1-\alpha$, we deduce that $\Pi_{Y\mid X=x}(\mathcal{R}_{z}(x;f_{\tau_{x,z}}(\varphi)))=1-\alpha$ almost surely.
Therefore, using the Markov’s inequality gives

$$ $\mathbb{P}^{\mathcal{T}}\left({\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in \mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)>1-\alpha+\rho}\right )\leq\frac{\mathbb{E}^{\mathcal{T}}\big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{P} _{Y\mid X};\Pi_{Y\mid X})}\big{]}+\mathbb{E}^{\mathcal{T}}\big{[}{q_{n+1}^{(X, Z)}}\big{]}}{\rho}.$ (25) $$

Moreover, applying Theorem A.7 shows that

$$ $q_{n+1}^{(x,z)}\leq n^{-1}+\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\hat{ \mu})<f_{\tau_{x,y}}^{-1}(\lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha+\epsilon_{n}} (\mu)\,|\,X_{n+1}=x,Z_{n+1}=z}\right).$ (26) $$

Let’s define $\underline{\gamma}$ by

$$ $\underline{\gamma}=\min(1,1-\alpha-\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y}; \textup{P}_{X}\times\Pi_{Y\mid X})).$ $$

We now show that $Q_{\underline{\gamma}}(\mu)\leq Q_{1-\alpha}(\hat{\mu})$.
By continuity of the cumulative density function of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$, we have

$$ $\displaystyle 1-\alpha$ $\displaystyle=\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X, \hat{Y},Z})\leq Q_{1-\alpha}(\hat{\mu})}\right)$ $\displaystyle\leq\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z}}^{-1}(\lambda_{X ,Y,Z})\leq Q_{1-\alpha}(\hat{\mu})}\right)+\mathrm{d}_{\mathrm{TV}}(\textup{P} _{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Hence, it follows that

$$ $\displaystyle\underline{\gamma}\leq\mathbb{P}^{\mathcal{T}}\left({f_{\tau_{X,Z }}^{-1}(\lambda_{X,Y,Z})\leq Q_{1-\alpha}(\hat{\mu})}\right).$ $$

Thus, we deduce that $Q_{1-\alpha}(\hat{\mu})\geq Q_{\underline{\gamma}}(\mu)$.
Using the continuity of the distribution of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$, we can write

$$ $\mathbb{P}^{\mathcal{T}}\left({Q_{1-\alpha}(\hat{\mu})<f_{\tau_{x,y}}^{-1}( \lambda_{x,Y_{n+1},z})\leq Q_{1-\alpha+\epsilon_{n}}(\mu)}\right)\leq\mathbb{P }^{\mathcal{T}}\left({Q_{\underline{\gamma}}(\mu)<f_{\tau_{x,y}}^{-1}(\lambda_ {x,Y_{n+1},z})\leq Q_{1-\alpha+\epsilon_{n}}(\mu)}\right)\\ =F\left({Q_{1-\alpha+\epsilon_{n}}(\mu)}\right)-F\left({Q_{\underline{\gamma}} (\mu)}\right)=\epsilon_{n}-1+\alpha-\underline{\gamma}\\ =n^{-1}+\sqrt{8\alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}}(\textup{ P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Plugging the previous inequality inside (26) yields

$$ $\mathbb{E}^{\mathcal{T}}\big{[}{q_{n+1}^{(X,Z)}}\big{]}\leq n^{-1}+\sqrt{8 \alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y}; \textup{P}_{X}\times\Pi_{Y\mid X}).$ $$

Therefore, (25) implies that

$$ $\mathbb{P}^{\mathcal{T}}\left({\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in \mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)<1-\alpha-\rho}\right )\\ \leq\frac{n^{-1}+\sqrt{8\alpha(1-\alpha)n^{-1}\log n}+\mathrm{d}_{\mathrm{TV}} (\textup{P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X})+\mathbb{E}^{\mathcal{T}} \big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X};\Pi_{Y\mid X})}\big{]}}{ \rho}.$ (27) $$

##### Step 3: Bound on 𝔼 𝒯 ⁢ [ d TV ⁢ ( P Y ∣ X ; Π Y ∣ X ) ] superscript 𝔼 𝒯 delimited-[] subscript d TV subscript P conditional 𝑌 𝑋 subscript Π conditional 𝑌 𝑋 \mathbb{E}^{\mathcal{T}}\big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X};%
\Pi_{Y\mid X})}\big{]} blackboard_E start_POSTSUPERSCRIPT caligraphic_T end_POSTSUPERSCRIPT [ roman_d start_POSTSUBSCRIPT roman_TV end_POSTSUBSCRIPT ( P start_POSTSUBSCRIPT italic_Y ∣ italic_X end_POSTSUBSCRIPT ; roman_Π start_POSTSUBSCRIPT italic_Y ∣ italic_X end_POSTSUBSCRIPT ) ] .

Let’s denote $\nu_{Y\mid X=x}=2^{-1}(\textup{P}_{Y\mid X=x}+\Pi_{Y\mid X=x})$.
Since $\textup{P}_{Y\mid X=x}\ll\nu_{Y\mid X=x}$ and $\Pi_{Y\mid X=x}\ll\nu_{Y\mid X=x}$, there exists two Radon–Nikodym derivatives $g_{1}(x,\cdot)$ and $g_{1}(x,\cdot)$ of $\textup{P}_{Y\mid X=x}$ and $\Pi_{Y\mid X=x}$ with respect to $\nu_{Y\mid X=x}$.
Moreover, $g_{1}$ and $g_{2}$ are also the Radon–Nikodym derivatives of $\textup{P}_{X,Y}$ and $\textup{P}_{X}\times\Pi_{Y\mid X}$ with respect to $\textup{P}_{X}\times\nu_{Y\mid X}$.
By definition of the total variation distance, we have

$$ $\displaystyle\mathbb{E}^{\mathcal{T}}\big{[}{\mathrm{d}_{\mathrm{TV}}(\textup{ P}_{Y\mid X};\Pi_{Y\mid X})}\big{]}$ $\displaystyle=\int\mathrm{d}_{\mathrm{TV}}(\textup{P}_{Y\mid X};\Pi_{Y\mid X}) \textup{P}_{X}(\mathrm{d}x)$ $\displaystyle=\frac{1}{2}\int\left\lvert{g_{1}(x,y)-g_{2}(x,y)}\right\rvert\nu _{Y\mid X=x}\textup{P}_{X}(\mathrm{d}x)$ $\displaystyle=\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y};\textup{P}_{X}\times \Pi_{Y\mid X}).$ (28) $$

##### Step 4: Combination.

Finally, using (24)-(27) and (28), it follows that

$$ $\mathbb{P}^{\mathcal{T}}\left({\left\lvert{\mathbb{P}^{\mathcal{T}}\left({Y_{n +1}\in\mathcal{C}_{\alpha}(X_{n+1})\,|\,X_{n+1},Z_{n+1}}\right)-1+\alpha} \right\rvert>\rho}\right)\\ \leq\frac{2n^{-1}+\sqrt{128\alpha(1-\alpha)n^{-1}\log n}+4\mathrm{d}_{\mathrm{ TV}}(\textup{P}_{X,Y};\textup{P}_{X}\times\Pi_{Y\mid X})}{\rho}.$ $$

Note that the proof assumes $n^{-1}\log n\leq 8^{-3}\alpha(1-\alpha)$. To ensure the validity of the previous bound even when this assumption does not hold, we increased the term $\sqrt{32\alpha(1-\alpha)n^{-1}\log n}$ to $\sqrt{128\alpha(1-\alpha)n^{-1}\log n}$.
∎

###### Theorem A.9 .

Assume 1-2-3 hold.
If the distributions of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ and $f_{\tau_{X,Z}}^{-1}(\lambda_{X,\hat{Y},Z})$ are continuous, then, $\forall\epsilon\in(0,1)$ there exists $(\Lambda_{n}^{(\epsilon)})_{n\in\mathbb{N}}$ such that $\liminf_{n\to\infty}\mathbb{P}((X_{n+1},Z_{n+1})\in\Lambda_{n}^{(\epsilon)})
\geq 1-\epsilon$ and also

$$ $\sup_{(x,z)\in\Lambda_{n}^{(\epsilon)}}\left\lvert{\mathbb{P}^{\mathcal{T}} \left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n+1})\,|\,(X_{n+1},Z_{n+1})=(x,z)} \right)-1+\alpha}\right\rvert=\operatorname{\mathrm{O}}_{\mathbb{P}}\left({ \sqrt{n^{-1}\log n}+r_{n}}\right).$ $$

###### Proof.

First of all, define the following variables

$$ $\displaystyle c_{n+1}(x,z)=\left\lvert{\mathbb{P}^{\mathcal{T}}\left({Y_{n+1} \in\mathcal{C}_{\alpha}(X_{n+1})\,|\,(X_{n+1},Z_{n+1})=(x,z)}\right)-1+\alpha} \right\rvert,$ $\displaystyle d_{n}=\mathrm{d}_{\mathrm{TV}}(\textup{P}_{X,Y};\textup{P}_{X} \times\Pi_{Y\mid X}^{(m_{n})}).$ $$

Applying Theorem A.8, we obtain

$$ $\mathbb{P}\left({c_{n+1}(X_{n+1},Z_{n+1})>\rho}\right)\leq\mathbb{P}\left({d_{ n}>r_{n}}\right)+\mathbb{P}\left({c_{n+1}(X_{n+1},Z_{n+1})>\rho;d_{n}\leq r_{n }}\right)\\ \leq\mathbb{P}\left({d_{n}>r_{n}}\right)+\mathbb{E}\left[{\mathds{1}_{d_{n} \leq r_{n}}\mathbb{P}^{\mathcal{T}}\left({c_{n+1}(X_{n+1},Z_{n+1})>\rho}\right )}\right]\\ \leq\mathbb{P}\left({d_{n}>r_{n}}\right)+\frac{2n^{-1}+\sqrt{128\alpha(1- \alpha)n^{-1}\log n}+4r_{n}}{\rho}.$ $$

Finally, using 3, we get $\lim_{n\to\infty}\mathbb{P}(d_{n}>r_{n})=0$.
Therefore, for any $\epsilon>0$, there exist $M_{\epsilon}>0$ and $\tilde{n}_{\epsilon}\in\mathbb{N}$ such that, $\forall n\geq\tilde{n}_{\epsilon}$, it holds

$$ $\mathbb{P}\left({c_{n+1}(X_{n+1},Z_{n+1})>M_{\epsilon}\times\left({\sqrt{n^{-1 }\log n}+r_{n}}\right)}\right)\leq\epsilon.$ (29) $$

Given $\epsilon\in(0,1)$, let’s consider the following set

$$ $\Lambda_{n}^{(\epsilon)}=\left\{{(X_{n+1}(\omega),Z_{n+1}(\omega))\colon\omega \in\Omega,c_{n+1}(X_{n+1},Z_{n+1})(\omega)\leq M_{\epsilon}\times\left({\sqrt{ n^{-1}\log n}+r_{n}}\right)}\right\}.$ $$

Equation 29 implies that

$$ $\liminf_{n\to\infty}\mathbb{P}\left({(X_{n+1},Z_{n+1})\in\Lambda_{n}^{( \epsilon)}}\right)\geq 1-\epsilon,$ $$

and by definition of $\Lambda_{n}^{(\epsilon)}$, we also have

$$ $\sup_{(x,z)\in\Lambda_{n}^{(\epsilon)}}c_{n+1}(x,z)=\operatorname{\mathrm{O}}_ {\mathbb{P}}\left({\sqrt{n^{-1}\log n}+r_{n}}\right).$ $$

∎

Note that, (29) also shows that

$$ $\left\lvert{\mathbb{P}^{\mathcal{T}}\left({Y_{n+1}\in\mathcal{C}_{\alpha}(X_{n +1})\,|\,X_{n+1},Z_{n+1}}\right)-1+\alpha}\right\rvert=\operatorname{\mathrm{O }}_{\mathbb{P}}\left({n^{-1/2}\sqrt{\log n}+r_{n}}\right).$ $$

###### Theorem A.9 .

###### Proof.

### A.4 Additional results

Let’s denote the conditional c.d.f of $f_{\tau_{X,Z}}^{-1}(\lambda_{X,Y,Z})$ by

$$ $F_{x,z}(\cdot)=\int_{\mathbb{R}^{d}\times\mathcal{Z}}\mathbb{P}\left({f_{\tau_ {x,z}}^{-1}(\lambda_{x,Y,z})\leq\cdot\mid(X,Z)=(x,z)}\right)\bar{\Pi}_{Z\mid X =x}(\mathrm{d}z)\textup{P}_{X}(\mathrm{d}x).$ $$

###### Lemma A.10 .

Assume that $Q_{1-\alpha}(\mu_{n})\to\varphi$ almost-surely as $n\to\infty$.
If $F_{X,Z}$ is continuous almost-surely, then $\lim_{n\to\infty}p_{n+1}^{(x,z)}=0$, $\bar{\Pi}_{Z\mid X}\times\textup{P}_{X}$-almost everywhere.

###### Proof.

First, define the following sets:

$$ $\displaystyle A=\left\{{\omega\in\Omega\colon\lim_{n\to\infty}Q_{1-\alpha}(\mu _{n}(\omega))=\varphi}\right\},$ $\displaystyle B=\left\{{\omega\in\Omega\colon F_{X(\omega),Z(\omega)}\text{ is continuous}}\right\}.$ $$

For all $\omega\in A\cap B$, it holds

$$ $\lim_{n\to\infty}F_{X(\omega),Z(\omega)}\left({Q_{1-\alpha}\left({\mu_{n}( \omega)}\right)\wedge\varphi}\right)=F_{X(\omega),Z(\omega)}\left({\varphi} \right).$ $$

Moreover, note that we can write

$$ $p_{n+1}^{(x,z)}=F_{x,z}(\varphi)-F_{x,z}(\varphi\wedge Q_{1-\alpha}(\mu_{n})).$ $$

Hence, we deduce that

$$ $\displaystyle 1=\mathbb{P}\left({A\cap B}\right)$ $\displaystyle\leq\mathbb{P}\left({\omega\in\Omega\colon\lim_{n\to\infty}F_{X( \omega),Z(\omega)}\left({Q_{1-\alpha}\left({\mu_{n}(\omega)}\right)\wedge \varphi}\right)=F_{X(\omega),Z(\omega)}\left({\varphi}\right)}\right)$ $\displaystyle=\mathbb{P}\left({\lim_{n\to\infty}p_{n+1}^{(X,Z)}=0}\right)$ $\displaystyle=\int_{\mathbb{R}^{d}\times\mathcal{Z}}\mathbb{P}\left({\lim_{n \to\infty}p_{n+1}^{(x,z)}=0\,\big{|}\,(X,Z)=(x,z)}\right)\,\textup{P}_{Z\mid X =x}(\mathrm{d}z)\,\textup{P}_{X}(\mathrm{d}x).$ $$

The last line implies that $p_{n+1}^{(x,z)}\to 0$ almost $\textup{P}_{Z\mid X}\times\textup{P}_{X}$-everywhere.
∎

The prediction set, defined in (4), is derived from the $(1-\alpha)$-quantile of the conformity scores $\{f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})\}_{k=1}^{n}\cup\{\infty\}$.
However, $\{\infty\}$ can be removed from these conformity scores.
Inspired by Romano et al. (); Sesia & Candès (), we prove a corollary of Theorem 3.1. Its result demonstrates the marginal validity of the prediction set defined as

$$ $\bar{\mathcal{C}}_{\alpha}(x)=\mathcal{R}_{z}\left({x;f_{\tau_{x,z}}\Big{(}{Q_ {(1-\alpha)(1+n^{-1})}\big{(}{\textstyle\frac{1}{n}\sum_{k=1}^{n}\delta_{f_{ \bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})}}\big{)}}\Big{)}}\right).$ (30) $$

While the prediction set $\bar{\mathcal{C}}_{\alpha}(x)$ relies on the quantile of the distribution $\frac{1}{n}\sum_{k=1}^{n}\delta_{f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})}$, its proof reveals that this prediction set is equivalent to $\mathcal{C}_{\alpha}(x)$.

###### Corollary A.11 .

Under the same assumptions as in Theorem 3.1, for any $\alpha\in[1/(n+1),1]$, we have

$$ $1-\alpha\leq\mathbb{P}\left({Y_{n+1}\in\bar{\mathcal{C}}_{\alpha}(X_{n+1})} \right)<1-\alpha+\frac{1}{n+1},$ $$

where the upper bound only holds if the conformity scores $\{f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k})\}_{k=1}^{n+1}$ are almost surely distinct.

###### Proof.

Let $\alpha\in\mathbb{R}$ such that $(n+1)^{-1}\leq\alpha\leq 1$, and recall that

$$ $\mu_{n}=\frac{1}{n+1}\sum_{k=1}^{n}\delta_{f_{\bar{\tau}_{k}}^{-1}(\bar{ \lambda}_{k})}+\frac{1}{n+1}\delta_{\infty}.$ $$

Since $\alpha\geq(n+1)^{-1}$, the quantile $Q_{1-\alpha}(\mu_{n})$ is the $k_{\alpha}$th order statistic of $V_{1},\ldots,V_{n}$, where

$$ $V_{k}=f_{\bar{\tau}_{k}}^{-1}(\bar{\lambda}_{k}),\quad\text{and}\quad k_{ \alpha}=\lceil(1-\alpha)(n+1)\rceil.$ $$

However, $\forall\beta\in(\frac{k_{\alpha}-1}{n},\frac{k_{\alpha}}{n}]$, we have

$$ $Q_{\beta}\left({\textstyle\frac{1}{n}\sum_{k=1}^{n}\delta_{V_{k}}}\right)=V_{( k_{\alpha})}.$ $$

Since $\mathcal{C}_{\alpha}(X_{n+1})=\mathcal{R}_{Z_{n+1}}({X_{n+1};f_{\bar{\tau}_{n+
1}}({V_{(k_{\alpha})}})})$, Theorem 3.1 implies that

$$ $1-\alpha\leq\mathbb{P}\left({Y_{n+1}\in\mathcal{R}_{Z_{n+1}}\left({X_{n+1};f_{ \bar{\tau}_{n+1}}\left({Q_{\beta}\left({\textstyle\frac{1}{n}\sum_{k=1}^{n} \delta_{V_{k}}}\right)}\right)}\right)}\right)<1-\alpha+\frac{1}{n+1}.$ $$

Setting $\beta=(1-\alpha)(1+n^{-1})$ in the previous inequality and using the definition of $\bar{\mathcal{C}}_{\alpha}(X_{n+1})$ given in (30) concludes the proof.
∎

###### Lemma A.10 .

###### Proof.

###### Corollary A.11 .

###### Proof.

## Appendix B Experimental Setup and Results

This section aims to provide a comprehensive understanding of the CP^2 algorithm. We want to further explore the CP^2 approach and to better explain the key concepts.

##### Choice of f t subscript 𝑓 𝑡 f_{t} italic_f start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT .

We present examples of mappings $f_{t}$ and their inverses $f_{\tau}^{-1}$ in Table 2.
The choice of the mapping $f_{t}$ is crucial for the performance of the method, and we investigate their impact in Section 4.
For instance, choosing $f_{\tau}(\lambda)=\tau\lambda$ results in approximately conditionally valid prediction sets, as long as $\Pi_{Y\mid X=x}$ accurately estimates the conditional distribution $\textup{P}_{Y\mid X=x}$; see Theorem 3.2-Theorem 3.3. Initially we also considered other adjustment functions based on exponent, sigmoid and $\tanh$ functions, but they all performed worse than linear and sum. As we show later in 3, these two selected adjustment function perform similarly, showing only marginal differences on some datasets. Designing new adjustment functions is a possible future research direction.

### B.1 Highest predictive density (HPD) regions

##### CP 2 with Explicit Conditional Density estimate: CP 2 -HPD .

Assume that an estimator the conditional density function is known, denoted by $\gamma_{Y\mid X=x}$.
The confidence set is defined as $\mathcal{R}(x;t)=\{y\in\mathcal{Y}\colon\gamma_{Y\mid X=x}(y)\geq-t\}$. We omit the variable $z$ from the notation, as we do not consider exogenous randomization in this case.
The parameter $\tau_{x}$ is obtained by solving

$$ $\textstyle\tau_{x}=\operatornamewithlimits{\arg\min}\left\{\tau\in\mathbb{R} \colon\int_{\mathcal{R}(x;\tau)}\gamma_{Y\mid X=x}(y)\,\mathrm{d}y\geq 1- \alpha\right\}.$ (31) $$

We then compute $\lambda_{x,y}=-\gamma_{Y\mid X=x}(y)$ and derive the prediction set as

$$ $\textstyle\mathcal{C}_{\alpha}(x)=\left\{{y\in\mathcal{Y}\colon\gamma_{Y\mid X =x}(y)\geq-f_{\tau_{x}}\left({Q_{1-\alpha}\left({\mu_{n}}\right)}\right)} \right\}.$ $$

If we take $f_{\tau}(\lambda)=\lambda$ and $\varphi=1$, the method shares similarity with the CD-split method, proposed in Izbicki et al. ().
While CD-split uses $\lambda_{x,y}$ as the conformity score, our method uses $f_{\tau_{x}}^{-1}(\lambda_{x,y})$, which incorporates the information from $\tau_{x}$ to modify $\gamma_{Y\mid X=x}(y)$. The CP^2-HPD workflow is summarized in Algorithm 2.

Of course, the computation of (31) is in general highly non-trivial. Izbicki et al. () suggested to use binning, therefore approximating the conditional predictive distribution with histograms. The method is restricted to the case where the dimension of $\mathcal{Y}$ the response is small; see Izbicki et al. () for the case of $\mathcal{Y}=\mathbb{R}$. When the dimension becomes larger, then the estimation of HPD is typically based on Monte Carlo methods, thus requiring the introduction of auxiliary variables.

Figure: Algorithm 2 CP^2-HPD

**Table 2: Adjustment Functions $f_{t}$, their inverses $f_{\tau}^{-1}$ and $\varphi$ values used in our experiments.**
| Name | $f_{\tau}(\lambda)$ | $f_{\tau}^{-1}(\lambda)$ | $\varphi$ |
| --- | --- | --- | --- |
| Linear | $\tau\lambda$ | $\tau^{-1}\lambda$ | 1 |
| Difference | $\tau+\lambda$ | $\lambda-\tau$ | 0 |

### B.2 Details of the experimental setup

We use the Mixture Density Network (Bishop, ) implementation from CDE (Rothfuss et al., ) Python package(^2^22[https://github.com/freelunchtheorem/Conditional_Density_Estimation](https://github.com/freelunchtheorem/Conditional_Density_Estimation)) as a base model for CP, PCP and CP^2. The underlying neural network contains two hidden layers of 100 neurons each and was trained for 1000 epochs for each split of the data. Number of components of the Gaussian Mixture was set to 10 for all datasets.

For the CQR (Romano et al., ) and CHR (Sesia & Romano, ) we use the original authors’ implementation(^3^33[https://github.com/msesia/chr](https://github.com/msesia/chr)). The underlying neural network that outputs conditional quantiles consists of two hidden layers with 64 neurons each. Training was performed for 200 epochs for batch size 250.

We replicate the experiments for 50 random splits of all nine datasets. To lower noise in calculated performance metrics we reuse trained networks and samples across different top-level algorithms for each replication.

### B.3 Worst-slab coverage

Here we present some additional experiments related to conditional coverage achieved by different methods. We have used Worst Slab Coverage metric, which is sensitive to the set of labs considered during the search. Following Cauchois et al. (); Romano et al. (), recall that a slab is defined as

$$ $S_{v,a,b}=\left\{x\in\mathbb{R}^{p}:a<v^{T}x<b\right\},$ $$

where $v\in\mathbb{R}^{p}$ and $a,b\in\mathbb{R}$, such that $a<b$. Now, given the prediction set $\mathcal{C}(x)$ and $\delta\in[0,1]$, the worst-slab coverage is defined as:

$$ $\mathrm{WSC}(\mathcal{C},\delta)=\inf_{v\in\mathbb{R}^{p},a<b\in\mathbb{R}} \mathbb{P}\left(Y\in\mathcal{C}(X)|X\in S_{v,a,b}\right)\>s.t.\>\mathbb{P}(X \in S_{v,a,b})\geq 1-\delta.$ $$

In our experiments we follow Romano et al. () in our implementation of this metric. Namely, we use $25\%$ of the data to find the worst slab and the use the remaining $75\%$ to calculate the final value on this slab. We use $5000$ randomly sampled directions, that are the same for each algorithm and change for each replication.

### B.4 Extended results of real data experiments

Table 3 we summarize all metrics from our real-world data experiments. For conditional coverage we report worst-slab coverage with $(1-\delta)=0.1$. On six out of nine datasets CP^2 method achieves the best result in conditional coverage. In terms of interval width PCP method produces the narrowest intervals.As we can see, it happens at the expense of conditional coverage: PCP often achieves significantly lower values.

We also present a more detailed view of set size differences between the methods. In the main part we reported average rank of each method in Figure 8. We ranked the algorithms by their projected area at each test point and averaged the ranks. Here we show raw areas of the projections onto each pairs of axes for sgemm_small dataset in table 4. All targets were standardised to zero mean and unit standard deviation so that different projections will be in the same scale. We see that PCP produces smaller set sizes like in one-dimensional case. Quantile-regression based methods have the largest sets, even larger than the fixed-sized sets of CP. Our approach demonstrates only modest increase in prediction set size compared to PCP while achieving sharper conditional coverage.

**Table 3: Summary results of experiments on real data. “M. Cov.” stands for marginal coverage, “C. Cov.” is the worst-slab coverage (here $(1-\delta)=0.1$), and $w_{sd}$ is average total length of the prediction sets, scaled by standard deviation of $Y$. Nominal coverage level is set to $(1-\alpha)=0.9$. For $\Pi_{Y\mid X}$, PCP, CP^2-PCP we use the same underlying mixture density network model with $50$ samples. CHR and $\texttt{CQR}(2)$ also share the same base neural network model. We average results of $50$ random data splits. For each dataset, we highlighted the algorithm achieving conditional coverage closest to the nominal level.**
|  |  |  |  |  | $\texttt{CP}^{2}$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Dataset | Metric | CP | PCP | $\Pi_{Y|X}$ | PCP-L | PCP-D | CHR | CQR | CQR2 |
| bike | M. Cov. | 0.90 | 0.90 | 0.93 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.79 | 0.85 | 0.92 | 0.89 | 0.89 | 0.88 | 0.90 | 0.87 |  |
| $w_{sd}$ | 0.71 | 0.71 | 0.83 | 0.79 | 0.80 | 1.94 | 2.25 | 2.31 |  |
| bio | M. Cov. | 0.90 | 0.90 | 0.91 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.88 | 0.89 | 0.91 | 0.90 | 0.90 | 0.90 | 0.89 | 0.89 |  |
| $w_{sd}$ | 2.34 | 1.89 | 1.95 | 1.97 | 1.95 | 1.92 | 2.13 | 2.10 |  |
| blog | M. Cov. | 0.90 | 0.90 | 0.91 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.60 | 0.74 | 0.91 | 0.89 | 0.90 | 0.87 | 0.87 | 0.86 |  |
| $w_{sd}$ | 0.60 | 0.30 | 0.72 | 0.72 | 0.71 | 0.31 | 0.44 | 0.39 |  |
| fb1 | M. Cov. | 0.90 | 0.90 | 0.93 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.49 | 0.64 | 0.92 | 0.88 | 0.89 | 0.87 | 0.90 | 0.87 |  |
| $w_{sd}$ | 0.47 | 0.28 | 0.58 | 0.56 | 0.59 | 0.26 | 0.37 | 0.33 |  |
| fb2 | M. Cov. | 0.90 | 0.90 | 0.93 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.50 | 0.61 | 0.91 | 0.88 | 0.88 | 0.88 | 0.89 | 0.89 |  |
| $w_{sd}$ | 0.53 | 0.32 | 0.65 | 0.62 | 0.65 | 0.33 | 0.43 | 0.37 |  |
| meps19 | M. Cov. | 0.90 | 0.90 | 0.89 | 0.90 | 0.90 | 0.90 | 0.89 | 0.90 |
| C. Cov. | 0.54 | 0.78 | 0.89 | 0.90 | 0.89 | 0.90 | 0.88 | 0.89 |  |
| $w_{sd}$ | 1.05 | 0.73 | 1.02 | 1.19 | 1.07 | 0.76 | 1.14 | 1.19 |  |
| meps20 | M. Cov. | 0.90 | 0.90 | 0.89 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.58 | 0.80 | 0.89 | 0.90 | 0.90 | 0.91 | 0.88 | 0.89 |  |
| $w_{sd}$ | 1.06 | 0.75 | 0.98 | 1.15 | 1.04 | 0.77 | 1.09 | 1.17 |  |
| meps21 | M. Cov. | 0.90 | 0.90 | 0.89 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.54 | 0.81 | 0.89 | 0.89 | 0.89 | 0.90 | 0.89 | 0.88 |  |
| $w_{sd}$ | 1.04 | 0.72 | 0.99 | 1.16 | 1.04 | 0.79 | 1.13 | 1.21 |  |
| temp | M. Cov. | 0.90 | 0.90 | 0.82 | 0.90 | 0.90 | 0.90 | 0.90 | 0.90 |
| C. Cov. | 0.87 | 0.89 | 0.81 | 0.88 | 0.89 | 0.86 | 0.85 | 0.86 |  |
| $w_{sd}$ | 0.87 | 0.92 | 0.78 | 0.96 | 0.93 | 1.31 | 1.48 | 1.30 |  |

**Table 4: Prediction set size comparison for sgemm_small dataset. Rows correspond to different pairs of targets (dataset has 4 targets). For each method the reported value is the mean area of the 2D projection of the prediction set to the corresponding axes pair.**
|  |  |  |  | $\texttt{CP}^{2}$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Axes | CP | PCP | $\Pi_{Y|X}$ | PCP-L | PCP-D | CHR | CQR | CQR2 |
| (0, 1) | 2.137 | 0.435 | 0.517 | 0.576 | 0.560 | 2.290 | 2.550 | 2.436 |
| (0, 2) | 2.145 | 0.435 | 0.518 | 0.577 | 0.561 | 2.267 | 2.506 | 2.358 |
| (0, 3) | 2.145 | 0.436 | 0.519 | 0.578 | 0.561 | 2.086 | 2.366 | 2.172 |
| (1, 2) | 2.146 | 0.435 | 0.517 | 0.576 | 0.560 | 2.388 | 2.622 | 2.546 |
| (1, 3) | 2.146 | 0.435 | 0.517 | 0.576 | 0.560 | 2.166 | 2.461 | 2.314 |
| (2, 3) | 2.154 | 0.436 | 0.519 | 0.578 | 0.562 | 2.153 | 2.430 | 2.255 |

### B.5 Other perspective on conditional coverage

The worst-slab coverage metric used in the previous section is not always helpful: (1) it provides a single number for each method, and (2) the selected slab is different for each algorithm. In practice we might be interested in how sharp the coverage is along the portion of the input space spanned by the test data. To explore this, we used two approaches: dimensionality reduction and clustering. Results for clustering with HDBSCAN are presented in the main part in Figure 5, here turn to dimensionality reduction.

First we apply UMAP algorithm to project data to two dimensions and then construct a heatmap plot to show coverage in each bin of the histogram. Results for meps_19 dataset are presented in Figure 9. Nominal coverage is set to $(1-\alpha)=0.9$ and corresponds to gray part of the color scale. We can see that our method and baseline $\Pi_{Y\mid X}$ perform better than CP and PCP across the space.

Figure: Figure 9: Conditional coverage after dimensionality reduction, meps_21 dataset. Data projected to two dimensions using UMAP algorithm with Canberra metric, with the n_neighbors hyperparameter set to 2. Nominal coverage is set to $(1-\alpha)=0.1$, it corresponds to gray on the color scale.
Refer to caption: x11.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: No further questions. I might ask more once I got the authors response.
soundness: 3
strengths: - The issue of conditional coverage of conformal prediction methods is of particular interest in CP community. They expand the ideas of using an estimation of conditional coverage to enhance the covariate-conditional performance of CP to the more general setting of multi-dimensional prediction tasks.
- They also provide theoretical guarantees on the relation between the quality of the conditional distribution estimation and the improvements in the conditional coverage property of CP, which was a missing flavor for most of the previous works in this line of thinking.
summary: This paper looks at the central issue of covariate-conditional validity of conformal prediction methods. They propose a systematic and novel procedure of turning a black box estimation of conditional distributions (Y|X) into a valid conformity score which then can be used to construct prediction sets. They carefully design the procedure to make their prediction sets adaptive to the behavior of the conditional distribution. Their method effectively generalizes the previous literature on using conditional distribution estimation by 1) providing algorithmic and theoretical results targeting conditional coverage and 2) going beyond one-dimensional regression setting.
weaknesses: - In general, the whole line of work based on estimations of conditional density suffers from the curse of dimensionality (covariate dimension). Therefore, it is important to see some empirical evaluations of the proposed method on high dimensional tasks (in particular image classification), even if it doesn't work very well, it helps the follow-up works with having a more accurate understanding of your method. Another interesting experiment to see is a synthetic regression task that you can directly plot the (conditional) coverage properties vs. dimension of the input and showcase the degree of tolerance of your proposed method. 
- Focusing on the lower dimensional tasks, where conditional density estimation is more feasible, It is good to see the performance of your proposed method in comparison with recent methods like the one proposed by https://arxiv.org/abs/2305.12616. You can for instance set their function class to be a linear map form the data itself to a scalar threshold, or train a neural network regression model in the training time and then set the function class to be a linear head from the last layer of the trained neural network to a scalar threshold. In general, it is important to see whether using the conditional density estimation (when feasible) can actually beat the methods that does not rely on conditional density estimation.
- Taking a step back, intuitively, it is obvious that the usage of an estimate of the conditional density estimation should improve the conditional coverage behavior of CP. However, it is not entirely obvious whether all the information in the density function is needed to construct perfect prediction sets. In other words, the key question is what parts of the information in the conditional density function are actually needed for constructing prediction sets. This question is very important to look at as the answer to it might reveal a way to beat the curse of dimensionality, as then you can only try to learn the useful information instead of the whole conditional density function. To this end, 1) it might be good (as a suggestion) to discuss this point of view (maybe as a future work) in your paper. In particular, how fundamental conditional density estimation is for your pipeline. 2) It is good to discuss the methods that are based on the described idea in your paper, for instance https://arxiv.org/abs/2404.17487, as that helps the reader to have a better understanding of the whole methodology.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: It would be nice if the authors added in some analysis of how likely or under what conditions these assumptions tend to hold, or rather how restrictive they are. It seems they are not that restrictive, but some analysis here would be nice.
soundness: 4
strengths: 1. the underlying approach of estimating the conditional distribution via sampling and proposing the conformal set using these samples (to my knowledge) is unique to this paper
2. This paper has very nice asymptotic guarantees of performance under different underly data distributions
3. The approach is very simple and applicable, making it generally useful in many fields. 
4. The empirical validation is very impressive
summary: This paper achieves powerful conditional coverage via a technique that relies on sampling from an empirical estimate of the underlying conditional data distribution. Their method achieves impressive empirical benefits alongside powerful provable guarantees. 

Overall this paper has clear contributions and is very sound. I vote to accept this paper.
weaknesses: 1. Some of the figures are not accessible. It is difficult to see for someone who is red-green colorblind. This includes Figure 1. 
2. Some of the presentation of the theoretical ideas, especially on page 6, is difficult to parse. Better formatting would increase readability.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. About the theoretical guarantees, there is a lack of rigorous theoretical discussion on the superiority of the size (length) of intervals. Is it possible to provide any properties about the superity in terms of size? Including some theoretical investigations would be appreciated.

2. Regarding the experiments, I notice that the authors didn't compare with some important baselines about this problem, for example, the LCP in [1], RLCP in [2] and conditional calibration method in [3] should be considered. I'm curious about the performance of the CP^2 method in comparison to these works in this area, and additional experiments would be beneficial to provide further insights.

References:

[1]  Guan, L. (2023). Localized conformal prediction: A generalized inference framework for conformal prediction. Biometrika, 110(1):33–50.  

[2]  Hore, R. and Barber, R. F. (2023). Conformal prediction with local weights: randomization enables local guarantees. arXiv preprint arXiv:2310.07850.  

[3]  Gibbs, I., Cherian, J. J., and Candes, E. J. (2023). Conformal prediction with conditional guarantees. arXiv preprint arXiv:2305.12616.
soundness: 3
strengths: This paper provides a specific approach to leverage information from the conditional distribution to improve conditional coverage. Besides estimating conditional distributions directly, the authors also consider alternatives using a conditional generative model when identifying HPDs is challenging.
summary: The authors propose a $\text{CP}^2$ framework for conformal prediction that aims to achieve approximately conditional coverage. The key difference from classical split CP is the introduction of a parameter, $\tau$, defined through either the estimation of the conditional distribution $P_{Y|X}$ or the use of a conditional generative model, allowing for adaptive threshold adjustments. Theoretically, they provide a finite-sample marginal coverage guarantee as well as an asymptotic guarantee through a upper bound on conditional coverage. The framework's efficiency is demonstrated through experiments on both synthetic and real-world data.
weaknesses: 1. The paper's writing could be improved, particularly in terms of clarity and presentation. For example, the literature review lacks clarity, and could benefit from a more comprehensive and topic-oriented discussion of related work.

2. Regarding the experimental section, it is noted that several key baselines from existing literature have been overlooked. Additionally, it seems that the superiority of CP² in terms of worst-slab coverage is not particularly evident comparing with CQR/CQR2 in Figure 3.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: It would be great to provide some intuition on when PCP-D is preferred over PCP-L (and vice versa).
soundness: 3
strengths: 1. A conformal framework to construct conformal prediction sets with improved adaptivity.

2. Thorough theoretical study (marginal, conditional, asymptotic analysis for the coverage).

3. Thorough experiments, comparing CP^2 to leading methods except the multi-dim. case; see comment below.
summary: The paper offers a new score function for constructing conformal prediction sets to improve adaptivity and approximated conditional coverage.

Instead of explicitly relying on a score function, they use a confidence set function---an idea that is known in the literature. The proposed construction of the prediction set suggests scaling the conformal scores \lambda_k by \tau_k, extracting the 1-alpha quantile Q, and then applying the "inverse scaling" of Q by using f_{\tau}.  

Since their method requires estimating Y|X which is generally hard---especially when Y is high dim.---they offer an approach to using a sampling generative model for constructing the prediction sets (Algorithm 1); inspired by Wang et al.

The authors provide a marginal guarantee, approximated conditional coverage (depending on the quality of the estimator), and an asymptotic conditional coverage guarantee.
weaknesses: 1. While the paper offers a general construction for prediction sets, it shares ideas with existing methods such as (1) local conformal---the scaling part; (2) HPD split or distributional conformal prediction by Chernozhukov et al.---both using the conditional distribution estimator; and (3) sampling using generative models, inspired by Wang et al. Of course, the paper nicely combines these ideas into one framework, which is great, but I am not sure it brings a novel concept. I would be happy to get the author's feedback here: what is the concrete edge of the proposed method compared to prior work that also builds on conditional distribution estimation to promote adaptive/conditional coverage?

2. I suggest better explaining the role of Z. It was not clear to me what the variable Z is until reaching line 252. For example, in line 152, it is not clear at all what you mean by "sampled according to a kernel Z|X". Consider improving the exposition in line 152, e.g., by providing a concrete example.

3. Section 4.1: how does the CP^2-HPD perform compared to HPD-split by Izbicki et al. (2022)? Such a comparison is missing, and it is important to include it to understand the advantages of the proposed method. Please clarify why you expect the proposed method to improve performance.

4. The experiments with multi-dimensional response lack stronger baseline methods, such as the one by Feldman et al (reference below) and/or related techniques. If a full experimental comparison is not feasible, it is crucial to discuss the conceptual differences between the methods.

S. Feldman, S. Bates, and Y. Romano, "Calibrated Multiple-Output Quantile Regression with Representation Learning," JMLR, 2022.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. Since the proposed method is built on many existing works, could the authors explain more concretely their main contribution about the methodology (like in which step)? This would be helpful for the readers to get the main idea that do not exist in previous literature.

2. As I have stated in the Strengths part, I believe the method will outperform most conditional conformal methods in multi-modal settings, which is shown in the synthetic experiment.  However, I'm still curious about the comparison result under the common unimodal case. It would be better to include some related experiments.

3. I have noticed that an important reference [1] is missing. The LCP method [2] is cited but has not been considered in experiments. As important progress in conditional conformal prediction, these two methods should act as good baselines in the experiment. I will be willing to raise my score with either more experiments or more theoretical results showing superiority to the above two methods. 


[1]Gibbs, I., Cherian, J., and Cand`es, E. (2024), “Conformal Prediction With Conditional Guarantees,” arXiv preprint arXiv:2305.12616.

[2]Guan, L. (2023), “Localized conformal prediction: A generalized inference framework for conformal prediction,” Biometrika, 110, 33–50.
soundness: 3
strengths: The nonconformity score function proposed is similar to dividing the residual score by a conditional standard deviation estimator, which is replaced by a factor relating the conditional distribution of Y|X by the authors. In my opinion, the proposed method is especially suitable for multi-modal settings, in which directly characterizing the overall upper and lower prediction bound is not effective.
summary: This paper proposes a conformal prediction method to enhance conditional coverage while still ensuring a finite-sample marginal coverage guarantee. The authors consider estimating the conditional distribution of $Y|X$ on the training data and use this estimator to define a localization parameter $\tau$ for each calibration and test sample. These parameters are then used to adjust non-conformity score values. The authors prove the marginal validity of their method and give an inflation bound for the conditional coverage. They also prove the asymptotic conditional validity of their proposed interval in the oracle sense. Experimental results show slight superiority of the proposed method compared with other conditional conformal prediction methods.
weaknesses: Despite reasonable, the proposed method only shows slight empirical improvement compared with existing methods. The theoretical results are also achievable by methods like the LCP (Guan, 2023). There is a lack of theoretical results regarding properties like optimality or superiority, etc.

Meta-review: The authors leverage an estimated conditional distribution $P_{Y|X}$ to construct prediction sets that adapt to the local structure of the data, achieving approximate conditional validity. They define confidence sets $R_z(x; t)$, calibrated using a parameter $\tau_{x,z}$ derived from $P_{Y|X}$, and use conformity scores $\lambda_{x,y,z}$ to ensure these sets dynamically adjust to heteroscedasticity and multimodality. For complex distributions, the method supports unions of high-density regions, enabling effective prediction even in multimodal cases. The framework guarantees approximate conditional coverage, with errors bounded by the total variation distance between the true and estimated conditional distributions, making it robust and adaptable to various machine learning tasks.

As a core technical novelty, the conformity score $\lambda_{x,y,z}$ is defined as the minimum "effort" i.e. size $t$ of the confidence set $R_z(x; t)$ needed to include the observed $y$. Intuitively, This score captures the relationship between the observation $y$ and the local structure of $P_{Y|X}$, allowing the prediction sets to adapt not only to the estimated distribution but also to the specific data point being predicted. 

The method is innovative but heavily dependent on high-quality conditional distribution estimators. In a strict sense, the proposed method does not provide fundamentally stronger conditional guarantees than existing methods. Asymptotically, conditional coverage has been established in numerous previous conformal prediction paper as well as consistency analysis for estimation of conditional distribution. So it is pretty difficult to grasp what new problem is being solved in this paper and what new it brings to the current literature. 

The paper suffers from significant clarity issues, primarily due to its dense and cluttered presentation. It overloads the reader with mathematical rigor, complex notations, and auxiliary variables without first providing an intuitive understanding of the method, making it difficult to grasp the core ideas. And at the end, we can hardly extract a significant novelty. This is really concerning and I encourage the authors to not obfuscate the readers. All the marginal coverage established in this paper hold directly from well known coverage results in conformal prediction but are presented as a *"new result"*.

---
*While I acknowledge that the paper would benefit greatly from simplifying the exposition and presenting the contributions in a more transparent way, the reviewers' post-rebuttal feedback indicates that the contributions are sound and sufficient for acceptance. The program chairs can weigh these merits alongside the remaining concerns to make the final decision.*

Decision: Accept (Poster)