Posterior predictive distribution and posterior predictive checks

Question

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.

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.