PAC PREDICTION SETS UNDER COVARIATE SHIFT
Year: 2021

Abstract:
An important challenge facing modern machine learning is how to rigorously quantify the uncertainty of model predictions. Conveying uncertainty is especially important when there are changes to the underlying data distribution that might invalidate the predictive model. Yet, most existing uncertainty quantification algorithms break down in the presence of such shifts. We propose a novel approach that addresses this challenge by constructing probably approximately correct (PAC) prediction sets in the presence of covariate shift. Our approach focuses on the setting where there is a covariate shift from the source distribution (where we have labeled training examples) to the target distribution (for which we want to quantify uncertainty). Our algorithm assumes given importance weights that encode how the probabilities of the training examples change under the covariate shift. In practice, importance weights typically need to be estimated; thus, we extend our algorithm to the setting where we are given confidence intervals for the importance weights. We demonstrate the effectiveness of our approach on covariate shifts based on DomainNet and ImageNet. Our algorithm satisfies the PAC constraint, and gives prediction sets with the smallest average normalized size among approaches that always satisfy the PAC constraint.

heading: 1 INTRODUCTION
text: A key challenge in machine learning is quantifying prediction uncertainty. A promising approach is via prediction sets, where the model predicts a set of labels instead of a single label. For example, prediction sets can be used by a robot to navigate to avoid regions where the prediction set includes an obstacle, or in healthcare to notify a doctor if the prediction set includes a problematic diagnosis.
Various methods based on these approaches can provide probabilistic correctness guarantees (i.e., that the predicted set contains the true label with high probability) when the training and test distributions are equal—formally, assuming the observations are exchangeable (Vovk et al., 2005; Papadopoulos et al., 2002; Lei et al., 2015) or i.i.d. (Vovk, 2013; Park et al., 2020a; Bates et al., 2021). However, this assumption often fails to hold in practice due to covariate shift—i.e., where the input distribution changes but the conditional label distribution remains the same. These shifts can be both natural (e.g., changes in color and lighting) (Hendrycks & Dietterich, 2019) or adversarial (e.g., `∞ attacks) (Szegedy et al., 2014).
We consider unsupervised domain adaptation (Ben-David et al., 2007), where we are given labeled examples from the source domain, but only unlabeled examples from the target (covariate shifted) domain. We propose an algorithm that constructs probably approximately correct (PAC) prediction sets under bounded covariate shifts (Valiant, 1984)—i.e., with high probability over the training data (“probably”), the prediction set contains the true label for test instances (“approximately correct”).
Our algorithm uses importance weights to capture the likelihood of a source example under the target domain. When the importance weights are known, it uses rejection sampling (von Neumann, 1951) to construct the prediction sets. Often, the importance weights must be estimated, in which case we need to account for estimation error. When only given confidence intervals around the importance weights, our algorithm constructs prediction sets that are robust to this uncertainty.
We evaluate our approach in two settings. First, we consider rate shift, where the model is trained on a broad domain, but deployed on a narrow domain—e.g., an autonomous car may be trained on both daytime and nighttime images but operate at night. Even if the model continues to perform well, this covariate shift may invalidate the prediction sets. We show that our approach constructs valid
prediction sets under rate shifts constructed from DomainNet (Peng et al., 2019), whereas several existing approaches do not. Second, we consider support shift, where the model is trained on one domain, but evaluated on another domain with shifted support. We train ResNet101 (He et al., 2016) using unsupervised domain adaptation (Ganin et al., 2016) on both ImageNet-C synthetic perturbations (Hendrycks & Dietterich, 2019) and PGD adversarial perturbations (Madry et al., 2017) to ImageNet (Russakovsky et al., 2015). Our PS-W algorithm satisfies the PAC constraint, and gives prediction sets with the smallest average normalized size among approaches that always satisfy the PAC constraint. See Figure 1 for a summary of our approach and results.
Related work. Our work is related to results on conformal prediction (Vovk et al., 2005; Balasubramanian et al., 2014) without a shift, where the goal is to construct models that predict sets of labels designed to include the true label with high probability (instead of predicting individual labels) under the assumption that the source and target distributions are the same. In particular, our setting is related to inductive (or split) conformal prediction (Papadopoulos et al., 2002; Lei et al., 2015), where the training set is split into a proper training set used to train a traditional predictive model, and a calibration set used to construct the prediction sets. The general approach is to train a model f : X ×Y → R that outputs conformity scores, and then to choose a threshold τ ∈ R that satisfies a correctness guarantee, where the corresponding prediction sets are C(x) = {y ∈ Y | f(x, y) ≥ τ} (Park et al., 2020a; Gupta et al., 2021); this notation is defined in Section 2.
Several kinds of correctness guarantees under no shift have been considered. One possibility is conditional (Barber et al., 2020) or object conditional (Vovk, 2013) validity, which ensures correctness for all future covariates x, with high probability only over the conditional label distribution p(y | x). This guarantee is very strong, and hard to ensure in practice. A weaker notion is approximate (Barber et al., 2020) or group (Romano et al., 2019) conditional validity, which ensures correctness with high probability over p(y | x) as well as some distribution p(x) over a subgroup. Finally, unconditional validity ensures only correctness over the joint distribution p(x, y). We focus on unconditional validity, though our approach extends straightforwardly to group conditional validity.
A separate issue, arising in the no shift setting, is how to condition on the calibration set Z. A conventional approach is to include the randomness in Z with the above, e.g., for unconditional validity, the high-probability guarantee is now over p(x, y, Z). We refer to this strategy as fully unconditional validity. However, the guarantee we consider uses a separate confidence level for the training data, which is called a training conditional guarantee (Vovk, 2013); this correctness notion is equivalent to a PAC correctness guarantee (Park et al., 2020a), and is also related to the notion of tolerance regions from statistics (Wilks, 1941). We build on Park et al. (2020a), which formulates the problem of choosing τ as learning a binary classifier where the input and parameter spaces are both one-dimensional; thus, the correctness guarantee corresponds to a PAC generalization bound. This approach is widely applicable since it can work with a variety of objectives (Bates et al., 2021).
Recent work has extended inductive conformal prediction to a setting with covariate shift (Tibshirani et al., 2019; Lei & Candès, 2020); similarly, Podkopaev & Ramdas (2021) considers conformal prediction under label shift (Lipton et al., 2018), i.e., assuming the conditional probabilities p(x | y) do not change. These approaches assume that the true importance weights (IWs) are known (Tibshirani et al., 2019; Podkopaev & Ramdas, 2021), or assume some properties of the estimated IWs (Lei & Candès, 2020), whereas our approach considers uncertainty in the estimated IWs. Furthermore, they are focused on fully unconditional validity, whereas we obtain PAC prediction sets. In addition, Cauchois et al. (2020) designs confidence sets that are robust to all distribution shifts with bounded f -divergence; in contrast, we consider the unsupervised learning setting where we have unlabeled examples from the target distribution. We provide additional related work in Appendix A.

heading: 2 BACKGROUND ON PAC PREDICTION SETS
text: We give background on PAC prediction sets (Park et al., 2020a); we also re-interpret this approach using Clopper-Pearson confidence intervals (Clopper & Pearson, 1934), which aids our analysis.

heading: 2.1 PAC PREDICTION SETS ALGORITHM
text: Let x ∈ X be covariates and y ∈ Y be labels. We consider a source distribution P over X ×Y with probability density function (PDF) p(x, y).1 A prediction set2 is a set-valued function C : X → 2Y . Inputs. We are given a held-out calibration set Sm of i.i.d. samples (xi, yi) ∼ P for i ∈ [m] := {1, . . . ,m}, written as Sm ∼ Pm, and a score function f : X × Y → R≥0. For example, f(x, y) can be a prediction for the probability that y is the label for x. The score function can be arbitrary, but score functions assigning higher scores to likely labels yield smaller prediction sets.
PAC prediction set. A PAC prediction set is a set-valued function C : X → 2Y satisfying the following two conditions. First, C is approximately correct in the sense that its expected error (failing to contain the true label) is bounded by a given ε ∈ (0, 1), i.e.,
LP (C) := E(x,y)∼P [1(y /∈ C(x))] = P(x,y)∼P [y /∈ C(x)] ≤ ε. Second, noting that CSm is constructed based on a calibration set Sm ∼ Pm, the condition that CSm is approximately correct must be satisfied with high probability—i.e., for given δ ∈ (0, 1), we have
PSm∼Pm [LP (CSm) ≤ ε] ≥ 1− δ. Our goal is to devise an algorithm for constructing a PAC prediction set C. Letting C(x) = Y for all x ∈ X always satisfies both conditions above, but this extreme case would be uninformative if used as a prediction set. Therefore, we additionally want to minimize the expected sizeE[S(C(x))] of the prediction sets C(x), where S : 2Y → R≥0 is a size measure, which is application specific (e.g., the cardinality of a set in classification); however, we only rely on the monotonicity of the size measure with respect to the prediction set parameterization in construction.
Algorithm. To construct C, we first define the search space of possible prediction sets along with the size measure S. We parameterize C by a scalar τ ∈ T := R≥0 as
Cτ (x) = {y ∈ Y | f(x, y) ≥ τ}, i.e., τ is the threshold on f(x, y) above which we include y in C(x). Importantly, τ controls the tradeoff between size and expected error. The reason is that if τ1 ≤ τ2, then Cτ2(x) ⊆ Cτ1(x) for all x ∈ X . Thus, size is monotonically decreasing in τ—i.e., S(Cτ2(x)) ≤ S(Cτ1(x)) for all x ∈ X , and error is monotonically increasing in τ—i.e., LP (Cτ1) ≤ LP (Cτ2). See Figure 1b for an illustration, and Park et al. (2020a) and Gupta et al. (2021) for details.
As a consequence, a typical goal is to construct Cτ that provably contains the true label with high probability, while empirically minimizing size (Vovk et al., 2005; Gupta et al., 2021). In our setting, we want to maximize τ (equivalently, minimize expected size) subject to the constraint that
1All quantities that we define are measurable with respect to a fixed σ-algebra on X × Y; for instance, p is the density induced by a fixed σ-finite measure. To be precise, we consider a probability measure P defined with respect to the base measureM on X ×Y; then, p = dP/dM is the Radon-Nykodym derivative of P with respect to M . For classification, M is the product of a Lebesgue measure on X and a counting measure on Y .
2We use the term “prediction set” to denote both the set-valued function and a set output by this function.
Cτ is PAC. Let L̄Sm(Cτ ) := ∑
(x,y)∈Sm 1(y /∈ Cτ (x)) be the empirical error count on Sm, and F (k;m, ε) = ∑k i=0 ( m k ) εi(1−ε)m−i be the cumulative distribution function (CDF) of the binomial distribution Binom(m, ε) with m trials and success probability ε. In prior work, Park et al. (2020a) constructs Cτ̂ by solving the following problem:
τ̂ = max τ∈T τ subj. to L̄Sm(Cτ ) ≤ k(m, ε, δ), (1)
where
k(m, ε, δ) = max k∈N∪{0}
k subj. to F (k;m, ε) ≤ δ.
Intuitively, the PAC guarantee via this construction is related to the Binomial distribution; L̄Sm(C) has distribution Binom(m,LP (C)) (given a fixed C), since 1(y 6∈ C(x)) has a Bernoulli(LP (C)) distribution when (x, y) ∼ P . Thus, k(m, ε, δ) defines a “confidence interval” such that if L̄Sm(C) ≤ k(m, ε, δ), then LP (C) ≤ ε with probability at least 1 − δ. Below, we formalize this intuition by drawing a connection to the Clopper-Pearson confidence interval.

heading: 2.2 CLOPPER-PEARSON INTERPRETATION
text: We interpret (1) using the Clopper-Pearson (CP) upper bound θ(k;m, δ) ∈ [0, 1] (Clopper & Pearson, 1934; Park et al., 2021). This is an upper bound on the true success probability µ, constructed from a sample k ∼ Binom(m,µ), which holds with probability at least 1 − δ, i.e., Pk∼Binom(m,µ)[µ ≤ θ(k;m, δ)] ≥ 1− δ, where
θ(k;m, δ) := inf {θ ∈ [0, 1] | F (k;m, θ) ≤ δ} ∪ {1}. (2)
Intuitively, for a fixed C, L̄Sm(C) ∼ Binom(m,LP (C)), so the true error LP (C) is contained in the CP upper bound θ(LSm(C);m; δ) with probability at least 1 − δ. Intuitively, we can therefore choose τ so that this upper bound is ≤ ε. However, we need to account for generalization error of our estimator. To this end, we have the following (see Appendix D.1 for a proof and Algorithm 2 in Appendix E for implementation details on the corresponding algorithm):
Theorem 1 Let UCP(C, Sm, δ) := θ(L̄Sm(C);m; δ), where θ is defined in (2). Let τ̂ be the solution of the following problem3:
τ̂ = max τ∈T
τ subj. to UCP(Cτ , Sm, δ) ≤ ε. (3)
Then, we have PSm∼Pm [LP (Cτ ) ≤ ε] ≥ 1− δ for any τ ≤ τ̂ .

heading: 3 PAC PREDICTION SETS UNDER COVARIATE SHIFT
text: We extend the PAC prediction set algorithm described in Section 2 to the covariate shift setting. Our novel approach combines rejection sampling (von Neumann, 1951) with Theorem 1.

heading: 3.1 PROBLEM FORMULATION
text: We assume labeled training examples from the source distribution P are given, but want to construct prediction sets that satisfy the PAC property with respect to a (possibly) shifted target distributionQ. Both of these are distributions over X × Y; denote their PDFs by p(x, y) and q(x, y), respectively. The challenge is that we are only given unlabeled examples from Q.
Preliminaries and assumptions. We denote the likelihood ratio of covariate distributions by w∗(x) := q(x)/p(x), also called the importance weight (IW) of x. Our main assumption is the covariate shift condition, which says the label distributions are equal (i.e., p(y | x) = q(y | x)), while the covariate distributions may differ (i.e., we can have p(x) 6= q(x)). Inputs. We assume a labeled calibration set Sm consisting of i.i.d. samples (xi, yi) ∼ P (for i ∈ [m]) is given, and a score function f : X × Y → R≥0. For now, we also assume we have
3We consider a trivial solution τ = 0 when (3) is not feasible, but omitting here to avoid clutter.
the true importance weights w∗i := w ∗(xi) for each (xi, yi) ∼ P , as well as an upper bound b ≥ maxx∈X w∗(x) on the IWs. In Sections 3.3 and Appendix B, we describe how to estimate these quantities in a way that provides guarantees under smoothness assumptions on the distributions.4
Problem. Our goal is to construct CSm that is a PAC prediction set under Q—i.e.,
PSm∼Pm [LQ(CSm) ≤ ε] ≥ 1− δ,
where LQ(C) := P(x,y)∼Q[y 6∈ C(x)] is the error of C on Q, while empirically controlling its size.

heading: 3.2 REJECTION SAMPLING STRATEGY
text: Our strategy is to use rejection sampling to convert Sm consisting of i.i.d. samples from P into a labeled calibration set consisting of i.i.d. samples from Q.
Rejection sampling. Rejection sampling (von Neumann, 1951; Owen, 2013; Rubinstein & Kroese, 2016) is a technique for generating samples from a target distribution q(x) based on samples from a proposal distribution p(x). Given importance weights (IWs) w∗(x) and an upper bound b ≥ maxx∈X w
∗(x), it constructs a set of i.i.d. samples from q(x). Typically, rejection sampling is used when the proposal distribution is easy to sample compared to the target distribution. In contrast, we use it to convert samples from the source distribution into samples from the target distribution.
In particular, our algorithm takes the proposal distribution to be the source covariate distribution PX , and the target distribution to be our target covariate distributionQX . Let Vi ∼ U := Uniform([0, 1]) be i.i.d., V := (V1, . . . , Vm), and ~w∗ := (w∗1 , . . . , w ∗ m) with w ∗ i = w
∗(xi) as defined before. Then, given Sm, it uses rejection sampling to construct a randomly chosen set of N samples
TN (Sm, V, ~w ∗, b) := { (xi, yi) ∈ Sm ∣∣∣∣ Vi ≤ w∗ib }
(4)
from QX . The expected number of samples is E[N ] = m/b; thus, rejection sampling is more effective when the proposal distribution is similar to the target.
Rejection sampling Clopper-Pearson (RSCP) bound. Given TN , our algorithm uses the CP bound to construct a PAC prediction set C for Q. Let σi := 1(Vi ≤ w∗i /b) indicate whether example (xi, yi) ∈ Sm is accepted—i.e., TN (Sm, V, ~w∗, b) = {(xi, yi) ∈ Sm | σi = 1},where |TN (Sm, V, ~w∗, b)| = ∑m i=1 σi. Then, it follows that URSCP bounds the error on Q, where
URSCP(C, Sm, V, ~w ∗, b, δ) := UCP(C, TN (Sm, V, ~w ∗, b), δ). (5)
Thus, we have the following (see Appendix D.3 for the proof, and Algorithm 4 in Appendix E for a detailed description of the PS-R algorithm that leverages this bound):
Theorem 2 Define URSCP as in (5). Let τ̂ be the solution of the following problem:
τ̂ = max τ∈T
τ subj. to URSCP(Cτ , Sm, V, ~w∗, b, δ) ≤ ε. (6)
Then, we have PSm∼Pm,V∼Um [LQ(Cτ ) ≤ ε] ≥ 1− δ for any τ ≤ τ̂ .
Note that V is sampled only once for Algorithm 4, so the randomness in V can contribute to the failure of the correctness guarantee; however, the failure rate is controlled by δ.

heading: 3.3 APPROXIMATE IMPORTANCE WEIGHTS
text: So far, we have assumed that the true importance weight w∗(x) is known. Since in practice, it needs to be estimated, we relax this assumption to only needing an uncertainty set of possible importance weights. This allows us to handle estimation error in the weights.
Problem. Let SXm be unlabeled calibration examples from the source distribution (i.e., the covariates in Sm), TXn be n unlabeled calibration examples from the target distribution (denoted by T X n ∼
4In practice, we can improve stability of importance weights by considering source and target distributions induced by a feature mapping, e.g., learned using unsupervised domain-adaptation (Park et al., 2020b).
QnX ), and ~w ∗ := (w∗1 , ..., w ∗ m) ∈ Rm be the vector of true importance weights w∗i := w∗(xi), for (xi, yi) ∈ Sm. Then, we assume an uncertainty setW ⊆ Rm that has ~w∗ with high probability, i.e., PSXm∼PmX ,TXn ∼QnX [~w
∗ ∈ W] ≥ 1− δw, (7) where δw ∈ (0, 1). We assumeW has the form W := {w ∈ Rm | ∀i ∈ [m] , wi ≤ wi ≤ wi} , for some wi and wi. See the discussion on our choice ofW in Appendix C.1. Robust Clopper-Pearson bound. To construct a PAC prediction set C for Q, it suffices to bound the worst-case error over w ∈ W , i.e., we have the following (see Appendix D.4 for the proof):
Theorem 3 SupposeW satisfies (7). Define URSCP as in (5). Let τ̂ be the solution of the following problem:
τ̂ = max τ∈T τ subj. to max w∈W URSCP(Cτ , Sm, V, w, b, δC) ≤ ε. (8)
Then, we have PSm∼Pm,V∼Um,TXn ∼QnX [LQ(Cτ ) ≤ ε] ≥ 1− δC − δw for any τ ≤ τ̂ .
A key challenge applying Theorem 3 is solving the maximum over w ∈ W . We propose a simple greedy algorithm that achieves the maximum.
Greedy algorithm for robust URSCP. The RSCP bound URSCP satisfies certain monotonicity properties that enable us to efficiently compute an upper bound to the maximum in (8). In particular, if C makes an error on (xi, yi) (i.e., yi 6∈ C(xi)), then URSCP is monotonically non-decreasing in w∗i = w
∗(xi); intuitively, this holds since a larger w∗i increases the probability that (xi, yi) is included in TN (Sm, V, w∗, b), which in turn increases the empirical error L̄TN (Sm,V,w∗,b)(C). Conversely, if C does not make an error on (xi, yi) (i.e., yi ∈ C(xi)), URSCP is non-increasing in w∗i . More formally, we have the following result (see Appendix D.5 for a proof):
Lemma 1 For any i ∈ [m], URSCP(C, Sm, V, w∗, b, δ) is monotonically non-decreasing in w∗i if yi /∈ C(xi), and monotonically non-increasing in w∗i if yi ∈ C(xi).
Our greedy algorithm leverages the monotonicity of URSCP. In particular, givenW and C, the choice ŵ := (ŵ1, . . . , ŵm), where ŵi = { wi if yi 6∈ C(xi) wi if yi ∈ C(xi) (∀i ∈ [m]) (9)
is the maximum value over w ∈ W of the constraint URSCP(Cτ , Sm, V, w, b, δC) in (8), i.e., max w∈W URSCP(C, Sm, V, w, b, δ) = URSCP(C, Sm, V, ŵ, b, δ). (10)
Thus, we have the following, which follows by (10) and the same argument as the Theorem 3:
Theorem 4 SupposeW satisfies (7). Define ŵτ,SXm,TXn as in (9), making the dependency on τ , S X m , and TXn explicit, and URSCP as in (5). Let τ̂ be the solution of the following problem: τ̂ = max
τ∈T τ subj. to URSCP(Cτ , Sm, V, ŵτ,SXm,TXn , b, δC) ≤ ε. (11)
Then, we have PSm∼Pm,V∼Um,TXn ∼QnX [LQ(Cτ ) ≤ ε] ≥ 1− δC − δw for any τ ≤ τ̂ . 5
Importance weight estimation. In general, to estimate the importance weights (IWs), some assumptions on their structure are required. A number of approaches have been proposed, with varying guarantees under different assumptions (Kanamori et al., 2009; Cortes et al., 2008; Nguyen et al., 2010; Lipton et al., 2018). We use a cluster-based approach (Cortes et al., 2008). In particular, we heuristically partition the feature space of the score function f into K bins by using a probabilistic classifier g that separates source and target examples, and then estimate the source and target covariate distributions p(x) and q(x) based on the smoothness assumption over the distributions, where the degree of the smoothness is parameterized by E. Then, we can construct the uncertainty setW that satisfies the specified guarantee in (7). See Appendix B for details.
Algorithm. Our algorithm, called PS-W, is detailed in Algorithm 1; it solves (20) and also performs importance weight estimation according to (17) in Appendix B. In particular, Algorithm 1 constructs a prediction set that satisfies the PAC guarantee in Theorem 4. See Section 4.1 for our choice of parameters (e.g., K, E, and grid search parameters).
5We assume b is given for simplicity, but our approach estimates it; see Appendix C.2 for details.
Algorithm 1 PS-W: an algorithm using the robust RSCP bound in (20)
1: procedure PS-W(Sm, TXn , f, g, T , ε, δw, δC ,K,E) 2: W ← ESTIMATEIWS(SXm , TXn , g, δw,K,E) (.) Compute uncertainty setW 3: b← maxi∈[K] wi (.) Compute maximum IW (see (19) in Appendix C.2) 4: return PS-ROBUST(Sm, f, T ,W, b, ε, δC) 5: procedure ESTIMATEIWS(SXm , TXn , g, δw,K,E) 6: Construct bins B1, . . . , BK using g as described in Appendix B.2 7: ConstructW using B1, . . . , BK , SXm , TXn , g, δw and, E as described in Appendix B.1 8: returnW 9: procedure PS-ROBUST(Sm, f, g, T ,W, b, ε, δC) 10: V ∼ Uniform([0, 1])m 11: τ̂ ← 0 12: for τ ∈ T do (.) Grid search in ascending order 13: Construct ŵ using (9) given τ , Sm, andW 14: if URSCP(Cτ , Sm, V, ŵ, b, δC) ≤ ε then 15: τ̂ ← max(τ̂ , τ) 16: else 17: break 18: return τ̂

heading: 4 EXPERIMENTS
text: We show the efficacy of our approach on rate and support shifts on DomainNet and ImageNet.

heading: 4.1 EXPERIMENTAL SETUP
text: Models. For each source-target distribution pair, we split each the labeled source data and unlabeled target data into train and calibration sets. We use a separate labeled test set from the target for evaluation. For each shift from source to target, we use a deep neural network score function f based on ResNet101 (He et al., 2016), trained using unsupervised domain adaptation based on the source and target training sets. See Appendix F for details.
Prediction set construction. To construct our prediction sets, we first estimate IWs by training a probabilistic classifier g using the source and target training sets. Next, we use g to construct heuristic IWs w(x) = 1/g(s = 1 | x)− 11, where s = 1 if x is from source. Then, we estimate the lower and upper bound of the true IWs using Theorem 5 with E = 0.001 and K = 10 bins (chosen to contain equal numbers of heuristic IWs), where we compute the lower and upper Clopper-Pearson interval using the source and target calibration sets. Furthermore, given a confidence level δ, we use δC = δw = δ/2. For the grid search in line 12 of Algorithm 1, we increase τ by 10−7 until the bound URSCP exceeds 1.5ε. Finally, we evaluate the prediction set error on the labeled target test set.
Baselines. We compare the proposed method in Algorithm 1 (PS-W) with the following:
• PS: The prediction set using Algorithm 2 that satisfies the PAC guarantee from Theorem 1, based on Park et al. (2020a), which makes the i.i.d. assumption. • PS-C: A Clopper-Pearson method accounting for covariate shift by a conservative upper bound on the empirical loss (see Appendix E.2 for details), resulting in Algorithm 3 that solves the following problem:
τ̂ = arg max τ∈T
τ subj. to UCP(Cτ , Sm, δ) ≤ ε/b.
• WSCI: Weighted split conformal inference, proposed in (Tibshirani et al., 2019). Under the exchangeability assumption (which is somewhat weaker than our i.i.d. assumption), this approach provides correctness guarantees only for a single future test instance, i.e., it includes an ε confidence level (denoted α in their paper) but not δ. Also, it assumes the true IWs are known; following their experiments, we use the heuristic IWs constructed using g. • PS-R: The prediction set using Algorithm 4 that satisfies the PAC guarantee from Theorem 2 with the heuristic IWs constructed using a probabilistic classifier g.
Metrics. We measure performance via the prediction set error and size on a held-out test set—i.e., error is the fraction of (x, y) such that y 6∈ C(x) and size is the average of |C(x)| over x. Rate shifts on DomainNet. First, we consider settings where the model is trained on data from a variety of domains, but is then deployed on a specific domain; we call such a shift rate shift. For instance, a self-driving car may be trained on both day and night images, but tested during the night time. While the model should still perform well since the target is a subset of the source, the covariate shift can nevertheless invalidate prediction set guarantees. We consider rate shifts on DomainNet (Peng et al., 2019), which consists of images of 345 classes from six domains (sketch, clipart, painting, quickdraw, real, and infograph). We split the dataset into 409,832 training, 88,371 calibration, and 88,372 test images; we use all domains as the source and each domain as a target.
Support shifts on ImageNet. Next, we consider support shifts, where the support of the target is different from the support of the source, but unsupervised domain adaptation is used to learn feature representations that align the two distributions (Ganin et al., 2016); then, the score function f is trained on this representation. First, we consider ImageNet-C (Hendrycks & Dietterich, 2019), which modifies the original ImageNet dataset (Russakovsky et al., 2015) using 15 synthetic perturbations with 5 severity levels. We use 13 perturbations, omitting “snow” and “glass blur”, which are computationally expensive to run. We consider the original ImageNet dataset as the source, and the all synthetic perturbations on all of ImageNet (denoted ImageNet-C13) as the target. We split ImageNet into 1.2M training, 25K calibration, and 25K test images, and ImageNet-C13 into 83M training, 1.6M calibration, and 1.6M test images. Second, we consider adversarial shifts, where we generate adversarial examples for ImageNet using the PGD attack (Madry et al., 2017) with 0.01 `∞-norm perturbations with respect to a pretrained ResNet101. We consider the original ImageNet as the source and the adversarially perturbed ImageNet as the target.

heading: 4.2 EXPERIMENTAL RESULTS
text: We summarize our results in Table 1, and provide more details in Figure 4 in Appendix G.2. Our PS-W algorithm satisfies the PAC constraint, and gives prediction sets with the smallest average normalized size among approaches that always satisfy the PAC constraint.
Rate shifts on DomainNet. We use ε = 0.1 and δ = 10−5. As can be seen in Table 1, the prediction set error of our approach (PS-W) does not violate the error constraint ε = 0.1, while all other comparing approaches (except for PS-C) violate it for at least one of the shifts. While PS-C satisfies the desired bound, it is overly conservative; its prediction sets are significantly larger than necessary, making it less useful for uncertainty quantification. For the shift to Infograph in Table 1, PS-W achieves relatively larger average prediction set size compared to other shifts; this is because the classification error of the score function f over Infograph is large—71.28%, whereas that over Sketch is 37.16%. Even with the poor score function, our proposed approach still satisfies the ε = 0.1 error constraint. Finally, while our approach PS-W generally performs best subject to satisfying the error constraint, we note that our ablations PS-R and PS-M also perform well, providing different tradeoffs. First, the performance of PS-W is significantly more reliable than PS-R, but in some cases PS-R performs better (e.g., it produces slightly smaller prediction sets on rate shifts but significantly larger sets on support shifts). Alternatively, PS-M consistently produces smaller prediction sets, though it sometimes violates the ε error constraints.
Support shifts on ImageNet. We show results for synthetic and adversarial shifts in Table 1 (and Figure 4 in Appendix G.2). As can be seen, the error of our approach (PS-W) is below the desired level. PS-R performs poorly, likely due to the uncalibrated point IW estimation—we find that calibrated importance weights mitigate these issues, though accounting for uncertainty in the IWs is necessary for achieving the desired error rate; see Appendix G.3 for details.
For adversarial shifts, the target classification error of the source-trained ResNet101 and the domainadapted ResNet101 is 99.97% and 28.05%, respectively. Domain adaptation can significantly decrease average error rate, but label predictions still do not have guarantees. However, our prediction set controls the prediction set error rate as specified by ε. As shown in Figure 4, our prediction set function outputs a prediction set for a given example that includes the true label at least 90% of the time. Thus, downstream modules can rely on this guarantee for further decision-making.
Ablation and sensitivity study. We conduct an ablation study on the effect of IW calibration and the smoothness parameter E for PS-W. We observe that the IW calibration considering uncertainty intervals around calibrated IWs is required for the PAC guarantee (e.g., Figure 2b). Also, we find that a broad range of E-s can be used to satisfy the PAC guarantee; we believe this result is due to the use of a domain adapted score function. See Appendices G.3 & G.5 for details.

heading: 5 CONCLUSION
text: We propose a novel algorithm for building PAC prediction sets under covariate shift; it leverages rejection sampling and the Clopper-Pearson interval, and accounts for uncertain IWs. We demonstrate the efficacy of our approach on natural, synthetic, and adversarial covariate shifts on DomainNet and ImageNet. Future work includes providing guarantees on the prediction set size, rigorous estimation of the hyperparameter E, and incorporating probabilistic IW uncertainty estimates.
Ethics statement. We do not foresee any significant ethical issues with our work. One possible issue is that end-users may overly trust the prediction sets if they do not understand the limitations of our approach—e.g., it is only a high probability guarantee.
Reproducibility statement. Algorithms used for evaluation, including ours, are stated in Algorithm 1, Algorithm 2, Algorithm 3, Algorithm 4, and Algorithm 5, along with hyperparameters for algorithms in Section 4.1 and Appendix F. Related code will be released once it is accepted. The proof of our theorems are stated in Appendix D; the required assumption is stated in Assumption 1.

heading: A ADDITIONAL RELATED WORK
text: Prediction sets under i.i.d. assumption. We built our proposed approach on the known PAC prediction set approach (Park et al., 2020a) due to its simplicity and sample efficiency. As other candidates, nest conformal prediction (Gupta et al., 2021) reinterprets multiple known conformal prediction approaches using a scalar parameteriztion of prediction sets; but different to Park et al. (2020a) it does not provide PAC guarantee. Kivaranovic et al. (2020) provides a PAC style guarantee, but their approach is limited to regression and sample-inefficient, i.e., it requires O(1/ 2) sample, while Park et al. (2020a) has a clearly better sample complexity as demonstrated in their paper.
Prediction set under various settings. Prior prediction set algorithms have considered several settings. First, traditional conformal prediction algorithms (Vovk et al., 2005) often consider the sequential setting, where labeled examples arrive sequentially from the same distribution; there has also been work extending conformal prediction to the setting where the distribution is timevarying (Gibbs & Candès, 2021). Alternatively, there has been work on constructing risk-controlling prediction sets for supervised learning setting (Vovk, 2013; Park et al., 2020a; Bates et al., 2021). Finally, there has been recent work on conformal prediction in the meta-learning setting (Fisch et al., 2021); in particular, given a few labeled examples drawn from the new task, their approach leverages labeled examples from previous tasks to construct a conformal predictor for the new task, assuming the tasks are exchangeable.
Calibration. An alternative way to quantify uncertainty is calibrated prediction (Murphy, 1972; DeGroot & Fienberg, 1983; Guo et al., 2017), which aims to ensure that among instances with a predicted confidence p, the model is correct a fraction p of the time. Techniques have been proposed to re-scale predicted confidences to improve calibration (Platt, 1999; Guo et al., 2017; Zadrozny & Elkan, 2001; 2002; Kuleshov & Liang, 2015; Kuleshov et al., 2018; Malik et al., 2019); including ones with theoretical guarantees (Kumar et al., 2019; Park et al., 2021) and ones that handle covariate shift (Park et al., 2020b; Wang et al., 2020). These approaches provide a qualitatively different form of uncertainty quantification compared to the one we consider.
Rejection sampling. Rejection sampling (sometimes accept-reject sampling) is a well-known technique (Owen, 2013; Rubinstein & Kroese, 2016) dating back at least to von Neumann (1951). We highlight that most work on covariate shift relies on importance weighting; our rejection sampling approach is relatively less common and more novel. In fact, to the best of our knowledge, only Pagnoni et al. (2018) has used rejection sampling for a different problem in this area.
IW estimation. There has been a long line of work studying the problem of estimating importance weights (IWs), also called likelihood ratios, in a way that provides theoretical guarantees. For instance, (Nguyen et al., 2010) provides a convergence rate analysis of IW estimators—under smoothness assumptions, they provide a finite sample bound on the Hellinger distance between the true IW and estimated IW. Next, (Kanamori et al., 2009) shows a similar finite sample guarantee, assuming the true IW can be represented as a linear combination of kernels. Finally, (Cortes et al., 2008) proposes non-parametric IW estimators, modeling the source and target distribution by histograms over clusters in sample space.
The IW estimation approaches can be used in conjuction with prediction set construction. Compared to Candès et al. (2021), which only guarantees asymptotic validity under certain assumptions on the estimated IWs, Theorem 4 provides a significantly stronger finite-sample correctness guarantee (in addition to providing PAC-style guarantees). Furthermore, we explicitly describe an algorithm for approximate IWs, required for the validity of Theorem 4, in Appendix B.
B IMPORTANCE WEIGHT ESTIMATION
In general, to estimate the importance weights (IWs), some assumptions on their structure are required. A number of approaches have been proposed, with varying guarantees under different assumptions (Kanamori et al., 2009; Cortes et al., 2008; Nguyen et al., 2010; Lipton et al., 2018). We use a cluster-based approach (Cortes et al., 2008); our approach is compatible with any of these strategies if they can be modified to provide uncertainty estimates of the IWs.
In our approach, given a partition X = ⋃K j=1Bj into bins, we can estimate the IWs based on the fractions of source and target samples in each bin. If the partition is sufficiently fine, then we can
obtain confidence intervals around the estimated IWs with finite-sample guarantees. However, this strategy requires the number of bins in the partition to be exponential in the dimension of X . Thus, in practice, we use a heuristic to construct the partition. We describe the cluster-based approach and our partition construction heuristic below.

heading: B.1 CLUSTER-BASED APPROACH
text: We assume given unlabeled calibration sets SXm and T X n , where S X m consists of i.i.d. samples xi ∼ p(x) for i ∈ [m], and TXn consists of i.i.d. samples xi ∼ q(x) for i ∈ [n], respectively6. Roughly speaking, the cluster-based strategy estimates the average IW in each bin Bj ; assuming p and q are roughly constant in each bin, these accurately estimate the true IWs. Let j(x) be the bin containing x (i.e., x ∈ Bj(x)), and let
pB(x) := pj(x) s.t. pj = ∫ Bj p(x′) dx′ and qB(x) := qj(x) s.t. qj = ∫ Bj q(x′) dx′
be the (unnormalized) approximations of the densities p and q, respectively, that are constant on each bin. We assume that pB and qB are accurate approximations:
Assumption 1 Given E ∈ R≥0, the partition satisfies∫ Bj |p(x)− p(x′)|dx′ ≤ E and ∫ Bj |q(x)− q(x′)|dx′ ≤ E (j ∈ [K],∀x ∈ Bj). (12)
Thus, p and q are roughly constant on the partitions. In general, (12) can hold for any E ∈ R>0 if p and q are Lipschitz continuous and each Bj is sufficiently small (see Appendix B.3 for discussion). Then, under Assumption 1, it can be verified that
|v(x) · p(x)− pB(x)| ≤ E and |v(x) · q(x)− qB(x)| ≤ E (∀x ∈ X ), (13) where v(x) = vj(x), and vj = ∫ Bj
dx′ is the volume of bin Bj (see the proof of Theorem 5 for the validity of (13)). Next, we have the following empirical estimates of pB and qB , respectively:
p̂B(x) := p̂j(x) s.t. p̂j = 1
m ∑ x′∈SXm 1 (x′ ∈ Bj) and (14)
q̂B(x) := q̂j(x) s.t. q̂j = 1
n ∑ x′∈TXn 1 (x′ ∈ Bj) . (15)
Now, 1(x′ ∈ Bj) has distribution Bernoulli(pj) when x ∼ P , thus m · p̂j has distribution Binom(m, pj), and pj is contained in a Clopper-Pearson interval around p̂j with high probability; in particular, let θ be the Clopper-Pearson lower bound corresponding to the Clopper-Pearson upper bound defined in Section 2.2, i.e., Pk∼Binom(m,µ)[µ ≥ θ(k;m, δ)] ≥ 1− δ. Then, we have
θ(m · p̂j ;m, δ′) ≤ pj ≤ θ(m · p̂j ;m, δ′) (16)
with probability at least 1 − δ′ with respect to the samples SXm . Combining (13) and (16), we have the following result (and see Appendix D.6 for a proof):
Theorem 5 Letting δ′ = δw/(2K) and [v]+ := max{0, v} for all v ∈ R, we have
w(x) := [θ (n · q̂B(x);n, δ′)− E]+ θ (m · p̂B(x);m, δ′) + E ≤ w∗(x) ≤ w(x) := θ (n · q̂B(x), n, δ ′) + E [θ (m · p̂B(x),m, δ′)− E]+ (∀x ∈ X )
(17)
with probability at least 1− δw over SXm and TXn .
6We can use the same calibration set to construct IWs and prediction sets due to the union bound in Theorem 3.
We use these upper and lower bounds on the IWs as the inputs w and w to our algorithm—i.e.,
W = {w : X → R | ∀x ∈ X , w(x) ≤ w(x) ≤ w(x)} ,
and use b = maxx∈X w(x) = maxj∈[K] wj as the maximum IW. If b is known, we need importance weights associated with source calibration samples SXm ; thus we use a simpler form ofW as follows:
W = {w ∈ Rm | ∀i ∈ [m], wi ≤ wi ≤ wi} ,
where wi := w(xi) and wi := w(xi) for xi ∈ SXm . Theorem 5 under Assumption 1 is one way to construct uncertainty setW that contains the true IWs.
Corollary 1 Given K ∈ N, E ∈ R≥0, and δw ∈ (0, 1), suppose Assumption 1 is satisfied andW is constructed using Theorem 5. Then, we have
PSXm∼PmX ,TXn ∼QnX [w ∗ ∈ W] ≥ 1− δw.

heading: B.2 PARTITION CONSTRUCTION HEURISTIC
text: In general, exponentially many bins are needed to guarantee Assumption 1. Instead, we consider an intuitive heuristic for constructing these bins, so that the importance weights w(x)—rather than the density functions p(x) and q(x) individually—are roughly constant on each bin, inspired by (Park et al., 2021; 2020b); a standard heuristic for estimating IWs is to train a probabilistic classifier g(s | x) to distinguish source and target training examples, and then use these probabilities to construct the IWs. In particular, define the distribution
g∗(x, y) = 1 2 p(x) · 1(s = 1) + 1 2 q(x) · 1(s = 0).
Then, letting g∗(y | x) be the conditional distribution, we havew∗(x) = 1/g∗(s = 1 | x)−1 (Bickel et al., 2007). Thus, we train g(s | x) ≈ g∗(s | x) and construct bins according to w(x) = 1/g(s = 1 | x)− 1, i.e.,
Bj = {x ∈ X | w(x) ∈ [wj , wj+1)},
where 0 = w1 ≤ w2 ≤ ... ≤ wK+1 =∞. Finally, we describe how to train g. Let SXm′ and TXn′ be unlabeled training examples from a source and target, respectively; then, the set
RXm′,n′ = {(x, 1) | x ∈ SXm′} ∪ {(x, 0) | x ∈ TXn′ }
consists of i.i.d. samples (x, y) ∼ s∗(x, y). Thus, we can train s on RXm′,n′ using supervised learning. In practice, the corresponding IW estimates w(x) can be inaccurate partly since w is likely overfit to RXm′,n′ , which is why re-estimating the IWs in each bin according to Theorem 5 remains necessary.

heading: B.3 DENSITY ESTIMATION SATISFYING ASSUMPTION 1
text: We consider the following binning strategy that satisfies Assumption 1 for Lipschitz continuous PDFs. In particular, assume the PDFs p(x) and q(x) are L-Lipschitz continuous for some norm ‖·‖. Then, for any given E, construct each bin Bj such that
‖x− x′‖ ≤ E L · vj
∀x, x′ ∈ Bj ,
where vj is the volume of bin Bj . Then, we know that
|p(x)− p(x′)| ≤ L‖x− x′‖ ≤ E vj ,
so ∫ Bj |p(x)− p(x′)|dx′ ≤ E.
Similarly, we have ∫ Bj |q(x)− q(x′)|dx′ ≤ E.
Thus, the bins satisfy Assumption 1. Finally, assuming ‖ · ‖ is the L∞ norm, we describe one way to construct bins. In particular, construct bins by taking an ε-net with ε = (E/L)1/(d+1), where d is the dimension of X . Then, we have
‖x− x′‖ ≤ ε = ε d+1
εd =
E
L · vj ,
as desired. In general, the Lipschitz constant L for the L∞ norm may scale exponentially in d.

heading: C ADDITIONAL DISCUSSION ON IMPORTANCE WEIGHTS

heading: C.1 APPROXIMATE IMPORTANCE WEIGHTS
text: In Section 3.3, we assumeW has the following form: W := {w ∈ Rm | wi ≤ wi ≤ wi} .
However, considering the fact that the expected importance weight is one, i.e., Ex∼P [w∗(x)] = 1, the uncertainty setW that contains the true importance weights with high probability can be further constrained as follows:
W ′ := {w ∈ Rm | ∀i ∈ [m] , wi ≤ wi ≤ wi , c ≤ ∑m i=1 wi ≤ c}
for some c and c. In particular, using the Hoeffding’s inequality for example, we can estimate c and c such that w1, . . . , wm can be the part of the true importance weight w∗ that satisfies the mass constraint Ex∼P [w∗(x)] = 1.
Recall that we have to find a maximizer in the uncertainty set W ′ as in (8); however, due to the additional constraint on ∑m i=1 wi inW ′, it is challenging to solve the maximization problem exactly.

heading: C.2 MAXIMUM IMPORTANCE WEIGHT ESTIMATION
text: To generalize our approach to estimate the maximum importance weight b, we redefine the uncertainty setW over importance weights as follows:
W := {w : X → R | ∀x ∈ X , w(x) ≤ w(x) ≤ w(x)} for some w : X → R and w : X → R such that it contains the true importance weight w∗ with high probability—i.e.,
PSXm∼PmX ,TXn ∼QnX [w ∗ ∈ W] ≥ 1− δw. (18)
Given this, the maximum importance weight is obtained as follows:
b̂ = max x∈X w(x).
Considering that W can be estimated using binning as in Appendix B, the maximum importance weight is rewritten as follows:
b̂ = max x∈X w(x) = max j∈[K] wj , (19)
where
wj := θ (n · q̂j , n, δ′) + E
[θ (m · p̂j ,m, δ′)− E]+ .
Here, θ, θ, m,n, δ′, E, p̂j , and q̂j are defined in Theorem 5 and Appendix B.
The guarantee (18) implies the guarantee (7). Thus, for the maximization over the uncertainty set in (8), we use the original definitionW and the greedy algorithm for ŵ in (9). Then, Theorem 4 with the estimated b using (19) still holds as follows:
Theorem 6 Suppose Assumption 1 is satisfied and W is estimated using Theorem 5. Define ŵτ,SXm,TXn as in (9) and b̂SXm,TXn as in (19), making the dependency on τ , S X m , and T X n explicit, and URSCP as in (5). Let τ̂ be the solution of the following problem:
τ̂ = max τ∈T τ subj. to URSCP(Cτ , Sm, V, ŵτ,SXm,TXn , b̂SXm,TXn , δC) ≤ ε. (20)
Then, we have PSm∼Pm,V∼Um,TXn ∼QnX [LQ(Cτ ) ≤ ε] ≥ 1− δC − δw for any τ ≤ τ̂ .

heading: C.3 CHOOSING HYPERPARAMETERS
text: In general, there is no systematic way to choose the smoothness parameterE and the number of bins K; we briefly discuss strategies for doing so.
Number of binsK. The bins are defined in one dimensional space as described in Appendix B.2, so we follow the standard practice in the calibration literature for binning (Guo et al., 2017; Park et al., 2021), where K is between 10 and 20. As we are using equal-mass binning, we choose the number of bins so that each bin contains sufficiently many source examples (in our case, 5000 examples) for the length of the Clopper-Pearson interval over IWs of each bin to be below some threshold (in our case, 10−3), which leads us to K = 10. We provide a sensitivity analysis in Appendix G.6.
Smoothness parameter E. Estimating E (i.e., computing the integral of the difference of source and target probabilities) is intractable in general as it requires that we perform density estimation in a high-dimensional space (i.e., 2048 in our case), and then integrate this density over each bin. To avoid this hyperparameter selection, we choose E equal to zero or close to zero (in our case, 0.001). Intuitively, we find that binning based on the source-discriminator scores is an effective way to group examples with similar IWs; thus, the main contribution to the uncertainty is the error of the point estimate. We provide a sensitivity analysis on E in Appendix G.5. Importantly, note that PS-W even with E = 0 satisfies PAC criterion. One direction for future work is devising better strategies for choosing E.

heading: C.4 COMPARISON WITH WSCI
text: Guarantee. The main difference between WSCI and PS-W lies in the guarantee provided (rather than prediction set size); PS-W provides a stronger guarantee than WSCI. In particular, WSCI does not provide a PAC guarantee. Instead, to satisfy the coverage probability guarantee, it requires a new calibration set for every new test example. However, in practice, we usually have a single held-out calibration set. Thus, our approach PS-W provides a guarantee that holds conditioned on this set. This difference is illustrated in Figure 3a, which compares WSCI and PS-W given the true importance weight. PS-W produces a larger set size than WSCI, but strictly satisfies the error constraint. The shortcoming of WSCI can also be observed in Figure 2 in Tibshirani et al. (2019), which empirically shows that the guarantee only holds on average over examples.
Usage of IWs. WSCI requires the true IWs, whereas our method can use approximate IWs. Also, IWs are used differently. Given the true importance weights, WSCI uses the importance weights to reweight examples in the calibration set, and PS-W uses the importance weights to generate a target labeled calibration set using rejection sampling.

heading: D PROOFS

heading: D.1 PROOF OF THEOREM 1
text: First, note that the constraint in (1) implies F (L̄Sm(Cτ );m, ε) ≤ δ; conversely, any value of τ satisfying F (L̄Sm(Cτ );m, ε) ≤ δ also satisfies L̄Sm(Cτ ) ≤ k(m, ε, δ). Thus, we can rewrite (1) as
τ̂ = arg max τ∈T τ subj. to F (L̄Sm(Cτ );m, ε) ≤ δ.
On the other hand, if τ satisfies (3), then by definition of UCP(Cτ , Sm, δ) = θ(k;m, δ), we have
inf{θ ∈ [0, 1] | F (L̄Sm(Cτ );m, θ) ≤ δ} ∪ {1} ≤ ε,
which implies that F (L̄Sm(Cτ );m, ε) ≤ δ (since the infimum is obtained within the set since the binomial CDF F is continuous in ε). Conversely, any value of τ satisfying F (L̄Sm(Cτ );m, ε) ≤ δ also satisfies UCP(Cτ , Sm, δ) ≤ ε. Thus, we can rewrite (3) as
τ̂ = arg max τ∈T τ subj. to F (L̄Sm(Cτ );m, ε) ≤ δ,
implying (1) and (3) are equal. Thus, the claim follows from Theorem 1 of (Park et al., 2020a) and the fact that the error LP (Cτ ) is monotonically increasing in τ .

heading: D.2 MONOTONICITY OF THE CLOPPER-PEARSON BOUND
text: The CP bound UCP enjoys certain monotonicity properties that we will need. Intuitively, the CDF decreases as the number of observations m increases while holding the number of successes k fixed, but increases if both m and k are increased by the same amount (i.e., holding the number of failures m− k fixed). In particular, we have the following:
Lemma 2 We have θ(k;m− 1, δ) ≥ θ(k;m, δ) and θ(k − 1;m− 1, δ) ≤ θ(k;m, δ).
Proof. Recall that F (k;m, θ) is the cumulative distribution function of a binomial distribution Binom (m, θ), or equivalently of the random variable ∑m i=1Xi, where Xi ∼ Bernoulli(θ) are i.i.d.
Decreasing case. If k ≤ m− 1, then we have
m∑ i=1 Xi ≤ k ⇒ m−1∑ i=1 Xi ≤ k,
hence
P [ m∑ i=1 Xi ≤ k ] ⊆ P [ m−1∑ i=1 Xi ≤ k ] ,
so F (k;m, θ) ≤ F (k;m− 1, θ). Then, we have
θ(k;m, δ) := inf {θ ∈ [0, 1] | F (k;m, θ) ≤ δ} ∪ {1} ≤ inf {θ ∈ [0, 1] | F (k;m− 1, θ) ≤ δ} ∪ {1} =: θ(k;m− 1, δ),
thus θ is monotonically non-increasing in m.
Increasing case. We have m−1∑ i=1 Xi ≤ k − 1⇒ m∑ i=1 Xi ≤ k,
hence
P [ m−1∑ i=1 Xi ≤ k − 1 ] ⊆ P [ m∑ i=1 Xi ≤ k ] ,
so F (k − 1;m− 1, θ) ≤ F (k;m, θ). Then, we have
θ(k;m, δ) := inf {θ ∈ [0, 1] | F (k;m, θ) ≤ δ} ∪ {1} ≥ inf {θ ∈ [0, 1] | F (k − 1;m− 1, θ) ≤ δ} ∪ {1} =: θ(k − 1;m− 1, δ),
thus θ is monotonically jointly non-decreasing in (m, k).

heading: D.3 PROOF OF THEOREM 2
text: The rejection sampling prediction set consists of two steps: (i) generate target samples, using source samples Sm, importance weightsw, and an upper bound on their maximum value b, and (ii) construct the Clopper-Pearson prediction set using the generated target samples.
From rejection sampling, we choose N := ∑m i=1 σi samples from Sm, denoting them by TN ; here, N ∼ Binom (m, 1/b), and 1/b is the acceptance probability (von Neumann, 1951)—i.e.,
P [ V ′ ≤ w(X)
b
] = 1
b ,
where V ′ ∼ Uniform([0, 1]). The samples in TN are independent and identically distributed, conditionally on the random number N of samples being equal to any fixed value n. The reason is that one can view the rejection sampling algorithm proceeding in stages, iterating through the samples one by one. The first stage starts at the very beginning, and then each stage ends when a datapoint is accepted, followed by starting a new stage at the next datapoint. The last stage ends at the last datapoint.
Based only on the source samples observed in one stage, rejection sampling produces a sample from the target distribution. Thus, within each stage, we produce one sample from the target distribution, and because each stage is independent of all the other ones, conditionally on any number of stages reached, our produced target samples are iid. Thus, we can use the Clopper-Pearson bound conditionally on each N = n.
To this end, let τ̂(Sm, V ) = τ̂ to explicitly denote the dependence on Sm and V , and let
τ̃(Tn) = arg max τ∈T
τ subj. to URSCP(Cτ , Tn, δ) ≤ ε.
Note that conditioned on obtaining n samples using rejection sampling (i.e., |Tn(Sm, V, w, b)| = n), we have τ̂(Sm, V ) D = τ̃(Tn), where D = denotes equality in distribution. Then, we have
PSm∼Pm,V∼Um [ LQ(Cτ̂(Sm,V )) ≤ ε ] =
m∑ n=0 PSm∼Pm,V∼Um [LQ(Cτ̂(Sm,V )) ≤ ε | N = n] ·P[N = n]
= m∑ n=0 PTn∼Qn [LQ(Cτ̃(Tn)) ≤ ε] ·P[N = n]
≥ m∑ n=0 (1− δ) ·P[N = n]
= 1− δ, where the inequality follows by Theorem 1. The claim follows.

heading: D.4 PROOF OF THEOREM 3
text: First, let
τ̃ = arg max τ∈T
τ subj. to URSCP(Cτ , Sm, V, ~w∗, b, δC) ≤ ε, (21)
which satisfies PSm∼Pm,V∼Um [LQ(Cτ̃ ) ≤ ε] ≥ 1 − δC by Theorem 2. Now, with probability at least 1− δw, we have ~w∗ ∈ W . Under this event, we have
URSCP(Cτ , Sm, V, ~w ∗, b, δC) ≤ max
w∈W URSCP(Cτ , Sm, V, w, b, δC),
so τ̂ satisfies the constraint in (21). Thus, we must have τ̂ ≤ τ̃ . By monotonicity of LQ(Cτ ) in τ , we have LQ(Cτ̂ ) ≤ LQ(Cτ̃ ), which implies that
PSm∼Pm,V∼Um,TXn ∼QnX [LQ(Cτ̂ ) ≤ ε] ≥ PSm∼Pm,V∼Um [LQ(Cτ̃ ) ≤ ε] ≥ 1− δC ,
where the last step follows by Theorem 2. The claim follows by a union bound, since ~w∗ ∈ W with probability at least 1− δw.

heading: D.5 PROOF OF LEMMA 1
text: Let w and v be IWs where w(xi) ≥ v(xi) and w(xj) = v(xj) for j 6= i. Additionally, we use the following shorthands:
nw := m∑ i=1 1 ( Vi ≤ w(xi) b ) ,
Tnw := { (xi, yi) ∈ Sm ∣∣∣∣ Vi ≤ w(xi)b } ,
kw := ∑
(x,y)∈Tnw
1 (y /∈ C(x)) ,
nv := m∑ i=1 1 ( Vi ≤ v(xi) b ) ,
Tnv := { (xi, yi) ∈ Sm ∣∣∣∣ Vi ≤ v(xi)b } , and
kv := ∑
(x,y)∈Tnv
1 (y /∈ C(x)) .
Here, nw ≥ nv since w(xi) ≥ v(xi). Finally, recall that F (k;m, θ) be the cumulative distribution function of a binomial random variable ∑m i=1Xi, where Xi ∼ Bern(θ).
Non-decreasing case. If yi /∈ C(xi), there are two cases to consider:
1. If v(xi)b < Vi ≤ w(xi) b , then we can verify that nw = nv + 1 and kw = kv + 1.
2. Otherwise, we can verify that nw = nv and kw = kv .
In both cases, kw ≥ nv and nw ≥ nv . Since θ is monotonically jointly non-decreasing in (m, k) as in Lemma 2, we have
URSCP(C, Sm, V, w, b, δ) := UCP(L̄Tnw (C), δ)
:= θ(kw;nw, δ)
≥ θ(kv;nv, δ) =: UCP(L̄Tnv (C), δ)
=: URSCP(C, Sm, V, v, b, δ),
thus URSCP is monotonically non-decreasing in w(xi).
Non-increasing case. If yi ∈ C(xi), then kw = kv . Since θ is monotonically non-increasing in m as in Lemma 2, we have
URSCP(C, Sm, V, w, b, δ) := UCP(L̄Tnw (C), δ)
:= θ(kw;nw, δ)
≤ θ(kv;nv, δ) =: UCP(L̄Tnv (C), δ)
=: URSCP(C, Sm, V, v, b, δ),
thus URSCP is monotonically non-increasing in w(xi).

heading: D.6 PROOF OF THEOREM 5
text: Recall that
p̂B(x) := K∑ j=1 1 (x ∈ Bj)  1 m ∑ x′∈SXm 1 (x′ ∈ Bj)  , q̂B(x) :=
K∑ j=1 1 (x ∈ Bj)  1 n ∑ x′∈TXn 1 (xi ∈ Bj)  , pB(x) :=
K∑ j=1 1 (x ∈ Bj) ∫ Bj p(x′) dx′,
qB(x) := K∑ j=1 1 (x ∈ Bj) ∫ Bj q(x′) dx′, and
v(x) := vj(x) = ∫ Bj(x) dx′.
Due to the assumption of (12), |v(x) · p(x)− pK(x)| is bounded for any x ∈ Bj as follows:
|v(x) · p(x)− pB(x)| = ∣∣∣∣∣ ∫ Bj p(x) dx′ − ∫ Bj p(x′) dx′ ∣∣∣∣∣ = ∣∣∣∣∣ ∫ Bj p(x)− p(x′)dx′ ∣∣∣∣∣
≤ ∫ Bj |p(x)− p(x′)|dx′
= E. (22)
Similarly,
|v(x) · q(x)− qB(x)| ≤ E. (23)
Observe that mp̂(x) ∼ Binom ( m, ∫ Bj p(x′)dx′ ) for any x ∈ Bj ; thus pK is bounded with probability at least 1− δ′ as follows due to the Clopper-Pearson interval (θ, θ):
θ(mp̂(x);m, δ′) ≤ pK(x) ≤ θ(mp̂(x);m, δ′). (24)
Similarly,
θ(nq̂(x);n, δ′) ≤ qK(x) ≤ θ(nq̂(x);n, δ′). (25)
From (22), (23), (24), and (25), the following holds:
θ(mp̂(x);m, δ′)− E ≤ v(x) · p(x) ≤ θ(mp̂(x);m, δ′) + E and θ(nq̂(x);n, δ′)− E ≤ v(x) · q(x) ≤ θ(nq̂(x);n, δ′) + E.
Therefore, for any x ∈ Bj , w∗(x) is bounded as follows:
θ(nq̂(x);n, δ′)− E θ(mp̂(x);m, δ′) + E ≤ w∗(x) = q(x) p(x) ≤ θ(nq̂(x);n, δ ′) + E θ(mp̂(x);m, δ′)− E .
Since we apply the Clopper-Pearson interval for K partitions for both source and target, the claim holds due to the union bound.
Algorithm 2 PS: an algorithm using the CP bound in (3) procedure PS(Sm, f, T , ε, δ)
τ̂ ← 0 for τ ∈ T do (.) Grid search in ascending order
if UCP(Cτ , Sm, δ) ≤ ε then τ̂ ← max(τ̂ , τ) else break
return τ̂
Algorithm 3 PS-C: an algorithm using the CP bound in (3) with ε/b procedure PS-C(Sm, f, T , b, ε, δ)
return PS(Sm, f, T , ε/b, δ)

heading: E ADDITIONAL ALGORITHMS

heading: E.1 PS ALGORITHM

heading: E.2 PS-C ALGORITHM
text: We describe the PS-C algorithm, which uses a conservative upper bound on the CP interval. Let
LP (C) := E (x,y)∼P
[1 (y /∈ C(x))]
LQ(C) := E (x,y)∼Q
[1 (y /∈ C(x))]
w∗(x) := q(x)
p(x)
b := max x∈X
w∗(x).
Then, we have
LQ(C) = E (x,y)∼Q
[1 (y /∈ C(x))]
= E (x,y)∼P
[w∗(x)1 (y /∈ C(x))]
≤ E (x,y)∼P [b · 1 (y /∈ C(x))]
= b E (x,y)∼P
[1 (y /∈ C(x))]
= b · LP (C).
Thus, LQ(C) ≤ ε if b · LP (C) ≤ ε. Equivalently, LQ(C) ≤ ε if LP (C) ≤ ε/b. As a consequence, we can choose C based on the CP bound for the i.i.d. case (i.e., Algorithm 2), except using the desired error of ε/b (instead of ε). The algorithm is described in Algorithm 3.

heading: E.3 PS-R ALGORITHM

heading: E.4 PS-M ALGORITHM

heading: F EXPERIMENT DETAILS

heading: F.1 DOMAIN ADAPTATION
text: We use a fully-connected network (with two hidden layers, where each layers has 500 neurons followed by ReLU activations and a 0.5-dropout layer) as the domain classifier (recall that the input of this domain classifier is the last hidden layer of ResNet101). We use the last hidden layer of
Algorithm 4 PS-R: an algorithm using the RSCP bound in (6) procedure PS-R(Sm, f, T , w, b, ε, δ)
V ∼ Uniform([0, 1])m τ̂ ← 0 for τ ∈ T do (.) Grid search in ascending order
if URSCP(Cτ , Sm, V, w, b, δ) ≤ ε then τ̂ ← max(τ̂ , τ) else break
return τ̂
Algorithm 5 PS-M: an algorithm using the RSCP bound in (6) along with IWs rescaling
procedure PS-M(Sm, TXn , f, T , w, b, ε, δ) ŵ(x)← q̂B(x)p̂B(x) for x ∈ X , where p̂B and q̂B are defined in (14) and (15), respectively return PS-R(Sm, f, T , ŵ, b, ε, δ)
the model as example space X , where its dimension is 2048. For neural network training, we run stochastic gradient descent (SGD) for 100 epochs with an initial learning rate of 0.1, decaying it by half once every 20 epochs. The domain adaptation regularizer is gradually increased as in (Ganin et al., 2016). We use the same hyperparameters for all experiments.
F.2 IMAGENETC-13
To train a model using domain adaptation, due to the large size of the target training set, we subsample the target training set to be the same size as the source training set on for each random trials.

heading: G ADDITIONAL RESULTS

heading: G.1 SYNTHETIC RATE SHIFT BY TWO GAUSSIANS
text: We demonstrate the efficacy of the proposed approaches (i.e., PS-R with the known IWs and PS-W with the estimated IWs) using a synthetic dataset consisting of samples from two Gaussian distributions.
Dataset. We consider two Gaussian distributions N (µ,Σ) and N (µ,Σ′) over 2048-dimensional covariate space X . Here, µ = 0; Σ and Σ′ are diagonal where Σ1,1 = 52, Σi,i = 10−1, Σ′1,1 = 1, and Σ′i,i = 10
−1 for i ∈ {2, . . . , 2048}. We consider the “flat” Gaussian N (µ,Σ) as the source and the “tall” GaussianN (µ,Σ′) as the target. Intuitively, there is a rate shift from the source to the target—i.e., the target examples are a subset of the source, but occur with higher frequencies. We use the following labeling function: p(y | x) = σ(5x1), where σ is the sigmoid function. Finally, we generate 50,000 labeled examples for each training, calibration, and test.
Results. We consider two different setups: 1) the true IW is known, and 2) the true IW is unknown. In Figure 3a, we demonstrate the prediction set errors given the true IW. As expected PS-R satisfies the PAC guarantee—i.e., the error is below ε. However, as shown in Figure 3b, when we need to estimate IWs, using just the point-estimate of the IW results in PS-R performing poorly in terms of prediction set error; it still satisfies the ε constraint, but the error is close to zero, indicating that the prediction set size is too large to be useful for an uncertainty quantifier. In contrast, PS-W (i.e., rejection sampling based on interval estimates of IWs) produces a larger, more reasonable error rate while still satisfying the PAC condition. These experiments demonstrate that PS-R works well when given the true IW, but accounting for IW uncertainty is important when using estimated IWs.

