Conformal Prediction under Lévy-Prokhorov Distribution Shifts: Robustness to Local and Global Perturbations

We are very grateful for your time and effort in reviewing our paper. Thank you for the positive assessment of our work. Below, we first provide answers to your questions and then to the highlighted weaknesses.

---

Question 1. We proceed to answer to the three weaknesses.

Weakness 1. We appreciate the suggestion to elaborate on the connection with [9]. Our work was indeed inspired by their formulation (as noted in the introduction) and, in particular, by their treatment of the fact that only samples from the calibration distribution are available. This influence is reflected in Theorem 4.1, whose proof adapts ideas from their Theorem 1 to our setting.

That said, there are key differences in ideas, methodology, theoretical formulation, and proof techniques. In terms of ideas, we replace the $f$‑divergence ambiguity sets in [9] with Lévy–Prokhorov (LP) sets. $f$‑divergences cannot accommodate distribution shifts with disjoint support, which occur in practice under adversarial perturbations, label shifts, or severe data contamination. The LP framework addresses this limitation and offers a natural decomposition into local and global perturbations, enabling separate reasoning about small‑magnitude shifts and large‑magnitude deviations.

Methodologically, our shift model is defined via an optimal transport (OT) discrepancy, placing it in a different mathematical class than $f$‑divergences. OT‑based sets have a distinct geometry and constraint structure, requiring tools from OT theory rather than the convex‑analytic duality used for $f$‑divergences. This also shapes the proof strategy: in Propositions 3.4 and 3.5, we explicitly construct worst‑case distributions within LP sets and derive closed‑form characterizations of worst‑case quantiles and coverage guarantees. These arguments, tailored to the LP geometry, differ fundamentally from the duality‑based derivations in [9] and, we believe, are more broadly applicable to other OT‑driven robustness models.

Weakness 2. Thank you for the suggestion. Space permitting, we will try to move as much of the discussion on the estimation of $(\varepsilon,\rho)$ as possible to the main text.

Weakness 3. We agree that experiments on regression settings would further strengthen the paper. In response to your suggestion, we are adding experiments on the PovertyMap dataset introduced by Yeh et al. (“Using publicly available satellite imagery and deep learning to understand economic well-being in Africa”). This dataset contains multispectral Landsat imagery paired with a continuous-valued asset wealth index, with natural distribution shifts arising from geographic (cross-country) and urban–rural differences. We have set up access via Google Earth Engine and are awaiting approval; results will be appended once available.

This setting is analogous to iWildCam but with continuous labels, for which we take the residual $s(X_n, Y_n) = |Y_n - f(X_n)|$ as the score. Based on our theoretical framework, we expect LP-Robust Conformal to perform well in this domain.

---

Questions 2+3. Thank you for the suggestion to use a different letter for the local perturbation in Proposition 2.3. We have double‑checked and confirm that Definition 3.2 is correct as stated and does not contain a typo.

---

Question 4. Thank you for the question. For smaller values of $\epsilon$, we observe a sharp increase in the corresponding $\rho$. This still yields valid coverage (consistent with the trade‑off described in Remark A.1, where smaller $\epsilon$ typically corresponds to larger $\rho$). However, because $\rho$ shifts the quantile level directly (see Proposition 3.4), even moderate increases in $\rho$ have a stronger effect on coverage than changes in $\epsilon$. As a result, parameter pairs with very small $\epsilon$ and large $\rho$ tend to be less efficient and are not selected under the coverage‑efficiency trade‑off in Corollary 4.2.

In the plots, we truncated the $\epsilon$ search range to $[0.5,\cdot]$ purely for illustration and readability, since smaller values would make $\rho$ spike sharply and obscure the relevant Pareto‑optimal region. In the revised version, we can enlarge the window and include an explanation of this phenomenon for completeness.

We are very grateful for your time and effort in reviewing our paper. Thank you for the positive assessment of our work. Below, we first provide answers to your questions and then to the highlighted weaknesses.

---

