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:
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)?
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)
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.