# Why does sampling from the posterior predictive distribution $p(x_{new} \mid x_1, \ldots x_n)$ work without having to average out the integral?

In a Bayesian model, the posterior predictive distribution is usually written as:

$$p(x_{new} \mid x_1, \ldots x_n) = \int_{-\infty}^{\infty} p(x_{new}\mid \mu) \ p(\mu \mid x_1, \ldots x_n)d\mu$$

for a mean parameter $\mu$. Then, inside most books, such as this link:

Sampling MCMC

It is claimed that it is often easier to sample from $p(x_{new} \mid x_1, \ldots x_n)$ using Monte Carlo methods. Commonly, the algorithm is to:

for $j=1 \ldots J$:

1) Sample $\mu^{\ j}$ from $p(\mu \mid x_1, \ldots x_n)$ then

2) Sample $x^{\ * j}$ from $p(x_{new} \mid \mu^{\ j})$.

Then, $x^{\ * 1}, \ldots, x^{\ * J}$ will be an iid sample from $p(x_{new} \mid x_1, \ldots x_n)$.

What confuses me is the validity of this technique. My understanding is that Monte Carlo approaches will approximate the integral, so in this case, why do the $x^{\ * j}$'s each constitute a sample from $p(x_{new} \mid x_1, \ldots x_n)$?

Why isn't is the case that the average of all those samples instead will be distributed as $p(x_{new} \mid x_1, \ldots x_n)$? I am under the assumption that I am creating a finite partition to approximate the integral above. Am I missing something? Thanks!

What you are actually doing with the two-step process you've outlined is sampling from the joint distribution $p(x_{new}, \mu \thinspace | \thinspace x_1, \dots, x_n)$, then ignoring the sampled values of $\mu$. It's not altogether intuitive, but, by ignoring the sampled values of $\mu$, you are integrating over it.

A simple example may make this clear. Consider sampling from $p_X(x \thinspace | \thinspace y) = 1/y \thinspace \text{I}(0,y)$, uniform over $(0,y)$, and $p_Y(y) = 1$, uniform over $(0,1)$. You should be able to see, intuitively, what $\int_0^1p_X(x \thinspace | \thinspace y)p_Y(y)dy$ will look like. We construct some simple, horribly inefficient, R code (written this way for expository purposes) to generate the samples:

samples <- data.frame(y=rep(0,10000), x=rep(0,10000))
for (i in 1:nrow(samples)) {
samples$y[i] <- runif(1) samples$x[i] <- runif(1, 0, samples$y[i]) } hist(samples$x)


samples is clearly a random sample from the joint distribution of $x$ and $y$. We ignore the $y$ values and construct a histogram of only the $x$ values, which looks like: If you think carefully about it, you will see that the samples of $x$ do not depend upon any particular value of $y$. Instead, they depend (collectively) on a sample of values of $y$. This is why ignoring the $y$ values is equivalent to integrating out $y$, at least from a random number generation perspective.

On the other hand, consider what happens if you average. You'll get just one number from your Monte Carlo run, namely, the average of the $x_{new}$ samples. This isn't what you want (in your case)!

• Thanks for your post, do you know if there's a mathematically rigorous way to think about it? – user1398057 Dec 6 '15 at 17:40

I think you definetly have to mix over the sampled values eventually. There are also lecture notes by Peter Hoff on "Introduction to Bayesian Statistics for the Social Sciences" saying so. Otherwise you wouldn't have taken into account the masses recieved from the posterior. So, you build the empirical distribution of your samled values x^{*j} and then sample again from this distribution.

As an example: If your posterior was discrete (only point masses on atoms) then some of your parameter samples are going to take on the same values. If you finally mix over them, you take into account "how often" such parameter emerged from the posterior - put differently, how likely it is. Then averaging according to this appearances gives the posterior predictive which should appro. be the same as doing above procedure with the eventual mixing, at least if the sample size(s) is(are) large.

I think that the existing answers, which are very good, might be enhanced by an example with discrete random variables. We have $$p(x_{new} \mid x_1, \ldots x_n) = \int_{-\infty}^{\infty} p(x_{new},\mu \mid x_1, \ldots x_n)d\mu = \int_{-\infty}^{\infty} p(x_{new}\mid \mu) \ p(\mu \mid x_1, \ldots x_n)d\mu$$

To simplify, consider a $$\mu$$ that is binary: $$p(\mu = 1 \mid x_1 \dots x_n) = p$$ and $$p(\mu = 0 \mid x_1 \dots x_n) = 1-p$$. Suppose further that $$x_{new}$$ is binary with $$p(X=1)=\mu-1$$ and $$p(X=0)=\mu$$. I won't use these probabilities going forward, but you can see that $$x_{new}$$ depends on $$\mu$$.

Say we then draw 14 samples using $$\mu \sim p(\mu \mid x_1,\dots, x_n)$$ and $$x_{new} \sim p(x_{new} \mid \mu )$$. We get the following. As mentioned by @jbowman, we are actually sampling from $$p(x_{new}, \mu \mid x_1 \dots x_n)$$.

    mu    x_new
1.  1       0
2.  1       1
3.  0       0
4.  1       1
5.  0       0
6.  0       0
7.  0       0
8.  1       1
9.  1       1
10. 0       1
11. 1       0
12. 1       1
13. 0       1
14. 1       1


We can illustrate the fact that we are sampling from the joint $$p(x_{new}, \mu \mid x_1,\dots, x_n)$$ more explicitly by first constructing a table of counts.

        x_new
0    1
-----------

0    6    1
mu
1    2    5


Dividing each entry by the total (6 + 1 + 2 +5 = 14) gives

        x_new
0    1
-----------

0    0.43    0.07
mu
1    0.14    0.36


Which is the empirical joint distribution. Eg, our estimate of $$p(x_{new}=0, \mu=0)=.43$$. Hence our sampling procedure has given us the joint.

Finally, we will see why it is actually necessary to "evaluate" the integral (although not to average out the integral). This is implicit in @jbowman's answer when they said

It's not altogether intuitive, but, by ignoring the sampled values of 𝜇, you are integrating over it.

To obtain $$p(x_{new} \mid x_1 \dots x_n)$$, we simply sum over rows.

        x_new
0    1
-----------

.57  .43


This is what's implied by "ignoring the sampled values of $$\mu$$" and this is the marginalization step. Another way this is commonly done is by constructing a histogram (by summing over rows, we have kind of constructed a histogram here).

So, the sampling procedure does not give us the marginal - in other words, it doesn't "work" according to your definition in the question. Rather, it gives us the joint, and we commonly (by ignoring $$\mu$$, by constructing a histogram, or by getting quantiles) marginalize over $$\mu$$.