Year: 2021

Abstract:
Accurate uncertainty quantification is a key building block of trustworthy machine learning systems. Uncertainty is typically represented by probability distributions over the possible outcomes, and these probabilities should be calibrated, e.g. the 90% credible interval should contain the true outcome 90% of the times. In the online prediction setup, existing conformal methods can provably achieve calibration assuming no distribution shift; however, the assumption is difficult to verify, and unlikely to hold in many applications such as time series prediction. Inspired by control theory, we propose a prediction algorithm that guarantees calibration even under distribution shift, and achieves strong performance on metrics such as sharpness and proper scores. We compare our method with baselines on 19 time-series and regression datasets, and our method achieves approximately 2x reduction in calibration error, comparable sharpness, and improved downstream decision utility.

heading: 1 INTRODUCTION
text: Accurate uncertainty quantification is crucial for machine learning predictions used in high-stakes decision making. Typically, uncertainty is represented by probability distributions over the possible outcomes, and these probabilities should be calibrated. In the regression setup, for example, the true label should be below the predicted 95% quantile for 95% of the samples (Gneiting et al., 2007). Calibration can convey confidence to decision makers because extreme values (e.g. true label above 95% quantile) are guaranteed to be rare. In addition to calibration, these probabilities should be sharp (i.e. concentrated and have low variance). Sharp probabilities are useful to decision makers because they are informative.
If data is i.i.d., then recalibration (Kuleshov et al., 2018) and conformal prediction (Vovk et al., 2020) algorithms can achieve low calibration error and good sharpness in the regression setup. However, the i.i.d. assumption is unlikely to hold in most time-series prediction tasks. Under distribution drift, it is possible to adapt regret minimization algorithms (Cesa-Bianchi & Lugosi, 2006; Kuleshov & Ermon, 2017) to achieve calibration. However, regret minimization calibration algorithms are designed for the asymptotic regime, and empirically we show that they are effective only with large sample size, and have very poor calibration and sharpness with short time series (e.g. 50 samples), limiting their practical utility.
Our goal is to design an algorithm to achieve good calibration and sharpness in the online regression setup even for short time series. We start with an existing prediction algorithm that has good calibration and sharpness for i.i.d. data (such as conformal prediction), and track the empirical frequency that the labels are below e.g. the 75% quantile. If the empirical frequency is significantly below 75%, we move future predictions up so more labels are below the prediction; vice versa. We design a very efficient adjustment method that can guarantee near perfect calibration with only tens of samples. We call this the “basic” prediction algorithm shown in Figure 1 (upper right).
Our main technical contribution is an improved algorithm that addresses two short-comings of the “basic” algorithm without breaking guaranteed calibration. First, the predictions should be feasible, e.g. the 75% quantile should never be smaller than the 25% quantile. The basic algorithm might violate this requirement because we adjust each quantile separately with no constraint on their relation. Second, the predictions should be stable if the distribution drift stops; empirically instability harms sharpness and proper scores. Our improved algorithm is based on an analogy between the
prediction task and a mechanical system. We encode the feasibility constraints as the feasible states of a mechanical system, and use proportional-integral-derivative (PID) control (Minorsky, 1922) to stabilize the system. Our final algorithm makes predictions that satisfy all three desirable properties: provable calibration, feasibility, and stability (Figure 1 bottom right).
We test our algorithm on two real world time series datasets (stock earnings and COVID cases) and a large benchmark of 17 common regression datasets. We reduce calibration error by approximately 2x compared too all baselines, while achieving comparable or better sharpness. We also simulate COVID response decisions based on predictions on the COVID dataset, and observe improved decision loss compared to baselines during periods of distribution drift (e.g. when case numbers surge).

heading: 2 PROBLEM SETUP AND BACKGROUND
text: We consider regression problems. Let X ∈ X denote the input feature and Y denote the label. We assume that the label is bounded by B, i.e. Y cannot take any value outside [−B,B] for some pre-specified B > 0.
We consider quantile predictions with a set of K equally spaced quantiles 0 < α1 < · · · < αK < 1. For example, when K = 9 we choose α1 = 10%, α2 = 20%, · · · , αK = 90%. A quantile prediction is a vector of K numbers (denoted by Z) where Zk should ideally predict the αk-th quantile of the label Y |X conditioned on the input feature X . We say a quantile prediction is feasible, if Yk < Yk+1,∀k, i.e. a lower quantile is always smaller than a higher quantile. In the limit of infinitely many quantiles K → ∞, a quantile prediction is equivalent to a probability distribution, however, for technical reasons that will become evident later, we only consider finitely many quantiles in this paper.
Our setup is an online prediction setup, where the forecaster sequentially makes the predictions. Consider a COVID prediction example: at time t, the forecaster observes some features Xt (such as current vaccination rate) and makes a prediction Zt about next week’s case number. After the forecaster makes a prediction, the true case number Y t is revealed; then we move on to time step
t + 1 and repeat this process. Following standard terminology we will use “nature” to refer to the process that generates the featuresXt and the true label Y t. Our setup is formalized by the following interaction between forecaster and nature: for t = 1, 2, · · · , T
1. Nature reveals feature Xt
2. Forecaster makes a quantile prediction Zt ∈ RK where Ztk is the αk-th quantile. 3. Nature reveals label Y t ∈ [−B,B]
We call a finite sequence of interactions a transcript, i.e. a sequence X1, Z1, Y 1, · · · , XT , ZT , Y T . A forecasting algorithm ψ is a function that maps a transcript and a new feature to a prediction:
ψ : X1, Z1, Y 1, · · · , Xt−1, Zt−1, Y t−1, Xt 7→ Zt. (1)
We will consider two types of assumptions on nature. We say that nature is i.i.d. if ∀t, Xt, Y t are drawn from some i.i.d. random variable (but we make no assumption on the probability law on this random variable). We say that nature is causal if nature chooses Xt, Yt without depending on the future. Specifically, Xt, Y t can only depend on the variables that precede it in the transcript (i.e. Xt can depend on X1, · · · , Y t−1 and Y t can only depend on X1, · · · , Y t−1, Xt, Zt). Causal is an extremely weak assumption. For example, nature can even adversarially choose the label Y t to increase prediction error after observing the forecaster’s prediction Zt.

