The posterior predictive density for replicated data $y_{rep}$ given observed data $y$ is $$ p(y_{rep} \mid y) = \int p(y_{rep} \mid \theta) \cdot p(\theta \mid y) \, \textrm{d}\theta, $$

where the parameter vector $\theta$ is marginalised out from the joint distribution $p(y_{rep}, \theta \mid y)$ (ref).

In graphical posterior predictive checks, $B$ draws are sampled from the posterior distribution of $\theta$ and $B$ outcome vectors ($y_{rep}^{(1)}, \ldots, y_{rep}^{(B)}$) of the same size of the analytic dataset are simulated from the model (one vector for each posterior draw).

Kernel density estimate (for example) of the observed outcome is then overlaid to $B$ KDEs, one for each of the $B$ simulated outcome vectors (see, for example, Figure 6 in Gabry et al. on arXiv).

My question is: why did they write that each of the $B$ outcome vectors ($y_{rep}^{(b)}$) is simulated from the posterior predictive distribution? It seems to me that $y_{rep}^{(b)}$ is simulated conditional on a specific value of $\theta$ from the posterior distribution. There's no marginalisation there (afaict), and therefore how can each given outcome vector $y_{rep}^{(b)}$ come from the posterior predictive distribution?

To obtain a random draw from the posterior predictive distribution I would concatenate the $(y_{rep}^{(1)}, \ldots, y_{rep}^{(B)})$ vectors, marginalising out the posterior distribution (see for example this question here on Cross Validated). At that point, I would overlay the KDE of the observed outcome with the KDE of (what I believe is) the posterior predictive distribution.


1 Answer 1


To get the posterior predictive distribution, you would need to calculate the integral

$$ \int p(y_{rep} \mid \theta) \cdot p(\theta \mid y) \, \textrm{d}\theta $$

Notice that the code snipped from your first link does, it draws samples for the particular values of $\theta = (\alpha, \beta, \sigma)$. The example uses MCMC for Monte Carlo integration, i.e. instead of calculating the integral directly, you sample from the distribution and average over the samples.

generated quantities {
  array[N] real y_rep = normal_rng(alpha + beta * x, sigma);

You need many samples because you approximate the distribution with the empirical distribution of the Monte Carlo samples.

  • $\begingroup$ Thank you, Tim! From that code, one gets a $S \times N$ matrix, where $S$ is the number of posterior draws and $N$ the number of observations. In Figure 6, they take 100 out of the $S$ rows and plot the corresponding 100 densities (the coloured lines). Each density is therefore based on $N$ observations simulated from the model using the same parameter vector drawn from the posterior. You wrote "You need many samples..." and I agree, but I would plot only one density from combining all the $S \times N$ sampled values, not by row as in Figure 6. What am I missing? $\endgroup$
    – boscovich
    Feb 10 at 13:31
  • $\begingroup$ @boscovich usually you are interested also in the variability of the predictions. $\endgroup$
    – Tim
    Feb 10 at 13:48
  • $\begingroup$ Thanks again. I don't want to take up any more of your time, but how each of the $S$ vectors of length $N$ is drawn from the PPD is still beyond me. If one wanted to assess variability, I would generate multiple $S \times N$ matrices (say, $M$ of such matrices) and plot the associated $M$ densities. But then again, it's clear I am missing something. $\endgroup$
    – boscovich
    Feb 10 at 13:58
  • $\begingroup$ @boscovich each of them is drawn from $p(y_{rep}|\theta)$, we use Monte Carlo integration to get the posterior predictive distribution. But the $p(y_{rep}|\theta)$ draws are interesting by themselves, so we often look at them as well. $\endgroup$
    – Tim
    Feb 10 at 14:39
  • $\begingroup$ Thanks! "Each of them is drawn..." now this makes a lot of sense to me! But calling the $p(y_{rep} | \theta)$ "simulations from the PPD ($p(y_{rep} \mid y)$)" is common. Other examples are: Stata's Bayesian analysis reference manual (page 444, stata.com/manuals/bayes.pdf) or Bayesian Data Analysis by Gelman et al. (Figure 6.2, stat.columbia.edu/~gelman/bayescomputation/bdachapter6.pdf). Personally, I find that quite confusing. $\endgroup$
    – boscovich
    Feb 10 at 15:08

Your Answer

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

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