heading: G.2 PREDICTION SET SIZE AND ERROR

heading: G.3 ABLATION STUDY ON CALIBRATION
text: G.4 COMPARISON WITH VARIOUS TOP-K PREDICTION SETS
G.5 VARYING A SMOOTHNESS PARAMETER
G.6 VARYING A NUMBER OF BINS

heading: G.7 PREDICTION SET VISUALIZATION

Review summaries:
Let $X$ be an instance space, $Y$ a set of labels, $D$ some underlying (hidden) distribution over $X \times Y$. This work studies a new method for converting the output of a probabilistic predictor (e.g. a deep net) to a good prediction set: that is a mapping from $C: X \to 2^Y$, such that for most samples $(x,y) \sim D$, $C(x)$ contains $y$. Formally, the authors study this problem in the PAC setting with covariate shift. The learner is given access to the score function, a labeled sample from the source distribution $P$ over $X \times Y$, and unlabeled samples from a target distribution $Q$ over $X \times Y$ whose marginal over $X$ may be shifted from the source. The goal is to output a prediction set which is as small as possible while still retaining PAC guarantees with respect to the shifted target distribution. This is a well-motivated model in practice. Probabilistic outputs of modern neural nets can be difficult to interpret, and small prediction sets may be useful in settings where one wishes to avoid or check more carefully a few marked outputs (e.g. a problematic medical diagnosis).  Assuming known, bounded importance weights, the authors provide an algorithm satisfying the PAC guarantee by maximizing a cutoff value for the score function that performs well over an empirical sample. They ensure their guarantee holds over the target rather than source distribution by rejection sampling to simulate the target distribution. The authors also extend this to settings where the importance weights are unknown but can be estimated from unlabeled samples. Finally, the authors provide experimental evidence over a couple common settings of covariate shift. They show that their algorithm outperforms baseline methods in the literature in the sense that it maintains PAC guarantees while outputting a smaller prediction set  in expectation. 