heading: 2.1 CALIBRATION AND ACHIEVING CALIBRATION WITH CONFORMAL METHODS
text: A natural property that we can request is (probabilistic) calibration (Gneiting et al., 2007; Kuleshov et al., 2018). Intuitively, ∀k, the label Y t should be below the predicted αk-th quantile for an αk proportion of the samples. Formally, given a transcript of any length, consider the empirical frequency that the label Y t is below the predicted quantile Ztk up to time T (for any T that’s less than the length of transcript)
FTk = ∑
1≤t≤T
I(Y t ≤ Ztk). (2)
Ideally, the label Y t should be below the αk-th quantile exactly αkT many times, i.e., ∀k, FTk /T = αk. Achieving this perfectly is difficult, so we define an approximation. For any function b : N → R+ we say that a transcript is b-calibrated if ∀k = 1, · · · ,K,∀T , |FTk /T − αk| ≤ b(T ). Correspondingly, a forecasting strategy is calibrated if the resulting transcript is calibrated:, Definition 1. For any function b : N → R+, a forecasting algorithm ψ is b-calibrated under i.i.d. (or causal) assumptions if for any nature’s strategy that is i.i.d. (or causal), the resulting transcript is b-calibrated almost surely.
Conformal Calibration Conformal calibration (Vovk et al., 2020) is a forecasting algorithm based on conformal prediction ideas (Vovk et al., 2005; Shafer & Vovk, 2008). It is a wrapper algorithm that transforms a initial predictor (such as an off-the-shelf predictor) into new predictions that are provably calibrated under i.i.d. assumptions. Conformal calibration is based on the following intuition: for example, if 75% of past labels are below the initial prediction’s mean, then we can use the initial prediction’s mean as the 75%-quantile prediction in the future. More generally, we identify a statistic of the initial prediction, such that αk-proportion of the past labels are below that statistic. We use this statistic as the αk-th quantile of the future prediction.
If the data is i.i.d., conformal calibration is a very useful algorithm. This is because it has extremely strong calibration properties (Proven in Proposition 1 of (Vovk et al., 2020)), where the random variables I(Y t ≤ Ztk), t = 1, 2, · · · is a sequence of i.i.d. Bernoulli random variables with mean αk. This is impressive because even if we have access to the oracle forecaster (i.e. Y tk is equal to the true conditional αk-th quantile of Y t|Xt), we cannot achieve better calibration — I(Y t ≤ Ztk) is also a sequence of i.i.d. Bernoulli random variables with mean αk. Empirically, conformal prediction also has very good sharpness if the initial prediction function is reasonable (Burnaev & Vovk, 2014).
A minor property (but needed for proofs) is that the predictions should bounded in [−B,B]. Conformal calibration satisfies this property when t > K, i.e. there are more samples than predicted quantiles. To bypass this limitation, we can initialize the conformal calibration algorithm with at least K offline samples. Alternatively, if we know B in advance, we can clip the prediction by B.

heading: 3 CALIBRATION UNDER DISTRIBUTION DRIFT
text: This section introduces our algorithm that guarantees calibration even when nature is not i.i.d. We start from a basic algorithm in Section 3.1, and discuss its shortcomings. Then we move on to more advanced algorithms in Section 3.2 and 3.3 to address the shortcomings.

heading: 3.1 A BASIC CALIBRATED PREDICTION ALGORITHM
text: Conformal calibration achieves good calibration and sharpness when nature is i.i.d. When we are unsure if nature is i.i.d., our basic idea is to use conformal calibration until it has failed to achieve calibration. Specifically, when nature is i.i.d., I(Y t ≤ Ztk) is a sequence of Bernoulli random variables, so FTk is the sum of these Bernoulli random variables (i.e. a binomial). We use b ∗ δ(T ) to denote the confidence interval of a binomial distribution, i.e. Pr[FTk /T ∈ αk ± b∗δ(T )] = 1 − δ. We fix some small δ (such as 0.05) and if under the conformal calibration algorithm FTk /T 6∈ αk ± b∗δ(T ), we know (with 0.95 confidence) that nature is not i.i.d, so we will make adjustments to salvage calibration.
From a high level, F tk > αkt + b ∗ δ(t)t implies that the true label is below the αk-th quantile too often, so our prediction Ztk has been too large and we should reduce it; vice versa; There are some design freedom in choosing how much to reduce or increase the prediction. We choose an adjustment that grows exponentially with larger calibration error. The reason for this choice is to make very large adjustments when the calibration error is large, so we can tightly bound the calibration error. Formally, let Z̃tk denote the initial prediction generated by the conformal calibration prediction algorithm, the prediction of our algorithm Y tk is given by
Zt+1k = Z̃ t+1 k + E t+1 k , where E t+1 k =  1− e β(F tk−αkt−b ∗ δ(t)t) F tk > αkt+ b ∗ δ(t)t eβ(αkt−b ∗ δ(t)t−F t k) − 1 F tk < αkt− b∗δ(t)t
0 otherwise (3)
Intuitively, if F tk is within the correct range we make no adjustments (i.e. Z t+1 k = Z̃ t+1 k ); when F tk falls outside the correct range, we make aggressive adjustments that increase exponentially with how much it falls outside the correct range. β > 0 is a hyper-parameter that controls how large the adjustments are (we will show in Section 3.4 that these hyper-parameters are easy to choose). Note that Et+1k is for the (t + 1)-th time step (instead of t). This is because any adjustments can only depend on past information, so we will only have access to F tk at the (t+ 1)-th time step. The following theorem shows that such adjustments indeed guarantee calibration.
Theorem 1. Let Z̃tk be generated by any forecasting algorithm and bounded in [−B,B], for any δ > 0, the prediction algorithm defined by Eq.(3) is b-calibrated where b(T ) = b∗δ(T )+ log(2B+1)+β βT .
Note that the guarantee of Theorem 1 does not require that conformal calibration generate the initial prediction Z̃tk. The use of conformal calibration is motivated by its strong empirical performance (Vovk et al., 2020) when the data is i.i.d. or close to i.i.d. (which is further supported by our experiments).
The calibration guarantee is strong because under conformal calibration with i.i.d. data (or the oracle forecaster) we can expect to see b∗δ(T ) calibration error (with 1 − δ probability); without i.i.d. assumptions, our algorithm can ensure that the calibration error only increases by O(logB/T ). Notably, the calibration error scales logarithmically with the assumed upper bound logB, so choosing a loose bound B does not significantly degrade the guarantee.
The basic idea has two shortcomings illustrated in Figure 1 that we will address in the next sections. First, in a feasible quantile prediction, the 75% quantile should never be smaller than the 25% quantile, but our adjustments might violate this requirement because we adjust each quantile separately with no constraint. Second, the basic idea might lead to oscillation and instability. Intuitively, when the empirical frequency is incorrect, we make adjustments to correct it, but once the empirical frequency improves, we reduce the adjustments, which will cause the empirical frequency to become incorrect again, leading to oscillation. Empirically instability can hurt performance metrics such as sharpness.

