In Bayesian linear regression, if we want to get confidence intervals for predictions of a new observation. I was thinking of the following two options.

  1. Use the quantiles from samples sampled from the coefficents in the posterior $\beta | \mathbf{y}; \mathbf{X}$

  2. Use the predictive posterior distribution, $y^* \, | \beta, \mathbf{y}; \mathbf{x}^*, \mathbf{X}$

Where $\mathbf{x}^*$ is a new observation. What are the advantages of 2. over 1.

If we have the setup:

\begin{gather} \mathbf{y} \, | \beta ; \mathbf{X} \sim \mathcal{N}( \mathbf{X} \beta, \mathbf{I} \sigma^2) \\ \beta \sim \mathcal{N}(0, \mathbf{I}\sigma_{\beta}^2) \nonumber \end{gather}

Then we can get a closed form solution for the posterior of $\beta$

\begin{equation} \beta \, | \mathbf{y}; \mathbf{X} \sim \mathcal{N}\big((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X} \mathbf{y}^T , \; (\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2\big) \end{equation}

where $\lambda := \frac{\sigma^2}{\sigma_{\beta}^2}$.

Therefore can see that for a new point $\mathbf{x}^*$ the variance of $(\mathbf{x}^*)^T\beta$ will be

$$ (\mathbf{x}^*)^T((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2)\mathbf{x}^* $$

But, if we compute the predictive posterior:

\begin{align*} p(y^* \, | \beta, \mathbf{y}; \mathbf{x}^*, \mathbf{X}) = \int p(y^* \, | \beta; \mathbf{x}^*, \sigma^2) p(\beta \, | \mathbf{y}; \mathbf{X}, \sigma^2, \sigma_{\beta}^2) \; d\beta \end{align*}

we can compute this integral and get that the variance is:

$$ \sigma^2 + (\mathbf{x}^*)^T((\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I})^{-1} \sigma^2)\mathbf{x}^* $$

  • 1
    $\begingroup$ I think what you're asking for is a posterior predictive interval, isn't it? As in, some chosen interval of $p(y_{new}|\theta,y,x_{new}) = \int_{\Theta} p(y_{new}|\theta,x_{new}) p(\theta|y) d\theta$ where $\theta = [ \beta \, \sigma^2 ]^{\mathrm{T}}$. I don't see how your option 1 would work there. $\endgroup$
    – Durden
    Commented Jun 2, 2023 at 16:46
  • $\begingroup$ Ok thank you, so for example in Stan doing something like this to sample from the posterior and estimate quantiles. $\endgroup$
    – Dylan Dijk
    Commented Jun 8, 2023 at 11:09

2 Answers 2


Your question is more involved than it seems. Let me explain.

First, you have to decide on a tentative model for regression. Well, every Bayesian model involves a likelihood function and a prior probability. The likelihood function usually is rather well understood. Broadly speaking, the likelihood function describes how the data agree with the tentative model.

But what about the prior? The prior requires you to specify everything you know about the solution. For example: images cannot have negative pixel values, and physical processes cannot change abruptly. All of this is part of your professional knowledge. Or should be. And YOU have to decide on your prior!

Second, after multiplying the prior probability distribution and the likelihood function, you have obtained an unnormalized posterior probability distribution. But your model is still tentative. There is no alternative than to compute the normalization constant. This constant is sometimes called evidence or partition function.

Third, the next step is "model selection." The evidence is a numerical indication of the quality of the tentative model. But with the evidence values of several alternative models at hand, you can make a well-motivated selection of the best model. Only after you have selected the best model can you start with predictions and the corresponding error bars.

I have written a tutorial about this subject What is your model?, which describes the many pitfalls. You don't need to buy it as you can read 80% for free on Amazon. Good luck.


When using MCMC to sample the predictions you would first take a sample from the posterior for the parameters $\beta$, and then plug in the parameters and the data $\tilde{\mathbf{x}}$.

Let's use the Stan code example to illustrate it.

// [...]
model {
  y ~ poisson_log(alpha + beta * x);
  { alpha, beta } ~ normal(0, 1);
generated quantities {
  array[N_tilde] int<lower=0> y_tilde
    = poisson_log_rng(alpha + beta * x_tilde);

What the code does, is it first samples the alpha and beta parameters, then it predicts the mean of the Poisson distribution in alpha + beta * x_tilde, and for each x_tilde it samples one value from the Poisson distribution with the predicted mean.

So it does use the samples from the posterior distribution but the samples alone do not produce the posterior predictive samples. I'm not sure how you would like to use them for this purpose other than by plugging them into the model to sample the predictions.

  • $\begingroup$ I think they differ in that the predictive posterior will take into account information in the likelihood. So for example in Bayesian linear regression the predictive posterior will have an extra term in the variance, as opposed to just the variance in the posterior of the coefficients. $\endgroup$
    – Dylan Dijk
    Commented Jun 4, 2023 at 10:33
  • $\begingroup$ @DylanDijk likelihood does not have any “extra variance”. $\endgroup$
    – Tim
    Commented Jun 4, 2023 at 16:39
  • $\begingroup$ I have added a section at the bottom of my question now, to illustrate what I mean $\endgroup$
    – Dylan Dijk
    Commented Jun 8, 2023 at 11:00
  • $\begingroup$ This section in Stan documentation, I just found, also highlights what I mean, "there are two forms of uncertainty in posterior predictive quantities". $\endgroup$
    – Dylan Dijk
    Commented Jun 8, 2023 at 11:11
  • $\begingroup$ Where are you getting these quotes from? I cannot see any of my comments saying that. $\endgroup$
    – Dylan Dijk
    Commented Jun 8, 2023 at 11:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.