Hello,

Thank you for providing these notebooks for conformal prediction, they have been immensely helpful.

Reading through the section 2.1 of the paper on "Classification with Adaptive Prediction Sets" and the associated notebook, I had some questions about the scoring function.

Namely, the paper provides the score function

$$s(x,y) = \sum_{j=1}^k \hat{f}(x)_{\pi_j(x)}$$

where $y = \pi_j(x)$. Why are we including $\hat{f}(x)_{y}$ in the sum? Doing so would lead to some possibly problematic scores. Consider, for example, a perfect predictor that assigns all of its mass to the correct label $y$, and a completely incorrect predictor that assigns all of its mass to some incorrect label $\ne y$. Both of these predictors would have the same score of 1. This breaks the assumption that a higher score corresponds to misalignment between the forecaster and the true label.

Investigating this issue further, I tried modifying the score function to greedily include all classes up to, but not including, the true label. Intuitively, a higher score would correspond to more probability mass assigned to incorrect labels, which is a better estimate of misalignment. Coding this up in the notebook for APS, this little fix increased the coverage slightly, but more importably it decreased the mean size of the confidence sets to 3.3 (compared to 187.5 in the original notebook). The confidence sets on the imagenet examples also seem to make more sense upon preliminary inspection. This could possibly address an issue raised previously.

Is there a typo/error in the score function of APS that would explain these results?

Thanks in advance!

Thanks a lot for writing this very interesting paper and for publishing the notebooks.

FYI, the interpolation argument appears to be deprecated in recent versions of numpy.quantile:

https://numpy.org/devdocs/reference/generated/numpy.quantile.html

Requires an external file. Maybe change to load it remotely from the repo so colab just works?

First of all, I would like to thank all contributors to this repository. Appreciate the great work that goes behind creating and maintaining this repository.

I was looking through the notebook, weather-time-series-distribution-shift.ipynb, and notice that in the last 4 lines of the second section, we have:
sort_idx = np.argsort(times)
pred_mean = pred_mean[sort_idx]
temperatures = temperatures[sort_idx]
times = times[sort_idx]

Should the sorting also be done for the uncertainty data, i.e., by adding the line pred_uncertainty=pred_uncertainty[sort_idx]?

Thank you for providing such a valuable resource. I have a few inquiries regarding the APS-Randomized algorithm.
To begin, I'd like to refer to the upper bound result for CP calibration, as stated in Theorem 2.2 of "Distribution-Free Predictive Inference For Regression":
$\mathbb{P}(Y_{test} \in C(X_{test}, U_{test}, \hat{q})) < 1 - \alpha + \frac{1}{n-1}$

Upon running the APS-Randomized algorithm for 100 trials, I observed a mean coverage of approximately 93%, consistent with the empirical coverage in the provided repo example (0.93020408163265). I can raise a rationale for this deviation: while the "split conformal algorithm" in the referenced paper operates with a deterministic model ($\mathcal{A}$), while in APS-Randomized, both the generated scores and the threshold are randomized, which may cause potential challenges.

Moreover, apart from the favorable over-coverage exhibited by this algorithm, its conditional coverage, quantifiable using metrics like SSCV, surpasses that of the APS algorithm outlined in the RAPS paper (which is RAPS with $\lambda = 0$), while maintaining identical set sizes. I'm interested in understanding the underlying rationale behind this algorithm and would appreciate insights into its origins, particularly if it was derived from a specific academic paper.
Thank you for your assistance!
Lahav.

Hi, thank you so much for the great work. I have a question regarding the notebook on conformal risk control.

in the notebook, you defined the risk optimization objective as

def lamhat_threshold(lam): return false_negative_rate(cal_sgmd>=lam, cal_gt_masks) - ((n+1)/n*alpha - 1/(n+1))

However, in section 4.3 of the paper, it is defined as $$\alpha - \frac{B-\alpha}{n}$$, which means the denominator should be $n$. Why is it $n+1$ in the code? Thanks.

While trying to understand the method for selective classification, I tried to run your code and to plot the curves of the selective risk and its upper bound.
This is the code that I added to plot these curves:

