3.4 Conformal prediction

3.4.1 Fundamentals

Suppose \(X_1, \ldots , X_n\) are iid samples from some probability distribution \(P\), and \(Y_1, \ldots , Y_m\) are iid samples from some other probability distribution \(Q\). Our goal is to test the hypothesis that \(P = Q\) even if we don’t know either of the distributions. The strategy we learn in early statistics is that we should choose a test statistic \(T\) (for example we can let \(T = |\overline {X} - \overline {Y}|\)) and reject the null based on whether \(T\) is unusually large. The question we need to ask is “how large does \(T\) need to be,” and typically this is where permutation tests are introduced: we compare \(T\) to how the statistic would look if we were to permute the data.

Formally, we introduce the following (randomized) permutation distribution: choose \(M\) uniform permutations \(\sigma _1, \ldots , \sigma _M\) from \(S_{n+m}\) (the set of permutations acting on a list of \(n + m\) objects), and we will have these permutations act on the vector \[ Z = (X_1, \ldots , X_n, Y_1, \ldots , Y_m). \] Specifically, the \(\sigma _i\)s will shuffle the entries of \(Z\), and then we evaluate the test statistic on the permuted entry. So we always take the average of the first \(n\) entries and the average of the last \(m\) entries after permuting, and we find their absolute difference; in other words, we compare \[ T(Z) = \left |\overline {Z_{1:n}} - \overline {Z_{n+1:n+m}}\right | \] to the corresponding values of \(T(Z_\sigma )\), and the \(p\)-value will essentially be the relative rank of \(T(Z)\) compared to \(T(Z_{\sigma _1})\) through \(T(Z_{\sigma _M})\): \[ p = \frac {1 + \sum _{i=1}^{M} \mathbf {1}\{T(Z_{\sigma _i}) \geq T(Z)\}}{1 + M}. \] Under the null, this is either uniformly distributed on \(\left \{\frac {1}{M+1}, \frac {2}{M+1}, \ldots , \frac {M+1}{M+1}\right \}\) or biased upward due to ties (because under the null we have exchangeability of the vector \(Z\)). Thus it is indeed a \(p\)-value, and notice that this doesn’t rely on us needing to know the distribution of \(T\) at all. Given a set of permutations \(S\), consider the quantity \[ p = \frac {1 + \sum _{\sigma \in S} \mathbf {1}\{T(Z_\sigma ) \geq T(Z)\}}{1 + |S|}. \] If the variables are exchangeable, then \(Z_\sigma \) has the same distribution as \(Z\) and thus \(T(Z_\sigma )\) has the same distribution as \(T(Z)\). We are interested in what sets of permutations \(S\) yield a valid \(p\)-value. Here are the answers:

Note that this is not a confidence interval – this is an observation we have not seen yet, and it’s interesting that we can solve this problem at all. In fact, the way we can do so is through permutation tests! We hypothesize a value of \(Y_{n+1}\) and test for exchangeability by computing \[ p^y = \frac {1}{(n+1)!} \sum _{\sigma \in S_{n+1}} \mathbf {1}\{T(Z^y_\sigma ) \geq T(Z^y)\}, \] where \(Z^y = \{(X_1, Y_1), \ldots , (X_n, Y_n), (X_{n+1}, y)\}\). Then setting \(y \in C(X_{n+1})\) if and only if \(p^y \geq \alpha \), we claim this is a valid prediction interval regardless of our choice of \(T\). Indeed, \(Y_{n+1}\) itself will be in the confidence interval if and only if \(p^{Y_{n+1}} \geq \alpha \), but we’ve already shown that \(p^{Y_{n+1}} \stackrel {\text {sto}}{\geq } U\) and thus we are done.

The motivation for conformal prediction is that we want some uncertainty in our prediction and some way of quantifying accuracy.

3.4.2 Approaches

Full conformal prediction is typically done as follows: we fit a model \(\hat {\mu }(\cdot ) = \mathcal {A}(Z^y)\) to \(Z^y = ((X_1, Y_1), \ldots , (X_n, Y_n), (X_{n+1}, y))\) which satisfies the symmetry assumption (for example in random forests, it doesn’t matter what order we pass in the data, and this is true of most algorithms), and we define \(T(Z^y) = |y - \hat {\mu }(X_{n+1})|\) to be the residual. Replacing \(T\) with what we have, we thus get a \(p\)-value \[ p^y = \frac {1}{n+1} \sum _{i=1}^{n+1} \mathbf {1}\{|Y_i - \hat {\mu }(X_i)| \geq |y - \hat {\mu }(X_{n+1})|\} = \frac {1}{n+1} \sum _{i=1}^{n+1} \mathbf {1}\{R^y_i \geq R^y_{n+1}\}. \] and again we include \(y\) in the interval if and only if \(p^y\ge \alpha \). This is all computationally intensive – every time we fit the model we might need to do a lot of computation, and in fact the prediction interval can be a disjoint set of intervals instead.

We introduce split conformal prediction then. This is a special case where we have \(n\) data points and we do sample splitting: we learn a model \(\hat {\mu }\) with the first split (also called a “fold”), and on the second split we calculate out-of-sample residuals (that is, learn the distribution of the residuals \(R_i = |Y_i - \hat {\mu }(X_i)|\)). Then the test residual relates to this second split by keeping track of quantiles, and the point is that we separately do training and calibration and form our interval from points that have all not been used for training.

Formally, we compute a score function \(S(x, y) = |y - \hat {\mu }(x)|\) by fitting a model on an independent training set. Once we have this, we use a distinct calibration set of size \(n\) to find typical size of residuals \[ S_i = S(X_i, Y_i), \quad S^y_{n+1} = S(X_{n+1}, y). \] The point is that if \(y = Y_{n+1}\) these points should all be indistinguishable (they’re from the same distribution), and now we include \(y\) if \[ p^y = \frac {1}{n+1} \sum _{i=1}^{n} \mathbf {1}\{S_i \geq S^y_{n+1}\} \geq \alpha . \] And this result now holds conditionally on the training set, since for all purposes \(\hat {\mu }\) is fixed.