I'm confused on how to evaluate the posterior predictive distribution for Bayesian linear regression, past the basic case described here on page 3, and copied below.

$$ p(\tilde y \mid y) = \int p(\tilde y \mid \beta, \sigma^2) p(\beta, \sigma^2 \mid y) $$

The basic case is this linear regression model:

$$ y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \sigma^2)$$

If we use either a uniform prior on $\beta$, with a scale-Inv$\chi^2$ prior on $\sigma^2$, OR the normal-inverse-gamma prior (see here) the posterior predictive distribution is analytic and is student t.

What about for this model? $$ y = X \beta + \epsilon, \hspace{10mm} y \sim N(X \beta, \Sigma)$$

When $y \sim N(X\beta, \Sigma)$, but $\Sigma$ is known, the posterior predictive distribution is multivariate Gaussian. Usually, you don't know $\Sigma$, but have to estimate it. Maybe you say its diagonal and make the diagonal a function of the covariates in some way. This is discussed in the linear regression chapter of Gelman's Bayesian Data Analysis.

Is there an analytic form for the posterior predictive distribution in this case? Can I just plug my estimate of it into a multivariate student t? If you estimate more than one variance, is the distribution still multivariate student t?

I am asking because say I have some $\tilde y$ already on hand. I want to know whether it is more likely to have been predicted by e.g. linear regression A, linear regression B

  • 1
    $\begingroup$ If you have posterior samples from the posterior distribution you can evaluate the predictive distribution with a Monte Carlo approximation $\endgroup$
    – niandra82
    Apr 10, 2015 at 21:04
  • $\begingroup$ Ah thanks, I could always do that. Is there no analytic formula in this case? $\endgroup$
    – bill_e
    Apr 10, 2015 at 21:19
  • $\begingroup$ The links are broken, by the way. Would be great if you were to incorporate the references in another way. $\endgroup$
    – Maxim.K
    Sep 21, 2018 at 13:51

2 Answers 2


If you assume a uniform prior on $\beta$, then the posterior for $\beta$ is $$ \beta|y \sim N(\hat{\beta},V_\beta). $$ with $$\hat{\beta} = [X'\Sigma^{-1}X]X'y \quad \mbox{and} \quad V_\beta =[X'\Sigma^{-1}X]^{-1}.$$ To find the predictive distribution, we need more information. If $\tilde{y} \sim N(\tilde{X}\beta,\tilde{\Sigma})$ and is conditionally independent of $y$ given $\beta$, then $$ \tilde{y}|y \sim N(\tilde{X}\hat{\beta}, \tilde{\Sigma} + V_\beta).$$ But typically for these types of models, $y$ and $\tilde{y}$ are not conditionally independent, instead, we typically have $$ \left( \begin{array}{c} y \\ \tilde{y} \end{array} \right) \sim N\left( \left[\begin{array}{c} X\beta \\ \tilde{X}\beta \end{array}\right], \left[ \begin{array}{cc} \Sigma & \Sigma_{12} \\ \Sigma_{21} & \tilde{\Sigma} \end{array} \right]\right). $$ If this is the case, then $$ \tilde{y}|y \sim N(\tilde{X}\hat{\beta} + \Sigma_{21}\Sigma^{-1}(y-X\hat{\beta}), \tilde{\Sigma} - \Sigma_{21}\Sigma^{-1}\Sigma_{12}). $$ This assumes $\Sigma, \Sigma_{12},$ and $\tilde{\Sigma}$ are all known. As you point out, typically they are unknown and need to be estimated. For the common models that have this structure, e.g. time series and spatial models, there generally won't be a closed form for the predictive distribution.


Under non-informative or multivariate Normal-Wishart priors, you have the analytical form as a multivariate Student's distribution, for a classical mutlivariate, multiple regression. I guess the developments in this document are related to your question (you may like the Appendix A :-) ). I typically compared the outcome with a posterior predictive distribution obtained using WinBUGS and the analytical form: they are exactly equivalent. The problem only becomes difficult when you have additional random effects in mixed-effect models, especially in unbalanced design.

In general, with classical regressions, y and ỹ are conditionally independent (residuals are i.i.d) ! Of course if it is not the case, then the proposed solution here is not correct.

In R,(here, solution for uniform priors), assuming you made a lm model (named "model") of one of the response in your model, and called it "model", here is how to obtain the multivariate predictive distribution

Y = as.matrix(datas[,c("resp1","resp2","resp3")])
X =  model.matrix(delete.response(terms(model)), 
           data, model$contrasts)
XprimeX  = t(X) %*% X
XprimeXinv = solve(xprimex)
hatB =  xprimexinv %*% t(X) %*% Y
A = t(Y - X%*%hatB)%*% (Y-X%*%hatB)
F = ncol(X)
M = ncol(Y)
N = nrow(Y)
nu= N-(M+F)+1 #nu must be positive
C_1 =  c(1  + x0 %*% xprimexinv %*% t(x0)) #for a prediction of the factor setting x0 (a vector of size F=ncol(X))
varY = A/(nu) 
postmean = x0 %*% hatB
nsim = 2000
ysim = rmvt(n=nsim,delta=postmux0,C_1*varY,df=nu) 

Now, quantiles of ysim are beta-expectation tolerance intervals from the predictive distribution, you can of course directly use the sampled distribution to do whatever you want.

  • 1
    $\begingroup$ Unfortunately, the document is not accesible anymore. Could you give the title of the publication or provide a new link? $\endgroup$
    – Jannis
    Aug 10, 2023 at 14:03

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.