# How to generate 95% prediction interval around predictions from ML model?

I have predictions from an ML model and would like to generate 95% prediction intervals around each prediction generated from the model such that I can claim that these are the plausible range of values one might expect for each observation.

I have seen some resources indicate it should be:

Yhat +- 1.96 * std(residuals)

I have seen another resource indicate it should be yhat that is used instead of residuals for the range:

Yhat +- 1.96 * std(yhat)

And I have seen yet another resource indicate it should be standard error and not standard deviation

Yhat +- 1.96 * se(Yhat)

I have also seen other resources indicate that you need to fit the Sklearn model with a quantile parameter set to .05 to generate upper and lower bounds, suggesting that I need to train the model again and it’s something you do during training.

Overall, I am now a bit confused and uncertain about the best or most appropriate way to generate a prediction interval around predictions to reflect the plausible range of values around each prediction. Does anyone have any recommendations?

• You may want to look up "Conformal prediction" Commented Jul 13 at 17:59
• I don’t understand the downvote. Sure, the title sounds like a request for the function argument or lines of code that give the prediction intervals, but the question itself is about the statistics.
– Dave
Commented Jul 13 at 19:05
• Obviously, all these suggestions make the assumption that the prediction error is normally distributed. Is this a reasonable assumption for your model and are your (leave-one-out) residuals actually normally distributed? Commented Jul 14 at 11:42
• @qwr How do you propose to use bootstrap to solve this? Bootstrap is typically used for confidence intervals, which are not the same as the prediction intervals that seem to be the subject of this question.
– Dave
Commented Jul 15 at 17:37
• @Dave bootstrapping can be used to estimate almost anything about a model. One approach saattrupdan.com/posts/2020-03-01-bootstrap-prediction
– qwr
Commented Jul 15 at 17:47

THIS SEEMS TO BE AN OPEN PROBLEM. Let's look at some possible solutions and their drawbacks.

First, you propose this Yhat +- 1.96 * std(residuals). Let's put that in more standard statistics notation and then break down what the pieces mean.

$$\hat y_i \pm 1.96\hat\sigma_{\hat\varepsilon}$$

$$\hat y_i$$ is the point prediction.

$$1.96$$ comes from the $$95\%$$ bounds of a Gaussian distribution.

$$\pm$$ means that you go $$1.96\hat\sigma_{\hat\varepsilon}$$ above and $$1.96\hat\sigma_{\hat\varepsilon}$$ below the point prediction.

$$\hat\sigma_{\hat\varepsilon}$$ is the (constant) standard deviation of an implicit error term, estimated by the standard deviation of the residuals.

Thus, this formulation implies that the way to get the $$95\%$$ prediction interval is to use the point estimate $$\hat y_i$$ as the mean of a conditional Gaussian and the residual variance as the variance of that conditional Gaussian. We then take the $$0.025$$ and $$0.975$$ quantiles of this $$N\left(\hat y_i, \hat\sigma^2_{\hat\varepsilon}\right)$$ conditional distribution.

Under the assumption that there is a Gaussian error term with constant variance, this makes some sense. However, that is extremely restrictive. What about if the error term is not Gaussian? What about if the error term lacks a constant variance? What if the non-Gaussian error term lacks constant variance? Under this quite plausible deviations from the $$N\left(\hat y_i, \hat\sigma^2_{\hat\varepsilon}\right)$$ ideal, the prediction interval could wind up being terrible.

Second, you propose Yhat +- 1.96 * std(yhat). I interpret this to mean $$\hat y_i\pm 1.96\sigma_{\hat y_i}$$.

$$\hat y_i$$ is the point prediction (again).

$$1.96$$ comes from the $$95\%$$ bounds of a Gaussian distribution (again).

$$\pm$$ means that you go $$1.96\hat\sigma_{\hat y_i}$$ above and $$1.96\hat\sigma_{\hat y_i}$$ below the point prediction.

This is somewhat less restrictive than the earlier formulation, as it allows for each $$i$$ to have a different variance. However, the $$1.96$$ still implies a Gaussian conditional distribution that might not be reasonable. Additionally, estimating the variance at each $$i$$ is not entirely straightforward.

Third, you propose Yhat +- 1.96 * se(Yhat). I think this is the same as the previous formulation but with some confusion about the difference between standard deviation and standard error. The standard error should shrink toward zero as the sample size increases. However, we do not want such behavior. That's why prediction intervals, such as the first proposed formulation, do not shrink to have zero width.

Fourth, you propose to fit the Sklearn model with a quantile parameter set to .05 to generate upper and lower bounds. My suspicion is that this implements the first of the proposed formulations, complete with the caveats about non-Gaussian conditional distributions and non-constant error-term variance.

Another idea is to use explicit quantile modeling, that is, quantile regression. This has appeal, since quantile modeling is not bound to a particular distribution of the error term the way the above techniques are bound to a Gaussian error term. Sure, you could replace that $$1.96$$ with, say, the quantiles of a t-distribution if you are comfortable assuming t-distributed errors with a particular number of degrees of freedom, but even that is quite restrictive unless you really believe that is conditional distribution.

By using explicit quantile modeling, you let the data tell you what the conditional quantiles are, rather than making assumptions about particular distribution families that might be wildly inaccurate.

However, quantile modeling is not a silver bullet! Flexibility in the model is needed to allow for prediction intervals from the quantile models to vary in width. You might want wide prediction intervals in certain regions and narrowed intervals in others. For instance, in the (granted, contrived) example below, there would have to be flexibility in the modeling to capture the oscillating error variance.

set.seed(2024)
N <- 1000
x <- runif(N, -pi, pi)
Ey <- 1 - x
e <- rnorm(N, 0, sin(x)^2)
y <- Ey + e
plot(x, y)


"Fine, I'll be throwing a bunch of flexibility at the problem by using deep learning, anyway," you say? That flexibility can result in what I remember an arXiv paper referring to as the "embarrassing" problem of quantile crossover, where high predicted quantiles like $$0.975$$ are lower than low predicted quantiles like $$0.025$$. That is, you predict the lower limit of the prediction range to be $$7$$ while the upper limit of the prediction range to be $$1$$. Maybe you just flip them in that situation to say the prediction interval is $$(1, 7)$$, but then you are not using the explicit quantiles, and the prediction interval endpoints lack their usual interpretation and meaning. They might not wind up being very good prediction interval endpoints!

That arXiv paper made attempts to keep quantile crossover from happening in their quantile modeling. I recall these attempts not being entirely successful.

Thus, I call this an open problem that lacks a definite solution.

• I find this to be a fairly decent discussion of using quantile regression neural networks to produce prediction intervals, despite my reservations about most content on Toward Data Science.
– Dave
Commented Jul 15 at 17:22

Simple linear regression is parametric, since you model $$Y = \beta X + \epsilon$$ with the assumption of $$\epsilon \sim N(0, \sigma^2)$$. The prediction interval for a new $$\hat y_h$$ is bigger than the confidence interval for estimated mean $$\hat \mu$$, because there's two sources of variance: $$\sigma^2$$ from the error term and $$\sigma^2(Y_h)$$ from estimating $$\mu_Y$$. So the prediction interval will be $$\hat y_h \pm 1.96 \sqrt{\sigma^2 + \sigma^2(Y_h)}$$ See notes on prediction intervals for a new response. Even without modeling the distributions, the two sources of noise in prediction intervals carry through to nonparametric models.

What you are looking for is nonparametric prediction intervals, which has lots of existing literature already. Major techniques: