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