How to get a confidence interval on linear regression at prediction time where parameters are distributions? I have a simple linear model $y = \beta+ \theta_1x_1 + \theta_2x_2$ where I have obtained the parameters through a bayesian MC approach, and so $\beta, \theta_1, \theta_2$ are all distributions (not fixed point estimates).
I have two questions here; at predict time (predicting test data the model has not been trained on)

*

*How should we compute the estimate for $y$ given the parameters are distributions with different means and variances? The prediction for $y$ should also be a distribution.


*When predicting, how should we obtain a confidence interval using the properties of the parameter distributions?
 A: If I understand correctly, what you are looking for is posterior predictive distribution of your test data. Let $\mathbf y_{tr}$ be your training set data and $\mathbf y_{ts}$ be your test data.
You need:
$$p(\mathbf y_{ts} \mid \mathbf y_{tr}) = \int_\theta p(\mathbf y_{ts}\mid \theta)p(\theta\mid \mathbf y_{tr})d\theta$$
where $\theta = (\beta, \theta_1, \theta_2)'$.
Here we are basically marginalizing over posterior density of coefficients. How to do this practically:
You already have samples drawn from $p(\theta\mid \mathbf y_{tr})$ - say of size $n$. Further, you must have assumed a family for $p(y\mid \theta)$. Let's say $y \mid \theta \sim N(X\theta,\sigma^2)$ with known $\sigma$ (for simplicity).
All you need to do is draw one sample from a normal distribution with each of the values of each value of $\theta$ for each of the data point in test data.
This is usually straight forward in software (for example, in R it would be just rnorm(n, mean = theta %*% t(X_test[1,]), sd = sig) - here theta is a $n \times 3$ matrix of coefficient draws and X_test is a $m \times 3$ test data of explanatory variables having $m$ observations). From this you will get $n$ values for each observation in your test data. Therefore you have samples drawn from $p(y_i \mid \mathbf y_{tr})$ for each of $i=1, \dots,m$.
Since for each $y_i$ you have a full sample, you can easily construct a confidence interval by finding quantiles.