heading: 3.2 IMPROVING THE BASIC ALGORITHM TO ENSURE FEASIBILITY
text: In this section, we address the feasibility problem. We first define some convenient notation shorthand. For our prediction Zt, we denote the distance between (neighboring) quantiles as ∆Ztk = Z t k+1 − Ztk for k = 0, · · · ,K. We can make a similar definition for the initial prediction ∆Z̃tk = Z̃ t k+1− Z̃tk. Under these new definitions, requiring feasibility is equivalent to requiring positive distance between neighboring quantiles, i.e. ∆Ztk > 0,∀k. The main challenge is that the basic adjustment algorithm Eq.(3) might not satisfy this new requirement, so we need to find a compromise. We draw inspiration from a mechanical system illustrated in Figure 2 where we imagine that our prediction quantiles Zt1, · · · , ZtK are a sequence of points in space. We imagine that each requirement (satisfying Eq.(3) or feasibility constraint) can be thought of as a “force” that “push” the quantiles to satisfy the requirement. These forces will contend with each other (because the requirements are conflicting), so the equilibrium of these forces is the natural compromise between the conflicting requirements.
Specifically, to enforce the requirement that the interval size ∆Ztk cannot be close to 0 or negative, we design a force that pushes the interval to be larger whenever it is small. Intuitively, we can think of this as inserting a “spring” (illustrated in Figure 2) into the interval. The spring will resist being compressed, and exert an opposing force whenever the interval becomes small. The smaller the interval, the more compressed the spring will be, so the larger the opposing force; in the limit where the interval size tends to 0, the opposing force will tend to∞. This prevents the interval size from going to zero (or negative) and ensures feasibility.
We add another desiderata: no force is generated only when the interval size ∆Ztk equals the initial prediction’s size ∆Z̃tk. Intuitively, we want to use the initial conformal prediction as much as possible (which we argued in Section 3.1 is desirable). This can be achieved by imaging that the spring has initial length ∆Z̃tk and resists both compressing and stretching.
To ensure that no prediction exceeds the bounds [−B,B], we also create two additional quantiles α0 = 0 and αK+1 = 1, and define Z0 = −B and ZK+1 = +B, so that assuming feasibility we will always have Z1 > Z0 = −B and ZK < ZK+1 = B. We formalize the above intuition as the following list of forces on each quantile:
1. Forces from miscalibration Etk defined by Eq.(3) and deviation from initial prediction Dtk := Z t k − Z̃tk. If we only have these two forces, we recover the basic algorithm. This is
because if we solve the equilibrium state of the systemDtk+E k t = 0, we get Z t k = Z̃ t k+E t k which is exactly Eq.(3).
2. Force from the spring in the interval above Atk := η(∆Z t k/∆Z̃ t k − ∆Z̃tk/∆Ztk) where η
is a hyper-parameter (we will show in Section 3.4 that setting these parameters is easy). Intuitively, if the interval above is larger than its initial size (i.e. ∆Ztk > ∆Z̃ t k) then the
spring “stretches“ which generates an upward force Atk > 0; conversely if the interval above is smaller than its initial size (i.e. ∆Ztk < ∆Z̃ t k) then the spring above “compresses“ which generates a downward force Ak < 0. If the interval has the same size, then no force is generated because Atk = 0.
3. Force from the spring in the interval below Btk := η(∆Z̃ t k−1/∆Z t k−1 −∆Ztk−1/∆Z̃tk−1).
We define the final prediction Ztk as the solution to the following equilibrium equations (i.e. all forces sum to 0 for all the quantiles)
Etk +D t k +A t k +B t k = 0,∀k = 1, · · · ,K. (4)
The resulting forecasting algorithm is provided as Algorithm 1.
Algorithm 1 Basic + Feasibility Forecasting Algorithm 1: For each k = 1, · · · ,K, initialize F 0k = 0 2: for t = 1, 2, · · · do 3: Observe Xt and use initial prediction algorithm (e.g. conformal calibration) to predict Z̃t. 4: For each k = 1, · · · ,K, compute Etk according to Eq.(3). 5: Use a physics simulation to find Ztk that solves Eq.(4) given Z̃ t and Etk, output Z t k.
6: Observe the true label Y t. For each k = 1, · · · ,K, set F tk = F t−1 k + I(Y t ≤ Ztk). 7: end for
Calibration Properties We will show that Algorithm 1 achieves calibration except when the “likelihood” is exponentially high (which is also very desirable). Specifically, we define a notion similar to likelihood
L(Zt, Y t) = (αk+1 − αk)/∆Ztk, for the k that satisfies Y t ∈ [Ztk, Ztk+1]. (5)
Intuitively, L(Zt, Y t) is large if the true label Y t belongs to a small interval (i.e. ∆Ztk is small). In the limit of infinitely many quantiles K → ∞, the quantile prediction becomes equivalent to a probability distribution, and Eq.(5) becomes the likelihood of the probability distribution. Theorem 2 (Simplified). Under the same condition as Theorem 1, for any c > 2/β, Algorithm 1 is b-calibrated where b = O(K/ √ T + c/T ) on any transcript such that ∀t, L(Zt, Y t) ≤ eΘ(c).
Theorem 2 is a simplified version of the full theorem (we replace all the constants with asymptotics). For a non-asymptotic version and proof see Appendix A. Intuitively, for any transcript generated by Algorithm 1, if the “likelihood” is never exponentially high, then the transcript is O(K/ √ T )- calibrated (c will be a bounded constant). When T is large, the calibration error will also go to 0 at a rate of O(1/ √ T ). On the other hand, if the likelihood is exponentially large, we are “almost” perfectly predicting the label, so there is little uncertainty. Uncertainty quantification metrics such as calibration become less relevant.

