Bayesian inference: numerically sampling from the posterior predictive

Question

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?

Answer 1

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.

Answer 2

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
  }
}