# How to create a prediction interval with the fact that the residuals follow a specific distribution (in python)

I am looking at a software development pipeline where I am predicting the lead time of different products flowing through the pipeline.

After applying a boxcox transformation on the lead time (target variable) and creating a XGBoost regressor model I can see that the residuals follow a t-locationScale distribution.

So now I looked at this guide which describes a method to create a prediction interval for any regression model assuming that the residuals are normally distributed. https://qucit.com/a-simple-technique-to-estimate-prediction-intervals-for-any-regression-model_en/

But I tried to tweak it to my distribution.

So a t-locationScale distribution has a $$\sigma$$, $$\mu$$ and $$\nu$$ parameter. The variance is only defined for $$\nu>2$$. My specific distribution has $$\nu = 2.56$$ and $$\mu = 0.04$$, $$\sigma = 0.97$$ So I could take the 95% interval of this distribution and say that for any $$\hat{y}$$, the prediction interval is the 95% interval of the residual distribution.

But I want to take into consideration that the prediction interval should change with different inputs. I created a regressor model, which I trained and then made predictions using the validation set. I then took the square of the error and trained an additional error model on this data. Such that the error model could predict the variance of the residuals distribution.

  xgb = XGBoostRegressor()
xgb.fit(X_train,y_train)
y_hat = xgb.predict(X_val)
val_error = (y_hat-y_val)**2

xgb_error = XGBoostRegressor()
xgb_error.fit(X_val, val_error)

variance_hat_residuals = xgb_error.predict(X_test)


The relationship between variance and $$\sigma$$ and $$\nu$$ for a t-locationScale distribution is

var = $$\sigma^2 *\frac{\nu}{\nu-2}$$

Now here is where I make an assumption which I am not sure makes sense.

I assume that the degrees of freedom $$\nu$$ is the same as for all residuals, $$\nu = 2.56$$ and then I solve for $$\sigma$$ through the following.

$$\hat{\sigma} = \sqrt{\frac{\hat{var}*(\nu-2)}{\nu}}$$

And estimate the lower and upper quantiles from this distribution.

 residual_distribution_lower_quantile = scipy.stats.t.ppf(q = 0.025, df = 2.56, scale = sigma)
residual_distribution_upper_quantile = scipy.stats.t.ppf(q = 0.0975, df = 2.56, scale = sigma)


I then predict the lead time $$\hat{y}$$ and say that the mean of the distribution is $$\hat{y}$$

   pred = xgb.prediction(X_test)
lower_interval = pred + residual_distribution_lower_quantile
upper_interval = pred + residual_distribution_upper_quantile


Does it make sense to make the claim of $$\nu$$ is static? My score for the prediction interval is now $$81\%$$ since I am clearly simplifying the problem.

Any suggestions for improving my method?