heading: 3.3 IMPROVING STABILITY WITH PID CONTROL
text: Finally we address the stability problem. As we argued, the key problem is the constant need to increase or decrease the adjustment Etk+1 based on whether the empirical frequency F t k is different from αkt. Our key insight is to think of this problem as a control problem: we would like F tk to match αkt and the error is measured by Etk (which we want to bring close to 0), so we can use control algorithms to stabilize the system. We show (heuristically) that a very popular control algorithm called proportional-integral-derivative (PID) is well suited for the task.
Our new algorithm replaces Etk in line 4 of Algorithm 1 with the weighted sum of three terms: proportional Etpk = Etk, integral Etik = ∑t τ=0E τ k and derivative Etdk = Etk − Etk−1. The proportional term Etpk is the same as Eq.(3), so if we only have the proportional term we would reduce to Algorithm 1. The integral term Etik is key to addressing the stability problem because it does not decrease even when theEtk goes to 0, so it applies a permanent adjustment when distribution drifts. The intuition of the integral term is demonstrated in Figure 2. Using the integral term is prone
to a well known problem called overshoot. Intuitively, the integral term will increase until Etk goes to 0, at which point the (permanent) adjustment might already be too large. The standard method to reduce overshoot is to use a derivative term Etdk (Minorsky, 1922), and we illustrate its effect in Figure 2. The overall adjustment is equal to
EtPID−k = kpE t pk + CLIP[−B,B] ( kiE t ik + kdE t dk ) (6)
where kp, ki, kd ≥ 0 are the hyper-parameters (again, we will show in Section 3.4 that choosing these hyper-parameters is not difficult). We also clip the integral and derivative terms to be within [−B,B] to ensure Corollary 1 even in the worst case scenario. For any real data in our experiments, the clipping never took effect.
We show that applying PID control does not break the calibration guarantee in Theorem 2. Corollary 1 (Simplfied). Under the same condition as Theorem 2, for any kp, ki, kd ≥ 0, Algorithm 1 with Etk replaced by E t PID−k in Eq.(6) is b-calibrated where b = O(K/ √ T + c/T ) on any transcript such that ∀t, L(Zt, Y t) ≤ eΘ(c).
A non-asymptotic version is also available in the appendix. With PID control Corollary 1 is essentially identical to Theorem 2 (in fact exactly the same for the simplified asymptotic version). Therefore, PID improves stability without degrading the calibration guarantees.

heading: 3.4 HYPER-PARAMETER SELECTION
text: Our algorithm consists of hyper-parameters which on first sight, seem difficult to choose for a practitioner. Specifically, the hyper-parameters are β, δ in Eq.(3), which specify how large is the adjustment Etk; the “elastic coefficient” η in the definition of A t k and B t k, and the PID coefficients kp, ki, kd.
Fortunately, we show that our algorithm is very robust to hyper-parameter choices. In fact, we tune these parameters on a few simple synthetic examples (similar to Figure 1), and fix them throughout all the experiments on real data. We show (in Figure 7 in the appendix) that even if we were to finetune these hyper-parameters specifically for each dataset, we would only see marginal improvements compared to using fixed hyper-parameters. Therefore, any practitioner can simply use the fixed hyper-parameters that we recommend, and generally do not have to worry about tuning them in practice.

heading: 4 RELATED WORK
text: In the regression setup, there are several notions of calibration. The most common one (which is studied in our paper) is probabilistic calibration (Gneiting et al., 2007; Kuleshov et al., 2018; Vovk et al., 2020). Other notions include marginal calibration (Gneiting et al., 2007), individual calibration (Zhao et al., 2020), distribution calibration (Song et al., 2019), threshold calibration, etc. Though not considered in our paper, calibration has also been extensively studied in the classification setup, calibration has been studied under (Brier, 1950; Murphy, 1973; Dawid, 1984; Cesa-Bianchi & Lugosi, 2006; Platt et al., 1999; Zadrozny & Elkan, 2001; 2002; Niculescu-Mizil & Caruana, 2005; Guo et al., 2017; Lakshminarayanan et al., 2017). There are also other notions of calibration such as multi-calibration (Hébert-Johnson et al., 2018), decision calibration (Zhao et al., 2021), classwise calibration (Kull et al., 2019), etc.
Many applications only require a confidence interval. For confidence intervals a notion very similar to calibration is coverage, i.e. the label should belong to the interval with high probability. Conformal prediction methods are the gold standard for predicting intervals with high coverage. (Vovk et al., 2005; Shafer & Vovk, 2008; Papadopoulos, 2008; Romano et al., 2019).
The above research primarily focus on the i.i.d. setup (or exchangeable setup), while many real world data are not i.i.d. In the online setup with distribution shift, regret minimization algorithms can guarantee asymptotic calibration (Cesa-Bianchi & Lugosi, 2006; Kuleshov & Ermon, 2017) for classification. For conformal prediction, when the likelihood ratio of the distribution shift is known, (Tibshirani et al., 2019) can still guarantee coverage. A concurrent work (Gibbs & Candès, 2021) can guarantee asymptotic coverage for conformal confidence intervals. Our work differs in that we focus on calibration, and solve additional feasibility and stability challenges.

