Computing Bayesian Credible Intervals for Bayesian Regression

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.

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!

• Thank you very much for your answer. Could you recommend some literature that extensively explains and illustrates this concept? – Sean Jun 20 '17 at 19:43
• Hi Sean, depending on your prior knowledge I would recommend Statistical Rethinking by Richard McElreath, which is a great resource for people of all levels who are interested in Bayesian analysis, or Bayesian Data Analysis by Andrew Gelman (and many others ) which is great but might not be accessible if you have no affinity at all for mathematics. The latter will give you more in terms of formulas for this. Google should help you find them, otherwise I will provide links when I'm not on my phone. – Maurits M Jun 20 '17 at 20:20
• As promised, some links. You can have a look at en.wikipedia.org/wiki/Posterior_predictive_distribution In addition, here are the webpages of the books I mentioned: stat.columbia.edu/~gelman/book and xcelab.net/rm/statistical-rethinking Finally, I mentioned Stan before; there are several other similar tools (JAGS, BUGS, ...) that can also help. Depending on the type of model that you want to fit you may not even need such specialized tools. – Maurits M Jun 21 '17 at 7:31
• Thank you very much Maurits. I have started reading the book Bayesian Data Analysis and find it very comprehensible thus far. – Sean Jun 22 '17 at 7:54
• You're welcome! In BDA3, the chapters on Bayesian linear regression will give you a lot of pointers related to your question. Good luck! – Maurits M Jun 23 '17 at 8:50