Questions 1+3. Thank you for these questions. Intuitively, the local perturbation parameter $\epsilon$ captures small, dense, and continuous shifts (such as Gaussian noise, mild input drift, or other shape-preserving distortions), while the global perturbation parameter $\rho$ models sparse or structural changes (such as label flips, adversarial attacks, or support mismatches). This decomposition, formalized in Propositions 2.1 and 2.3, allows the LP-bounded shift model to represent both small-magnitude changes and large, isolated deviations. Many realistic shifts (bounded $\ell_{\infty}$ attacks, mild covariate drift, or mixed noise-plus-outlier scenarios) fall within this class.

In practice, $(\epsilon,\rho)$ can be approximated from a modest number of calibration and test samples using the procedure in Appendix B (Algorithm 1, which we will move to the main text for clarity). The method proceeds in two steps:

(i) grid search over $\epsilon$, and for each candidate $\epsilon$ estimate $\rho$ as the Lévy-Prokhorov (Wasserstein-$\infty$) distance between independent calibration and test score distributions, computed via 1D optimal transport with cost $\mathbb{1}(|x-y| \geq \epsilon)$;

(ii) among candidate pairs, select the one minimizing the worst‑case $(1-\alpha)$-quantile (Corollary 4.2) to balance robustness and prediction-set size. Calibration scores used for parameter estimation are kept independent from those used to compute the conformal quantile to preserve validity.

This process reveals a natural trade-off: smaller $\epsilon$ typically corresponds to larger $\rho$, and vice versa. As shown in Figures 5 and 6, Pareto-optimal pairs (small $\epsilon$ with moderate $\rho$) achieve the best robustness-efficiency balance. The parameter $\epsilon$ affects set size additively, while $\rho$ shifts the quantile level itself (Proposition 3.4), making underestimation of $\rho$ riskier for coverage than underestimating $\epsilon$.

---

Question 2. Thank you for pointing this out. In Figure 1, the red curve represents the nominal score distribution $F_{\mathbb{P}}$, while the blue curve shows an adversarial distribution in the LP ambiguity set $\mathbb{B}_{\epsilon,\rho}(\mathbb{P})$, chosen to either maximize the $\beta$-quantile (left) or minimize coverage at a fixed threshold $q$ (right). The figure illustrates exactly the construction of the worst‑case distribution used in the proofs of Propositions 3.4 and 3.5, where it is cited for intuition. We will clarify both the meaning of the curves and this connection in the caption.

---

Weakness. We agree that broader empirical validation would further strengthen the paper. While synthetic shifts are useful for isolating and visualizing the effects of local and global perturbations, we also evaluate on iWildCam, a benchmark dataset with naturally occurring covariate shift. Our estimated parameters align closely with the Pareto frontier in Figure 4 and yield tight prediction sets, supporting the method’s practical utility.

In response to your suggestion, we are adding experiments on the PovertyMap dataset introduced by Yeh et al. (“Using publicly available satellite imagery and deep learning to understand economic well-being in Africa”). This dataset contains multispectral Landsat imagery paired with a continuous-valued asset wealth index, with natural distribution shifts arising from geographic (cross-country) and urban-rural differences. We have set up access via Google Earth Engine and are awaiting approval; results will be appended once available.

This setting is analogous to iWildCam but with continuous labels, for which we take the residual $s(X_n, Y_n) = |Y_n - f(X_n)|$ as the score. Based on our theoretical framework, we expect LP-Robust Conformal to perform well in this domain.

We are very grateful for your time and effort in reviewing our paper. Thank you for the positive assessment of our work. Below, we first provide answers to your questions and then to the highlighted weaknesses.

---

Question 1. Thank you for highlighting this important aspect. Intuitively, the local perturbation parameter $\epsilon$ captures small, dense, and continuous shifts (e.g., Gaussian noise or input drift) that alter the distribution’s shape, while the global perturbation parameter $\rho$ accounts for sparse or structural changes (e.g., label flips, adversarial attacks, or support mismatches). This decomposition is formally described in Propositions 2.1 and 2.3.

For estimation in practice, Appendix B (Algorithm 1) describes an efficient two‑step procedure using a small number of calibration and test samples:

(i) Grid search over $\epsilon$: For each candidate $\epsilon$, estimate $\rho$ as the LP distance between independent calibration and test score distributions, computed via 1D optimal transport with cost $\mathbb{1}(|x - y| \geq \epsilon)$. Each $(\epsilon_i, \rho_i)$ defines a valid LP ambiguity set satisfying empirical coverage.