heading: 5 EXPERIMENTS
text: This experiment aims to show on real time series data and standard regression benchmarks that our algorithm achieves strong calibration, sharpness and proper scores.
Datasets We use two time series prediction tasks (stock and COVID) and a regression benchmark with 17 datasets. For times series we start from expert predictions (such as from stock analysts or COVID forecasting teams), use conformal calibration to transform the expert predictions, and finally apply our algorithm to adjust the output of conformal calibration. For the regression benchmark we split each dataset into train and test. We use the training set to learn a prediction function offline, and assume that the test data are available sequentially online. For each dataset we add the following types of drifts: linear shift: add a constant bias to the label that increases over time; cycle: add a bias to the labels that moves up or down repeatedly; scale shift: multiple the label by a constant that increases with time; jump: add a constant bias at some ∀t > t0 for some time step t0. More details about the datasets are provided in appendix.
Performance Metrics : Calibration loss We use a notion similar to the expected calibration error (Guo et al., 2017)Lece = 1K ∑ k |FTk −αkT |. Pinball loss We use the average pinball loss (aka.
hinge loss) Lpinball(Y,Z) = 1K ∑ k αkReLU(Z −Yk) + (1−αk)ReLU(Yk −Z). The pinball loss is a proper score, i.e. the oracle quantile prediction minimizes the expected pinball loss. Sharpness We use the average interval size to measure the sharpness Lsharpness(Z) = 1K ∑ k |Zk − ZK−k|.
Baselines: We compare three baselines. original is the original expert predictions. We use kernel density estimation (Van Kerm, 2003) to transform the expert predictions into a probability; conformal is conformal calibration only (without Algorithm 1); regret is an adaptation of the regret minimization algorithm in (Kuleshov & Ermon, 2017) (details in the appendix).
Qualitative Results In Figure 3 and Figure 8,9 (in the Appendix) we plot some example prediction sequences and reliability diagrams. Our algorithm is able adjusts the predictions up or down to ensure calibration, and does not hurt sharpness (compared to conformal calibration only) when there is no distribution drift. In Figure 4 we plot the reliability diagram for the stock dataset. Our method is always calibrated (i.e. within the tolerance shaded in yellow), conformal calibration is uncalibrated (without the tolerance) for about 10% of the samples, while the regret minimization baseline is generally uncalibrated (because all time series are very short with around 50 samples).
Quantitative Results In Table 1 we show the performance of our algorithm compared to several baselines. Our method achieves the best calibration error (ece) by far, and comparable or better pinball loss compared to baselines. Even though the original prediction achieves the best sharpness overall, the pinball loss and calibration error of the original prediction are so large, that sharpness is meaningless. If we exclude the original prediction, then our method also achieves comparable or better sharpness than other methods that have reasonable calibration.
Downstream Task For COVID we simulate a response strategy (i.e. whether to close businesses) based on forecasted cases. We use the median of the prediction to select the best action. For details see Appendix C.1. The results are shown in Figure 5. Both conformal calibration and our method improves decision loss compared to using the original prediction; our method has an edge when the distribution changes rapidly, e.g. when cases surge.

heading: 6 CONCLUSION
text: This paper introduces a novel type of online algorithms to ensure calibration, where we start with an algorithm that guarantees calibration in the i.i.d. setup, and make modifications whenever it fails to achieve calibration. Future work can explore its application to other calibration definitions or more generally, other online learning problems.

