Computing Bayesian Credible Intervals for Bayesian Regression

Question

I understand the Bayesian regression approach, in which I compute the posterior distribution of the parameters $\theta_i$. Now I am interested in how to compute the credible intervals of the regression curve: $$y=x_0+\theta_1*x_1+..+\theta_n*x_n$$

In a frequentist setting, a prediction interval of a linear regression prediction for $x_{new}$ states: $$\hat{\mu}±t_{(α/2,n-p)}*s*\sqrt{1+x_{new}(X^TX)^{-1}x^\top_{new}}$$

yet for the Bayesian regression, I am unsure how to obtain the credible interval of the regression curve from the credible intervals of the parameters $\theta_i$. What is the correct formula, and how does it come about? Thank you for your help in advance!

EDIT:

I have found out that the predictive distribution for Bayesian regression with a zero mean Gaussian prior with covariance matrix $\Sigma_p$, i.e.

$$w\sim N(0,\Sigma_p)$$

comes out to be:

$$\int p(f_*|x_*,w)p(w|X,y)dw=\int x_*^\top *w* p(w|X,y)dw$$

and this is distributed normally: $$N\sim (\frac{1}{\sigma_n^2}x_*^{\top}A^{-1}Xy , x_*^{\top}A^{-1}x_*)$$

where $A=\sigma^{-2}_nXX^\top+\Sigma_p^{-1}$

I have thus computed the standard deviation $std=x_*^{\top}A^{-1}x_*$ for a range of values of $x_*$ to plot the upper and lower prediction boundaries:

The polynomial regression is unregularized (as I do not know how to compute prediction intervals for the regularized case). I am left with the questions:

1. The posterior predictive distribution predicts an $f_*$. Am I right in the assumption, that the credible interval obtained from this distribution is an interval surrounding the mean prediction, and not a single prediction (thus not a prediction interval)?

2. If this is correct: To get the (Bayesian Regression) plot, I have added the MSE to the computed standard deviations, to account for the error term $\epsilon$. Is this a valid approach to get from a mean prediction interval to a single prediction interval? (It is my intuitive understanding, that the $\epsilon_s$ cancel out in the case of a mean prediction)

3. Lastly, to compute the Bayesian Regression curve, I have used the Python library Scikit-Learn (http://scikit-learn.org/stable/modules/linear_model.html#bayesian-ridge-regression). However, this library assumes a NIG-prior (which is, to my best understanding, conjugate). I have looked through a lot of literature, but as I am just beginning to understand this topic, I simply could not find out: What is the predictive distribution in this case of a (hierarchical) NIG prior? I would like to use it to compute prediction intervals, as I attempted to in the plots above.

Answer 1

Unfortunately, the interval that you are looking for is not uniquely determined. Essentially, what you need is the Posterior Predictive Density (PPD, see https://en.wikipedia.org/wiki/Posterior_predictive_distribution), which is the density function of new/unseen data given the observed data. This PPD depends on the posterior distribution of the parameters, $\theta_1, ..., \theta_n$ in your case. It can be written as

$p(y^* | y, x, x^*) = \int p(y^*, \theta | y, x, x^*) d\theta = \int p(y^* | \theta) p (\theta | y, x, x^*)d\theta$

where $y^*$ represents the unseen response data, $y$ represents the known response data, $x$ and $x^*$ represent the predictor values that correspond to $y$ and $y^*$, and $\theta$ represents the parameters. The last factor in the final integral is the posterior distribution of $\theta$ given $y, x, x^*$. As the PPD depends on the posterior distribution of the parameters, this, in turn, depends on the prior distribution of the parameters (and on the chosen data model / likelihood function). This means that for each prior you may choose, your posterior distribution changes (and, as a result, your interval as well).

Usually, when choosing completely uninformative (i.e. flat) priors, along with a Normal likelihood for the response values given the predictors, the results of a Bayesian analysis overlap with those of a frequentist analysis. Then again, flat priors are usually a poor choice for such a model.

When you know which priors you want to use for your analysis, it may be possible to compute the PPD analytically, but in many cases this is simply impossible. I'd recommend using a tool like Stan (http://mc-stan.org) to draw samples from the posterior distribution and then use those to determine a credible interval for your parameters and your new (simulated) data.

Hope this helps!