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?