This paper is concerned with learning prediction sets (in a PAC sense) under the covariate shift assumption, given a model $f(x,y)$. The form of the sets is restricted to $C_\tau(x) := \{y : f(x,y) \ge \tau\},$ and the problem is set up as learning a $\tau$, using a  labelled source dataset $S_m \sim P^{\otimes n}$ and an unlabelled target dataset $T_n \sim Q^{\otimes n}$ such that with high probability over $S_m, T_n,$ $Q(Y \in C_\tau(X)) \ge 1-\varepsilon,$ where $\varepsilon$ is a given target coverage level. It is desired that $\tau$ is as large as possible to minimise the size of the prediction sets learned.  The main scheme is presented modularly. The paper first describes how finding the maximum $\tau$ whilst ensuring that the number of points in $S_m$ it captures is large enough (as specified by using a Binomial tail inverse) gives valid sets when $Q = P$. Next, it is argued that when $Q$ is absolutely continuous with respect to $P$, and the derivative $w(x) = \frac{\mathrm{d}Q}{\mathrm{d}P}(x)$ is upper bounded and known, then one can importance sample the set $S_m$ to generate a sample from $Q$, which can then plug into the previous procedure. This step fundamentally uses the covariate shift assumption, and, to my understanding, is folklore. The following, then, constitute the main technical contributions.  Next, the assumption of exact knowledge of $w$ is relaxed, and it is argued that if instead for each $x_i \in S_m,$ bounds $\underline{w}_i \le w(x_i) \le \overline{w}_i$ are available, then one can produce a worst-case estimate of the coverage of $C_\tau$ for any $\tau$ (by taking points that were missed to have high weights, and points covered to have low weights), and if this pessimistic coverage also satisfies the constraint. So, for each $\tau$, the worst $w_\tau$ can be produced, which can then feed into the importance sampling procedure above.  Then, it is pointed out that such a confidence bound on the $w_i$s can be learned using a probabilistic classifier $s(\cdot|x)$ to separate data from $S_m$ and $T_n$, which leads to the main proposed algorithm. While it is roughly justified in the appendix that an accurate estimate can be obtained under smoothness assumptions with an appropriately fine gridding of the space (and a consequently huge sample complexity), the concrete proposal is to replace this step with a heuristic method for fitting bounds on the weights, and then plugging these into the above strategy. Note that this is not implemented as a direct optimisation over $\tau$ - instead a set $\mathcal{T}$ of possible values is pre-selected, and the procedure is executed for each $\tau$ (the algorithm recommends an increasing order on the same).  Finally, the paper presents experiments in the DomainNet dataset, and in the ImageNet dataset, where the shift in the latter corresponds to adversarial perturbations. The principal baselines are the weighted split conformal inference (WSCI) method, which in my opinion is an appropriate choice, and the "PS-C" method, which simply uses a prediction set that has $1-\varepsilon/b$ coverage on the source data, where $b$ is an estimated upper bound on $w$. The proposed method is ablatively presented, and it is seen that on just the DomainNet data, the method PS-R (which simply takes $(1-s)/s$ as an estimate of $w$) is both reliable and performs well, while the PS-M method, which further integrates samples in bins tends to be slightly optimistic on this dataset. Conversely, in the adversarially perturbed dataset, PS-R produces trivial prediction sets (since presumably this shift goes entirely outside the source domain's support), while PS-M performs well. The proposed method, PS-W, which further uses upper and lower bounds on $w$s after integrating them over bins, tends to be pessimistic (error-rates of $0.06-0.07$ are observed when only $0.1$ is demanded), but is not too much worse than either of these methods, and is effective in both types of shifts.

The paper presents algorithms for PAC prediction sets under the assumption of covariate shift. Using estimated/known importance weights that encode the shift, the method optimizes for the smallest subset of labels such that with high probability the error of this prediction set is low. A rejection sampling based strategy is then shown to satisfy the PAC constraints. The methodology is extended to the case where the importance weights are uncertain (due to estimation error). It is shown that the robust variant can be solved using the extreme case weights and the rejection sampling based algorithm. Simulations show the efficacy of the method.

This paper proposes a new method to construct approximately correct (PAC) prediction sets for uncertainty quantification in the presence of covariate shift. It is a natural and interesting extension of the previous works [Park et al. 2020a, 2021] on PAC prediction sets. The building blocks this paper took from these previous works are the optimization problem in Equation (1) and the Clopper-Pearson confidence intervals for the Binomial distribution in Section 2.2. The extension is based on a rejection sampling Clopper-Pearson bound given in Section 3.2.  The authors propose an algorithm with and without access to the true important weights. The algorithm is evaluated in the settings of “rate shift” and “support shift” on the DomainNet and ImageNet datasets. The experiments show that the algorithm gives the smallest prediction sets among approaches that always satisfy the PAC constraint.