(ii) Selection for efficiency: Among candidate pairs, choose the one minimizing the worst‑case $(1-\alpha)$-quantile (Corollary 4.2) to balance robustness and set size. Calibration scores for parameter estimation are kept independent from those used to compute the conformal quantile to ensure validity.

This yields a natural trade-off: smaller $\epsilon$ typically corresponds to larger $\rho$, and vice versa. As shown in Figures 5 and 6, Pareto-optimal pairs (small $\epsilon$ with moderate $\rho$) offer the best balance between robustness and efficiency. Importantly, $\epsilon$ contributes additively to the prediction set size, whereas $\rho$ shifts the quantile level itself (see the worst-case quantile characterization in Proposition 3.4) and therefore has a more pronounced effect on coverage. Consequently, underestimating $\rho$ poses a greater risk (potentially violating coverage) while underestimating $\epsilon$ primarily impacts efficiency.

---

Question 2. The LP‑bounded shift model is intentionally designed to be flexible: it only constrains most of the probability mass to remain within a small tolerance while allowing a limited fraction to be displaced arbitrarily. This enables it to capture both smooth, small‑magnitude perturbations and large, sparse deviations, covering a wide range of realistic shifts: from gradual covariate drift to adversarial perturbations and isolated but severe outliers. Many common real‑world shifts (e.g., bounded $\ell_{\infty}$ attacks, mild input drift, or mixed noise‑plus-outlier scenarios) fall within this class.

While exact verification would require the full test distribution, the parameters $(\varepsilon,\rho)$ can be approximated reliably from a modest number of calibration and test samples. In our experiments, we estimate them by comparing the empirical score distributions via a simple 1D optimal transport calculation. Conservative estimates ensure that the method’s coverage guarantee remains valid, and the results show that these estimates are accurate enough to retain efficiency.

We agree that parameter estimation is critical and have developed an explicit procedure (Algorithm 1 in Appendix B). In the revision, we will move this discussion to the main text (Section 4) for greater visibility, add further intuition and empirical support, and release our implementation to facilitate reproducibility.

---

Weaknesses. Answer to Weaknesses: We agree that broader empirical validation would further strengthen the paper. While synthetic shifts are useful for isolating and visualizing the effects of local and global perturbations, we also evaluate on iWildCam, a benchmark dataset with naturally occurring covariate shift. Our estimated parameters align closely with the Pareto frontier in Figure 4 and yield tight prediction sets, supporting the method’s practical utility.

In response to your suggestion, we are adding experiments on the PovertyMap dataset introduced by Yeh et al. (“Using publicly available satellite imagery and deep learning to understand economic well-being in Africa”). This dataset contains multispectral Landsat imagery paired with a continuous-valued asset wealth index, with natural distribution shifts arising from geographic (cross-country) and urban-rural differences. We have set up access via Google Earth Engine and are awaiting approval; results will be appended once available.

This setting is analogous to iWildCam but with continuous labels, for which we take the residual $s(X_n, Y_n) = |Y_n - f(X_n)|$ as the score. Based on our theoretical framework, we expect LP-Robust Conformal to perform well in this domain.

To provide additional results on real-world distribution shift, particularly for regression tasks, we aimed to include the PovertyMap dataset in our experiments. However, accessing PovertyMap requires an approval process that may take several weeks. To avoid delaying the response, we selected the Shifts Weather Prediction dataset instead, which offers a high-quality regression task under naturally occurring covariate and domain shifts. We believe this dataset serves as a strong and realistic benchmark, and complements our existing evaluation by further demonstrating the robustness and generalizability of our method in real-world settings.

Summary of the dataset and regression task:

