Simple constant-width prediction interval for a regression model

Question

Consider the following approach to generating prediction intervals for a regression problem:

1. Train a regression model on a training set. Let $f$ denote the fitted model, i.e. $f(x_i)$ is the model's prediction given inputs $x_i \in \mathbb{R}^k$.
2. Compute the model's predictions and prediction errors $e_i \equiv y_i - f(x_i)$ on a test set (i.e. not the training set used to fit $f$ in step 1).
3. Define a desired prediction interval coverage probability $p \in \left[0, 1\right]$. For example, we might set $p=0.8$ if we want 80% of future observations to fall inside our prediction intervals. Calculate the $\tau_{\text{lower}} \equiv (1-p)/2$ and $\tau_{\text{upper}} \equiv (1+p)/2$ quantiles of the test set errors. With $p=0.8$ we'd have $\tau_{\text{lower}} = 0.1$ and $\tau_{\text{upper}} = 0.9$, so we'd be calculating the 10th and 90th percentiles of our errors on the test set. Let $Q^\text{upper}$ and $Q^\text{lower}$ denote the error quantiles, and note that we typically expect $Q^\text{lower} < 0 < Q^\text{upper}$.
4. At prediction time (when making new predictions "in production"), construct prediction intervals $\left[f(x_i) + Q^\text{lower}, f(x_i) + Q^\text{upper}\right]$. The intent is for $y_i$ to be inside this interval with probability $p$.

Note that the intervals described above are "simple" in the sense that their width $Q^\text{upper} - Q^\text{lower}$ does not vary with $x_i$. We compute the width in step 3 using our test set errors, and it is fixed going forward to step 4.

Question: If we assume that the data in all steps (the training set, the test set, and the "production" data) are i.i.d. samples from the same data generating process, i.e. from the same joint distribution over $X$ and $Y$, will the prediction intervals described in step 4 have the desired coverage in expectation? Will we have $\Pr\left[y_i \in \left[f(x_i) + Q^\text{lower}, f(x_i) + Q^\text{upper}\right]\right] \approx p$?
(Note that this probability is not conditional on $x_i$. There could be values of $x_i$ for which $\Pr\left[y_i \in \left[f(x_i) + Q^\text{lower}, f(x_i) + Q^\text{upper}\right] \, | \, x_i\right] \neq p$, and that would be fine as long as we have the desired coverage on average, taking an expectation over all future observations).

If so, does this technique have a name? Would it work for any arbitrary regression model (linear regression, gradient boosting, neural nets, etc), and for any distribution of $Y \,|\, X$ (not necessarily Gaussian)?

---

Related questions:

• Difference between confidence intervals and prediction intervals
• Bootstrap prediction interval
• How do I calculate prediction intervals for random forest predictions?

---

Python 3.8 simulation (using CatBoostRegressor). I originally tried copying the code directly into a <pre> block, but I couldn't get it to render properly.

The technique I am describing appears to work in the simulation: 90% prediction intervals contain approximately 90% of "production" Y values. The simulation also shows that this is true only on average: if we condition on X, we can find subsets of the predictor space where the intervals have lower coverage.

Answer 1

As described in your question, you should expect that the prediction intervals in production to have the desired coverage.

During training you come up with a fitted model $f$, which need not be equal to, or even close to, the true data generating function.

Using the test set you get a sample $e_1,\dots ,e_{n_\text{test}}$. By Glivenko-Cantelli as $n_\text{test}\rightarrow \infty$ the empirical CDF of $e_1,\dots ,e_{n_\text{test}}$ should converge to the true CDF of $e_i$.

Because of the above result, using the empirical CDF you can come up with an interval such that for a new independent set of observations in production, for each $e_i^*=y_i^*-f(x_i^*)$, you will have $$P\left(e_i^*\in[e_{low},e_{high}]\right)\approx p \Rightarrow P\left(y_i\in[e_{low}+f(x_i),e_{high}+f(x_i)]\right)\approx p.$$

Answer 2

I think the fundamental thing to realize here is that if your test set is really a test set (used in no way for training) and your test data are really i.i.d., and test and production really have the same distribution, then $f$ can be considered as any other function determined independently of your data. The fact that $f$ was originally constructed using a training set is interesting biographical information about the function, but is irrelevant to the problem.

In particular, if $\{(x_i, y_i)\}$ are i.i.d. and drawn from the same distribution as production, then $\{e_i := f(x_i) - y_i\}$ are i.i.d. and from the same distribution as production. This is where we use the property that the test set was not used for training. If $f(x_i)$ were somehow dependent on $x_j$, $j \neq i$, as would be the case if the test set were used for training or model selection, then $\{e_i\}$ would not be mutually independent. And if $y_i$ were used to determine the value $f(x_i)$, again as would be the case if the test set were used for training, the distribution of $y - f(x)$ in the test set would be different than in production.

The empirical distribution function of $e$ converges to its true distribution function as the size of your test set grows, by the Glivenko-Cantelli theorem, and so your prediction intervals are correct.

Worth pointing out explicitly that the convergence has nothing to do with the any properties of the training set or the regression model, although maybe the speed of convergence does. The training set itself could be non-i.i.d., or even drawn from a different distribution than the test set and production. The only "condition" for the Glivenko-Cantelli theorem is that the sample used to construct the empirical distribution function is drawn independently from a single distribution.