# Simple constant-width prediction interval for a regression model

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:

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.

• I am inclined to say no it will not because $f$ is just an estimate of the true data generating process, so come 'production' time there will be an error between your model and the true data generating process which will throw off the coverage probabilities. May 21, 2021 at 19:31
• @David: couldn’t this be easily „fixed“ by reformulating the question as “does $|P[y \in […]] - p|$ converge to zero as the size of the training set increases, right? Maybe the formula is not the right one (I am far from being an expert in these convergence in probability stuff) but I think the direction is the right one… May 21, 2021 at 20:07
• @Fabian, so if we knew the true $f$ then the empirical distribution function of $e_i$ would converge almost surely to the true distribution function of the error by Glivenko-Cantelli. I have to think about it more, but I think if you assume as $n$ increases the $f$ we fit gets close to the true $f$ then you would get the same result (although maybe not as strong convergence?). May 21, 2021 at 20:46
• @David I am interested in the case where $f$ is an estimate learned from a training set (step 1 in my question). We definitely do not know the true regression function $E[Y | X]$. However, we have access to a test set from which we can calculate the empirical distribution of our prediction errors. May 21, 2021 at 22:10
• @Adrian the more I have thought about this, the more I would agree that what you are suggesting will work. Since what will happen is for large enough $n$ at each step of train/test/production you are getting iid samples of $e_i$ which are ultimately what you care about in your prediction intervals. May 25, 2021 at 2:02

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.$$

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.

• which is true as long as $f(x_i)$ does not depend on points in the test set other than $x_i$. I suppose there are no other points in the test set besides the ones indexed by $i$. Did you perhaps mean something like does not depend on variables in the test set other than $x$? May 28, 2021 at 12:03
• I meant that for any particular $i$, $f(x_i)$ should depend only on $x_i$, and not any other $x_j$, $j \neq i$, in the test set. May 28, 2021 at 12:10