# Bayesian inference: numerically sampling from the posterior predictive

In the Bayesian approach the posterior distribution $P(\theta|x)$ of some hyper-parameter $\theta$ is given by:

$$P(\theta|x) = \frac{P(x|\theta)P(\theta)}{P(x)}$$

Where $x$ are samples of a random variable $X$. In many cases the posterior $P(\theta|x)$ can not be analytically computed and is instead calculated numerically using Markov Chain Monte Carlo. A popular approach is the Metropolis-Hastings algorithm which draws samples from:

$$P(\theta|x) \propto P(x|\theta)P(\theta)$$

The posterior predictive which describes the distribution of a new i.i.d. data point $x_\text{new}$ is given by:

$$P(x_\text{new}) = \int_\theta P(x_\text{new}|\theta)P(\theta|x)\, d\theta$$

My question is what are the popular methods to numerically sample from the posterior predictive $P(x_\text{new})$ if the posterior $P(\theta|x)$ is computed using Markov chain Monte Carlo?

If you can simulate values from $P(x_{\text{new}}|\theta)$, you can simply use your $N$ samples from posterior predictive and generate $x_{\text{new},i}$ for each posterior sample from this model to get a sample from the posterior predictive $\{x_{\text{new}}\}_{i=1}^N$.

This amounts to obtaining a collection $\{x_{\text{new},i}, \theta_i\}$ and discarding the value $\theta_i$, thus marginalizing over the vector of model parameters.

• @Ibelzile Good answer. Not sure whether systematic or random sampling of $\theta_i$ is best. I.e. do I cycle through the vector of $\theta$ or choose values of the vector at random. – egg Jan 30 '17 at 10:58

Here is an instantiated example of the answer provided by lbelzile. The application is linear regression, and the goal is to find the posterior predicted distribution of $y$ values at a probed $x$ value: http://doingbayesiandataanalysis.blogspot.com/2016/10/posterior-predictive-distribution-for.html The key idea is that at every step in the MCMC chain, use that step's parameter values to randomly generate a $y$ value from the model.

Edit in response to comment: Below is an extended excerpt from the blog post.

Suppose you've done a (robust) Bayesian multiple linear regression, and now you want the posterior distribution on the predicted value of $y$ for some probe value of $⟨x_1,x_2,x_3,...⟩$. That is, not the posterior distribution on the mean of the predicted value, but the posterior distribution on the predicted value itself. I showed how to do this for simple linear regression in a previous post; in this post I show how to do it for multiple linear regression. (A lot of commenters and emailers have asked me to do this.)

The basic idea is simple: At each step in the MCMC chain, use the parameter values to randomly generate a simulated datum $y$ at the probed value of $x$. Then examine the resulting distribution of simulated $y$ values; that is the posterior distribution of the predicted $y$ values.

To implement the idea, the first programming choice is whether to simulate the $y$ value with JAGS (or Stan or whatever) while it is generating the MCMC chain, or to simulate the $y$ value after the MCMC chain has previously been generated. There are pros and cons of each option. Generating the value by JAGS has the benefit of keeping the code that generates the $y$ value close to the code that expresses the model, so there is less chance of mistakenly simulating data by a different model than is being fit to the data. On the other hand, this method requires us to pre-specify all the $x$ values we want to probe. If you want to choose the probed $x$ values after JAGS has already generated the MCMC chain, then you'll need to re-express the model outside of JAGS, in R, and run the risk of mistakenly expressing it differently (e.g., using precision instead of standard deviation, or thinking that y=rt(...) in R will use the same syntax as y~dt(...) in JAGS). I will show an implementation in which JAGS simulated the $y$ values while generating the MCMC chain.

[... example in original post not copied here ...]

The Jags model specification looks like the following. Notice at the very end the randomly generated $y$ values, denoted yP[i].

# Standardize the data:
data {
ym <- mean(y)
ysd <- sd(y)
for ( i in 1:Ntotal ) {
zy[i] <- ( y[i] - ym ) / ysd
}
for ( j in 1:Nx ) {
xm[j]  <- mean(x[,j])
xsd[j] <-   sd(x[,j])
for ( i in 1:Ntotal ) {
zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
}
# standardize the probe values:
for ( i in 1:Nprobe ) {
zxProbe[i,j] <- ( xProbe[i,j] - xm[j] ) / xsd[j]
}
}
}
# Specify the model for standardized data:
model {
for ( i in 1:Ntotal ) {
zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )
}
# Priors vague on standardized scale:
zbeta0 ~ dnorm( 0 , 1/2^2 )
for ( j in 1:Nx ) {
zbeta[j] ~ dnorm( 0 , 1/2^2 )
}
zsigma ~ dunif( 1.0E-5 , 1.0E+1 )
nu ~ dexp(1/30.0)
# Transform to original scale:
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
beta0 <- zbeta0*ysd  + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
sigma <- zsigma*ysd
# Predicted y values at xProbe:
for ( i in 1:Nprobe ) {
zyP[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zxProbe[i,1:Nx] ) ,
1/zsigma^2 , nu )
yP[i] <- zyP[i] * ysd + ym
}
}

• Thanks for your answer John (it's great to see your answers in here!), but we want out site to be self-contained and not to rely entirely on external links. Could you possibly summarize the link in your answer, or alternatively copy-and-paste some part of it for reference? – Tim Jan 27 '17 at 19:38
• Thanks, John: the answer is provided by Ibelzile, I only edited it. – Xi'an Jan 29 '17 at 10:46
• This is the usual answer, but I began to doubt it. The MCMC sample isn't iid: it's heavily autocorrelated. Even if the initial samples are discarded, so that I retain only samples from the (more or less) stationary part, still sample $\theta^{(i)}$ and $\theta^{(i+1)}$ are strongly correlated, no matter how large $i$ is. It's a Markov chain after all. So I think one should thin the chain before collecting the $y$ sample, or we risk having a non-iid sample, and thus for example the mean of $y$ computed with this sample won't be as accurate as we would expect. – DeltaIV Jan 29 '17 at 14:04
• There is a general result that tells you that thinning does not reduce the variance of your estimate – Xi'an Jan 30 '17 at 5:35
• @Xi'an ah, this is very interesting, thanks. So thinning does reduce the autocorrelation of the sample, but not its variance? I'm surprised: I would expect a more correlated sample to have higher variance, since consecutive errors "go together" instead than randomly compensating. Does your book with George Casella discuss this point? PS I would also like to ask you another question about your book. May we talk in chat? Otherwise I can send you a mail. – DeltaIV Jan 30 '17 at 22:59