heading: A PROOFS
text: Theorem 1. Let Z̃tk be generated by any forecasting algorithm and bounded in [−B,B], for any δ > 0, the prediction algorithm defined by Eq.(3) is b-calibrated where b(T ) = b∗δ(T )+ log(2B+1)+β βT .
Proof of Theorem 1. We prove that at any time t we cannot have |F tk − αkt| > b∗δ(t) + log(1 + 2B)/β + 1. This is equivalent to two statements that we prove separately by induction
F tk ≥ αkt− b∗δ(t)− log(1 + 2B)/β − 1 (7) F tk ≤ αkt+ b∗δ(t) + log(1 + 2B)/β + 1 (8)
To prove Eq.(7) suppose at time t we have F tk ≥ αkt− b∗δ(t)− log(1 + 2B)/β − 1, then there are two situations:
Case 1. F tk ≥ αkt− b∗δ(t)− log(1 + 2B)/β, then we naturally have
F t+1k ≥ F t k ≥ αkt− b∗δ(t)− log(1 + 2B)/β > αk(t+ 1)− b∗δ(t+ 1)− log(1 + 2B)/β − 1
(9)
Case 2. F tk < αkt− b∗δ(t)− log(1 + 2B)/β, then by Eq.(3) we have
Et+1k > 1 + 2B − 1 = 2B By the boundedness assumption,
Zt+1 ≤ B < −B + Et+1k < Z̃ t+1 k + E t+1 k = Y t+1 k
so we must have F t+1k = F t k + 1, so
F t+1k = F t k + 1 ≥ αkt− b∗δ(t)− log(1 + 2B)/β > αk(t+ 1)− b∗δ(t)− log(1 + 2B)/β − 1
(10)
Similarly to prove Eq.(8) suppose at time t we have F tk ≤ αkt + b∗δ(t) + log(1 + 2B)/β + 1 then at time t+ 1 we must have
F t+1k ≤ αk(t+ 1) + b ∗ δ(t+ 1) + log(1 + 2B)/β + 1 (11)
Theorem 2. Suppose the initial predictions are feasible (i.e. Z̃tk < Z̃tk+1,∀t, k). For some c > 2/β, Algorithm 1 is b-calibrated where b(T ) = (2Kb∗δ(T ) + Kc/T, 0) on any transcript that satisfies ∀t ≤ T, L(Y t, Zt) ≤ 12ηB(K+1) (e β(c−1)/2 − 4B).
Proof of Theorem 2. We first sketch our proof strategy. The proof will follow 3 steps:
1. We prove that if an interval contained the label too often, there will be a large force difference Etk+1 − Etk. 2. We prove that if there is a large force difference, then the interval will be small.
3. If at time t an interval contains the label N times too often, at some previous time step t′ ≤ t, it must have contained the label N − 1 times too often and contained the label once again. However, at t′ the interval must have been small, hence the likelihood of the label must be large. Therefore, unless the likelihood is large, no interval can contain the label too often.
First we show a simple Lemma based on boundedness assumptions. This Lemma is a simple consequence of the boundedness Y0 = −B, YK+1 = B. We do not use this Lemma directly but it will be used in the proof of the other Lemmas.
Lemma 1. ∀t, ∀k = 0, · · · ,K + 1, Atk < B,Btk > −B, |Dtk| < B
The following Lemma shows that intuitively, if the ground truth label belongs to an interval too often, then there Etk+1 must be significantly less than E t k
Lemma 2. For any c > 2/β, if F tk+1 − F tk = ∆αkt+ 2b∗δ(t) + c then E t+1 k+1 − E t+1 k < −eβc/2.
Next we show that if the forcesEtk+1 is significantly smaller thanE t k, then intuitively, this difference in force will “compress“ the interval, making it small.
Lemma 3. For any L > 0, if Etk+1 − Etk < −L we have ∆Z̃tk/∆Y tk > L−4B2η .
Based on these lemmas we can complete the proof. Suppose for some k ∈ 1, · · · ,K and c > 2 we have F tk+1 − F tk = ∆αkt+ 2b∗δ(t) + c, then consider last time that the label belongs to the interval, i.e. we want to find the largest t′ ≤ t such that Zt′ ∈ [Y t′k , Y t ′ k+1). For t ′ we have
F t ′−1 k+1 − F t′−1 k = F t k+1 − F tk − 1 (12)
> ∆αkt+ 2b ∗ δ(t) + c− 1 (13) = ∆αkt ′ + 2b∗δ(t
′) + (c+ ∆αK(t− t′) + 2b∗δ(t)− 2b∗δ(t′)− 1) (14) ≥ ∆αkt′ + 2b∗δ(t′) + c− 1 (15)
Where the last inequality is because b∗δ is monotonically non-decreasing in t. Therefore by Lemma 2 we have
Et ′ k+1 − Et ′ k < −eβ(c−1)/2 (16)
So by Lemma 3 we can conclude
∆Z̃t ′ k /∆Y t′ k > eβ(c−1)/2 − 4B
2η (17)
which implies whenever eβ(c−1)/2 > 4B, or c > 2 log(4B)/β + 1
∆Y t ′ k < 2η∆Z̃ t′ k (e −β(c−1)/2 − 4B) < 2ηB(e−β(c−1)/2 − 4B) = 2ηB(e−β(c−1)/2 − 4B) (18)
Hence. by the definition of the “likelihood“
L(Y t ′ , Zt ′ ) > αk+1 − αk 2ηB (eβ(c−1)/2 − 4B) = 1 2ηB(K + 1) (eβ(c−1)/2 − 4B) (19)
In other words, unless ∃t′ such that L(Y t′ , Zt′) > 12ηB(K+1) (e β(c−1)/2 − 4B), we cannot have F tk+1 − F tk = ∆αkt+ 2b∗δ(t) + c.
Finally we connect the bound on F tk+1 − F tk with the quantile calibration error.
Lemma 4. For any b > 0, if ∀k ∈ K, F tk+1 − F tk ≤ ∆αkt+ b, then
max k |F tk − αkt| ≤ Kb (20)
So combined we have if ∀t, L(Y t, Zt) ≤ 12B(K+1) (e (c−1)/2 − 4B), then
max k |F tk − αkt| ≤ 2Kb∗δ(t) +Kc (21)
Finally we prove the lemmas used in this theorem.
Proof of Lemma 2. Suppose F tk+1 − F tk = ∆αkt + 2b∗δ(t) + c We first use simple algebraic transformation to get
(F tk+1 − αk+1t)− (F tk − αkt) = c+ 2b∗δ(t) (22)
For notation simplicity denote the first term LHS as u = F tk+1 − αk+1t and second term as v = F tk − αkt, we can simplify the equation as so u− v = 2b∗δ(t) + c. First we consider the scenario where the interval is not the first interval α0, α1 or the last interval αK , αK+1, in other words, we consider the situation where k ∈ [1,K−1]. There are three different
possible situations: u ≥ 0, v ≥ 0 and u ≥ 0, v < 0 and u < 0, v < 0. Note that we can never have u < 0, v ≥ 0. We separately consider the three situations. 1. If u ≥ 0, v ≥ 0 then
Et+1k+1 − E t+1 k = −e
β(u−b∗δ(t)) + 1− (−eβ(v−b ∗ δ(t)) + 1) (23)
= −eβ(u−b ∗ δ(t)) + eβ(v−b ∗ δ(t)) ≤ 1− eβ(c+b ∗ δ(t)) (24)
where the final inequality achieves equality when u = c+ 2b∗δ(t), v = 0
2. If u ≥ 0, v < 0 then
Et+1k+1 − E t+1 k = −e
β(u−b∗δ(t)) + 1− (e−β(v+b ∗ δ(t)) − 1) (25)
= −eβ(u−b ∗ δ(t)) − e−β(v+b ∗ δ(t)) + 2 < −2eβc/2 + 2 (26)
where the final inequality achieves equality when u = (c+ 2b∗δ(t))/2, v = −(c+ 2b∗δ(t))/2. 3. If u < 0, v < 0 then
Et+1k+1 − E t+1 k = (e
−β(u+b∗δ(t)) − 1)− (e−β(v+b ∗ δ(t)) − 1) (27)
= e−β(u+b ∗ δ(t)) − e−(v+b ∗ δ(t)) ≤ 1− eβ(c+b ∗ δ(t)) (28)
where the final inequality achieves equality when u = 0, v = −c− 2b∗δ(t). Next we consider the situation where k = 0. In this case v = F0 − α0t = 0 by assumption. So
Et+11 − E t+1 0 = E t+1 1 = −eβ(c+2b
∗ δ(t)−b ∗ δ(t)) + 1 = −eβ(c+b ∗ δ(t)) + 1 (29)
Next we consider the situation where k = K. Similar to the previous argument, u = 0 by assumption, so
Et+1K+1 − E t+1 K = −E t+1 K = −(e
β(c+2b∗δ(t)−b ∗ δ(t)) − 1) = 1− eβ(c+b ∗ δ(t)) (30)
Combined we have ∀k, we take the maximum of the lower bounds and have
Et+1k+1 − E t+1 k < max(−2e
βc/2 + 2, 1− eβ(c+b ∗ δ(t))) (31)
Since c > 2/β we can further simplify this expression by observing −2eβc/2 + 2 < −eβc/2 and 1− eβ(c+b∗δ(t)) < 1− eβc < eβc/2, so we have
Et+1k+1 − E t+1 k < max(−2e
βc/2 + 2, 1− eβ(c+b ∗ δ(t))) < −eβc/2 (32)
Proof of Lemma 3. Observe the equilibrium equationDk+Bk+Ak+Ek = 0 andDk+1 +Bk+1 + Ak+1 + Ek+1 = 0; by Lemma 1
Ak < B,Bk > −B, |Dk| < B
so consequently we have
Ak −Bk+1 = Dk+1 + Ek+1 +Bk+1 −Dk − Ek −Ak < −L+ 4B
By the definition of Ak and Bk+1 we have
2η
( ∆Yk
∆Z̃k − ∆Z̃k ∆Yk
) < −L+ 4B (33)
which implies that
∆Z̃k ∆Yk > L− 4B 2η + ∆Yk ∆Z̃k > L− 4B 2η (34)
Proof of Lemma 4. We prove the contra-positive statement of the Lemma. If maxk |F tk−αkt| ≤ Kb is violated, there must be one of two situations:
1. there exists some k such that F tk − αkt ≥ Kb then
Kb ≤ F tk − αkt− F t0 + α0t F t0 ≡ 0, α0t ≡ 0 (35) = ∑
0≤j≤k
(F tj − F tj−1)− (αj − αj−1)t Telescope (36)
≤ (k + 1) max j≤k (F tj − F tj−1)− (αj − αj−1)t (37)
2. there exists some k such that F tk − αkt ≤ −b then
Kb ≤ −F tk + αkt+ F tK+1 − αK+1t F tK+1 ≡ t, αK+1t ≡ t (38) = ∑
k≤j≤K+1
(F tj − F tj−1)− (αj − αj−1)t Telescope (39)
≤ (K + 1− k) max j≥k (F tj − F tj−1)− (αj − αj−1)t (40)
In both cases, there must be some j such that
(F tj − F tj−1)−∆αj−1t ≥ b (41)
Corollary 1. Suppose the initial predictions are feasible (i.e. Z̃tk < Z̃tk+1,∀t, k). For some c > 2/β, the prediction algorithm in Section 3.3 is b-calibrated where b(T ) = 2Kb∗δ(T ) +Kc/T on any transcript that ∀t, L(Y t, Zt) ≤ 12ηB(K+1) (e β(c−1)/2 − 6B).
Proof of Corollary 1. To prove this corollary we only have to modify Lemma 2, and the rest of the proof follows as Theorem 2.
Lemma 5. For any c > 2β, if F tk+1 − F tk = ∆αkt + 2b∗δ(t) + c then E t+1 PID−k+1 − E t+1 PID−k < −eβc/2 + 2B.
Lemma 5 holds because Et+1PID−k differs from E t+1 k by at most B, so
Et+1PID−k+1 − E t+1 PID−k < E t+1 k+1 − E t+1 k + 2B < −e βc/2 + 2B (42)