selective_risk_values = np.array([selective_risk(lam) for lam in lambdas])
selective_risk_values_ub = np.array([selective_risk_ub(lam) for lam in lambdas])
plt.plot(lambdas, selective_risk_values, label='Selective Risk')
plt.plot(lambdas, selective_risk_values_ub, label='Selective Risk upper bound')
plt.axhline(y=alpha, color='red', linestyle='--', label='alpha')
plt.legend()

I was expecting to see the upper bound to be monotonic and descending, but as you can see in the image below, it is not.

From the paper "Gentle introduction to conformal prediction", I assumed that the upper bound should be monotonic because it was introduced to overcome the issue related to the fact that the selective risk is non monotonic (section 5.5 Selective classification).

When testing it with my own dataset, I get an extreme example of this behaviour.

Is this an expected behaviour or the assumption of the upper bound being monotonic is wrong?

Thank you!

First many thanks for this awesome repo and the tutorial on split conformal prediction. While going through the conformal prediction under distribution drift section and the corresponding example weather prediction with time-series distribution shift, I noticed that the naive implementation for determining $\hat{q}$ uses expanding window that takes all scores up to time $t$ and compute the quantile and iterates over $t$

naive_qhats = np.array( [np.quantile(scores[:t], np.ceil((t+1)*(1-alpha))/t, interpolation='higher') for t in range(K+1, scores.shape[0]) ] )

But one can think of another approach by using a rolling window of fixed window size $K$ (in the example you were using $K=1000$) and then compute the quantile on each of the windows -- I rewrote and tested the function to support both options below

def compute_qhats(scores, alpha, wsize, opt='fixed_window'):
    lst = []
    qlst = []
    K = wsize
    for t in range(K+1, scores.shape[0]):
        if opt == 'fixed_window':
            start = t - K - 1
            nsamples = K + 1 # = t - start = t - (t - wsize -1) = wsize + 1 = K + 1    QED
        elif opt == 'expanding_window':
            start = 0
            nsamples = t
        q_level = np.ceil((nsamples+1)*(1-alpha))/nsamples
        qlst.append(q_level)
        tmp = np.quantile(scores[start:t], q_level, interpolation='higher')
        lst.append(tmp)
        print('start:', start, 'end:', t)
    return np.array(qlst), np.array(lst)

Plot of the $\hat{q}_{expandingwindow}$ expanding window approach (i.e. naive implementation)

vs. the plot for e $\hat{q}_{rollingwindow}$ rolling window approach

If we compare the results to the weighted conformal prediction approach

it is very similar to the rolling window but with the cost of additional computation for finding infimum of q
$$\hat{q} = inf \{ q : \sum^{n}_{i=1} \tilde{w}_i \mathbb{1} \{s_i \leq q\}\geq 1- \alpha \}$$

that requires finding the roots of the expression above after moving the $1-\alpha$ to the left side of the inequality (i am using the generalized expression but in practice we are using the adaptation to window based version in section 5.3).

Lastly here is the comparison of coverage over time of the three approaches

In fact when computing the overall coverage, the rolling window version achieves the best coverage with score 0.900665 vs. 0.8995545 for the weighted version.

It might be the case that for this data adding the constraint of finding the infimum does not add much, but overall if the argument of weighted conformal prediction is based on weighting the recent observation in a window then a sensibly defined window would be sufficient to counter the drift especially that we are not "learning the weights $w$" but rather fixing them to a uniform across the window in both cases (unless I am missing something here :) )

On a separate note, one minor issue in the code is the size of the window K. The way it is coded now it is translated to be K+1 observations used, and in the case of weighted conformal prediction, when computing qhats you omit the first observation from the computation.

Thank you again for your work and effort to make conformal prediction accessible to the masses.

Hello,

Thank you so much for providing the conformal prediction tutorial & corresponding notebooks, they are super helpful!

I had a question regarding the size of the prediction sets returned using the APS methods. In the implementation provided in the notebooks, the prediction sets are far larger than reported on your paper than introduced RAPS. The notebook implementation returns sets that are on average >200 labels, whereas the paper reports an average set size of 10.4, on ResNet152.

I have not done extensive evaluation on RAPS, but it seems the notebook implementation also returns slightly larger sets (set size of ~3).

I was wondering if you have any ideas as to what might be causing this discrepancy, and what the best way to replicate the results in the paper might be.

Also, I wasn't sure which repo this issue should be opened in, so apologies if it doesn't fit here. Thanks in advance!