If the pp distribution is
p(yrep | y) = \int p(yrep | theta) p(theta | y) dtheta,
and we know just some (S) samples from the posterior then the pp distribution should be
p(yrep | y) = mean( p(yrep | theta_s) ).
If, for instance, the likelihood is gaussian, then the pp distribution is an average of S gaussians of parameters
theta_s. I would have made a random sample of, say, M numbers of every those gaussians to get a matrix of S rows and M columns, and then the pp distribution would be (with R):
colMeans(yrep) # A vector of M elements.
But in the many examples and documentation I read the sampling procedure is to sample M = 1 from the s-th likelihood and have the distribution of yrep as the vector of S samples.
What is wrong in my first assumptions for the distribution?