heading: B ADDITIONAL EXPERIMENT RESULTS
text: B.1 HYPER-PARAMETER TUNING
We tune hyper-parameters that minimize a weighted combination of the expected calibration error (ece) and the pinball loss.
Lmix(wmix) = wmixLece + (1− wmix)Lpinball (43)
where wmix is a weight coefficient we choose to be 0.4. We use Optuna (Akiba et al., 2019) to seek promising hyper-parameter candidates. We try 100 hyper-parameter candidates and pick the best hyper-parameter that minimizes the loss function Lmix(wmix). As we mentioned in Section 3.4, we use simple synthetic tasks for tuning. We first define the label of a task as follows:
Y t = a+ b t, t ∼ N (0, 1) (t = 1, . . . , T ) (44)
where coefficients a and b are sampled by a ∼ N (0, 1) and b ∼ U(0, 1). We also set predictions Z̃t = F−1(αk) where F is the cumulative distribution function of standard Normal distribution.
After tuning, we evaluate the best hyper-parameters for each wmix on other synthetic tasks. Table 2 compare performance metrics between different wmix. Large wmix leads to small calibration loss values and large pinball loss/sharpness values. We adopt wmix = 0.4 since the hyperparameter provides good calibration and pinball values. The obtained hyper-parameter values are δ = 0.47, β = 0.16, kd = 0.08, kp = 1.00, η = 0.96. As for ki, we assume ki depends on quantile level k as follows:
ki,k = ki,max − (ki,max − ki,min)|1− 2(k − 1)/(K − 1)|, ∀1 ≤ k ≤ K. (45)
The obtained hyper-parameters are ki,max = 0.09, ki,min = 0.04.
Next we verify that our algorithm is robust to hyper-parameter choice. In Figure 6 we compare the performance of a hyper-parameter on synthetic sequences and real data. They are highly correlated, which indicates that hyper-parameters that do well on synthetic sequences also perform well on real data.
Finally, we demonstrate that the adopted hyper-parameter in the paper works well compared to “optimal” hyper-parameters for each dataset. Concretely, we tune hyper-parameter on regression benchmark datasets with distribution drifts, and compare the tuned hyper-parameter and the recommended ones (which are tuned on synthetic sequences). Figure 7 shows that tuning on each dataset separately (which is only possible for off-line data) only improves performance marginally compared to using fixed hyper-parameters. This suggests that the recommended hyper-parameter is a reasonable choice in a variety of cases.
B.2 VISUALIZATIONS OF PREDICTIONS ON VARIOUS DATASETS
In addition to Figure 3 in Section 5, we visualize additional example predictions on COVID and regression benchmark datasets. Figure 8 shows example predictions on the COVID prediction dataset, and Figure 9 shows example predictions on regression benchmark datasets with distribution drifts. Our method adapts to distribution drifts.