The Shifts Weather Prediction dataset defines a scalar regression task: predicting surface temperature from 123 heterogeneous meteorological covariates and 4 metadata attributes (latitude, longitude, time, and climate type). It contains 10 million entries collected from September 2018 to September 2019 across diverse climate zones, reflecting real-world distribution shifts. A canonical split defines in-domain data (Tropical, Dry, Mild Temperate) and out-of-domain test sets (Snow, Polar), enabling robust evaluation under temporal and geographic variation. In our experiments, we predict the temperature using a CatBoostRegressor ensemble with 10 members as the forward model, and use the absolute error between the predicted and true temperature as the conformity score. To assess robustness in this setting, we applied our proposed LP-robust conformal prediction method, alongside five benchmark approaches including standard split conformal prediction and distribution shift–aware baselines. We estimated the $(\varepsilon, \rho)$ parameters in our method using the algorithm described in Appendix A.

Explanation of results:

Classical conformal prediction (SC), covariate shift–aware weighting (Weight), and the $f$-divergence–based method ($\chi^2$) all exhibit undercoverage, as they fail to account for global distribution shifts. Among the methods achieving valid coverage—RSCP, fine-grained CP (FG-CP), and our LP-robust approach—ours yields the smallest prediction intervals. RSCP produces wider sets potentially due to smoothed scores, and FG-CP slightly overcovers due to challenges in estimating conditional divergences. In contrast, our method achieves both valid and sharp coverage, suggesting that the learned $(\varepsilon, \rho)$ pair effectively captures the shift.

We are very grateful for your time and effort in reviewing our paper. Thank you for the positive assessment of our work. Below, we provide an answer to each of your questions (which also provide an answer to the weaknesses).

---

Question 1. Thank you for this insightful question. The standard split conformal prediction framework relies on the assumption that the $n$ calibration scores and the $(n+1)$‑th test score are exchangeable, which in particular requires that they are drawn from the same distribution (though they may be correlated). Under LP distribution shifts, however, the distribution of the test score necessarily differs from that of the calibration scores, violating exchangeability and rendering the standard method inapplicable. In order to address this more general setting, we require the assumption that the calibration scores are i.i.d. from the calibration distribution. This assumption seems technically essential, as it enables us to control the discrepancy between the empirical calibration distribution and its true, unknown counterpart, a property that is used explicitly in the final step of the proof of Theorem 4.1.

---

Question 2. Following the answer provided to the first question, the assumption of exchangeability between the $n$ calibration scores and the $(n+1)$‑th test score is strictly stronger than the LP distribution shift assumption that we impose (indeed, if the scores are exchangeable, then LP distance between the calibration and test distributions must be zero). In fact, since the test distribution and the calibration distribution are different, the standard proof reasoning based on rank statistics fails in the distribution shift case. We hope that this answers the reviewer's question, but are happy to provide more details in case our answer is not fully clear.

---

Question 3. Thank you for this sharp question. As noted in Remark A.1 (Appendix A of the submitted version), we discussed the sensitivity of the coverage guarantee in Theorem 4.1 to the parameters $\varepsilon$ and $\rho$. However, we acknowledge that this remark does not directly address the statement in question 3 that “to get a non‑trivial result, one needs to assume that $\rho < \alpha$.” We emphasize that $\rho < \alpha$ is not an assumption required by any of our technical results. Rather, it excludes the degenerate scenario in which the confidence interval covers the entire real line. Indeed, if $\rho \geq \alpha$, one can construct a sequence of distributions whose $(1-\alpha)$‑quantiles diverge to infinity (see Remark 3.3 in the submitted version). Thus, the condition $\rho < \alpha$ is not a limitation of our analysis, but a reflection of the phenomenon itself.

We also note that the condition $\rho < \alpha$ does not imply that the total variation distance is less than $\alpha$. Such an implication would hold only in the special case $\varepsilon = 0$, since any local perturbation also contributes to the total variation distance.

---

Question 4. Thank you for the opportunity to clarify this point. Among all conformal methods that achieve robustness in the LP sense, we claim that ours attains the highest possible efficiency. This follows directly from the exact (i.e., approximation-free) results in Propositions 3.4 and 3.5, which characterize the worst‑case quantile and coverage, respectively. These exact characterizations imply that the inequality in the first display equation of the proof of Theorem 4.1 becomes an equality in the worst case. The only source of approximation in Theorem 4.1 arises in its final step, where we relate the $\beta$‑quantiles of the empirical calibration distribution to those of the true calibration distribution. This approximation appears technically unavoidable in the absence of full knowledge of the true calibration distribution.