heading: C EXPERIMENT DETAILS
text: Stock Dataset For stock we predict quarterly earnings per share (EPS) for S&P 500 stocks from 2010/01 to 2021/03. We use analyst earnings forecasts provided by I/B/E/S as the original prediction. To remove stocks with low quality data, we apply the following filtering rule for each stock.
1. exclude analyst forecasts whose missing rates are over 20% 2. if the number of analysts is less than five, remove the stock 3. if all remaining analysts do not predict earnings at least a quarter, remove the stock
By the above procedure, 279 stocks are selected from 500 stocks in S&P 500. The average number of analysts per stocks is 9.3 and each time series length is 45 data points.
COVID dataset For COVID we predict COVID cases for the next week. We use expert forecasts from the COVID forecast hub (Ray et al., 2020). Because the data is heavy tailed, we transform all forecasts and labels by Z 7→ log(1 + Z) and predict the transformed label. We select the teams and locations that do not have missing predictions during the period between Aug 2020 and Aug 2021 (6 teams) as the original prediction, and randomly sample 100 locations. The time series length is 56 data points (so approximately one year of data).
Regression Benchmarks We use the following UCI datasets: blog, boston, concrete, crime, energy-efficiency, fbcomment-1, fbcomment-2, forest-fires, mpg, naval, power-plant, protein, superconductivity, wine, and yacht. In addition to UCI datasets we also use two common regression datasets in the uncertainty quantification literature, which are kin8nm and medical-expenditure. We use 100 samples for each test dataset, i.e., the length is 100.
Distribution drifts For the regression benchmarks we consider four types of artificial distribution drifts. The definition of the distribution drift is given by
• linear shift: add linear drifts to the labels, i.e., Y t = Y t + βt • scale shift: change the scales of the labels, i.e., Y t = Y t(1 + β √ t)
• jump: jump the mean at the halfway point, i.e., Y t = Y t (t < T/2), Y t = Y t + β (t ≥ T/2)
• cycle: add fluctuations to the labels, i.e., Y t = Y t + β sin(2πt/T )
In the experiments, we set coefficient β = 0.1, 1.0, 3.0, 3.0 for linear shift, scale shift, jump, and cycle, respectively.
C.1 DOWNSTREAM TASK DETAILS
COVID Decision Making We consider the following simple decision loss: for each county, there are four tiers of response 0,1,2,3. A higher tier response is needed whenever there are more forecasted cases. Let a denote the response, and z denote the (log) cases, we choose the loss function as the sum of COVID losses and shut-down losses, respectively defined as
a = 0 a = 1 a = 2 a = 3 COVID loss z ≤ 4 0 0 0 0
4 < z ≤ 6 1 0 0 0 6 < z ≤ 8 3 2 0 0 z > 8 6 5 3 0
Shutdown loss 0 0.3 0.7 1.0
C.2 ONLINE LEARNING BASELINE
We adopt Algorithm 1 in (Kuleshov & Ermon, 2017) based on internal regret minimization, we use exactly the same algorithm except that (Kuleshov & Ermon, 2017) minimizes the L2 loss or Brier loss (which are proper scoring rules for the conditional expectation), while we minimize the
pinball loss (which is the proper scoring rule for quantiles). In addition, we separately calibrate the different quantiles Y1, · · · , YK (i.e. for each k we run a separate instance of modified Algorithm 1 in (Kuleshov & Ermon, 2017)).
The original procedure in (Kuleshov & Ermon, 2017) can guarantee calibration in the binary classification setting; we give intuition on why our modification can guarantee calibration in the regression setting. we first bin the possible regression outputs into M equal sized bins −B = y0 < · · · < yM = B. As a notation convenience, denote the center of each bin as ycm = ym+ym−1 2 .
Consider the αk-th quantile, our goal is to minimize the interval regret under the pinball loss, define by ∀c ∈ [−B,B]
Rcinterval = ∑ t I(Y tk = c)Lpinball(c, Zt)− inf y ∑ t I(Y tk = c)Lpinball(y, Zt)
Intuitively, among the time that the true label equals c there is not a better alternative prediction y that achieves lower pinball loss. Because pinball loss is a proper scoring rule, Rcinterval = 0 implies that c is the αk-quantile of {Zt | Y tk = c}, i.e. for this subset of samples, the true label Zt is below the prediction Y tk = c for αk proportion of times. This is exactly the definition of calibration. To implement the algorithm realistically, we have to choose a suitable binning strategy, as well as a internal regret minimization algorithm. We use the same choice as the experiments in (Kuleshov & Ermon, 2017). Note that to choose the binning strategy we have to know the lower bound and upper bound of the data. We give the baseline an unfair advantage by choosing the lower bound as the lowest value in the data, and the upper bound as the highest value in the data (up to a 10% margin).

Review summaries:
This paper considers the problem of calibrating quantiles outputted by a regression model. The data is assumed to draw from a non-iid source (e.g. time series with distribution shift). The paper proposes a method that adjusts the output quantiles to be better calibrated for a less-restricted definition of b-calibrated. Some further methods are also proposed to ensure the quantiles are monotonic and stable. Experiments are conducted to compare the proposed method with uncalibrated models and conformal predictions.

For real-world applications where ML models are used, quantifying predictive uncertainty is an essential task as, for example, in safety-critical applications misclassification might have disastrous consequences. On the theory side, the analysis is typically performed under the i.i.d. assumption which could be unsatisfactory as deployed models inevitably encounter changes in the data generating distribution and/or certain dependencies that invalidate the results established for the i.i.d. setting. Focusing on regression, the authors design an adaptive (to distribution shifts) procedure that satisfies a certain calibration guarantee. The authors study the empirical performance of the proposed procedure on a collection of time-series and regression datasets.