Understanding Bayesian Predictive Distributions I'm taking an Intro to Bayes course and I'm having some difficulty understanding predictive distributions. I understand why they are useful and I'm familiar with the definition, but there are some things I don't quite understand. 
1) How to get the right predictive distribution for a vector of new observations
Suppose that we have built a sampling model $p(y_i | \theta)$ for the data and a prior $p(\theta)$. Assume that the observations $y_i$ are conditionally independent given $\theta$. 
We have observed some data $\mathcal{D} = \{y_1, y_2, \, ... \, , y_k\}$, and we update our prior $p(\theta)$ to the posterior $p(\theta | \mathcal{D})$. 
If we wanted to predict a vector of new observations $\mathcal{N} = \{\tilde{y}_1, \tilde{y}_2, \, ... \, , \tilde{y}_n\}$, I think we should try to get the posterior predictive using this formula
$$
p(\mathcal{N} | \mathcal{D}) = \int p(\theta | \mathcal{D}) p ( \mathcal{N} | \theta) \, \mathrm{d} \theta = \int p(\theta | \mathcal{D}) \prod_{i=1}^n p(\tilde{y}_i | \theta) \, \mathrm{d} \theta,
$$
which is not equal to
$$
\prod_{i=1}^n  \int p(\theta | \mathcal{D}) p(\tilde{y}_i | \theta) \, \mathrm{d} \theta,
$$
so the predicted observations are not independent, right?
Say that $\theta | \mathcal{D} \sim$ Beta($a,b$) and $p(y_i | \theta) \sim$ Binomial($n, \theta$) for a fixed $n$. In this case, if I wanted to simulate 6 new $\tilde{y}$, if I understand this correctly, it would be wrong to simulate 6 draws independently from the Beta-Binomial distribution that corresponds to the posterior predictive for a single observation. Is this correct? I don't know how to interpret that the observations are not independent marginally, and I'm not sure I understand this correctly.
Simulating from posterior predictives
Many times when we simulate data from the posterior predictive we follow this scheme:
For $b$ from 1 to $B$:
1) Sample $\theta^{(b)}$ from $p(\theta | \mathcal{D})$.
2) Then simulate new data $\mathcal{N}^{(b)}$ from $p(\mathcal{N} | \theta^{(b)})$. 
I don't quite know how to prove this scheme works, although it looks intuitive. Also, does this have a name? I tried to look up a justification and I tried different names, but I had no luck.
Thanks!
 A: Suppose that $X_1,\dots,X_n,X_{n+1}$ are conditionally independent given that $\Theta=\theta$. Then,
$$
  f_{X_{n+1}\mid X_1,\dots,X_n}(x_{n+1}\mid x_1,\dots,x_n) = \int f_{X_{n+1},\Theta\mid X_1,\dots,X_n}(x_{n+1},\theta\mid x_1,\dots,x_n)\,d\theta
$$
$$
= \int f_{X_{n+1}\mid\Theta,X_1,\dots,X_n}(x_{n+1}\mid\theta,x_1,\dots,x_n) f_{\Theta\mid X_1,\dots,X_n}(\theta\mid x_1,\dots,x_n) \, d\theta
$$
$$
= \int f_{X_{n+1}\mid\Theta}(x_{n+1}\mid\theta) f_{\Theta\mid X_1,\dots,X_n}(\theta\mid x_1,\dots,x_n) \, d\theta \, ,
$$
in which the first equality follows from the law of total probability, the second follows from the product rule, and the third from the assumed conditional independence: given the value of $\Theta$, we don't need the values of $X_1,\dots,X_n$ to determine the distribution of $X_{n+1}$.
The simulation scheme is correct: for $i=1,\dots,N$, draw $\theta^{(i)}$ from the distribution of $\Theta\mid X_1=x_1,\dots,X_n=x_n$, then draw $x_{n+1}^{(i)}$ from the distribution of $X_{n+1}\mid\Theta=\theta^{(i)}$. This gives you a sample $\{x_{n+1}^{(i)}\}_{i=1}^N$ from the distribution of $X_{n+1}\mid X_1=x_1,\dots,X_n=x_n$.
A: I'll try to go over the intuition behind generating the posterior predictive distribution step-by-step.
Let $y$ be a vector of observed data that come from a probability distribution $p(y|\theta)$ and let $ \tilde y$ be a vector of future (or out-of-sample) values we want to predict. We assume that $ \tilde y$ comes from the same distribution as $ y$. It might be tempting to use our best estimate of $\theta$---such as the MLE or MAP estimate---to obtain information about this distribution. However, doing so would inevitably ignore our uncertainty about $ \theta$. Thus, the appropriate way to procede is to average over the posterior distribution of $ \theta$, namely $ p(\theta|y)$. Notice also that $ \tilde y$ is independent of $ y$ given $ \theta$, as it is assumed to be an independent sample drawn from the same distribution as $ y$. Thus,

$ \displaystyle
p(\tilde y| \theta, y) = \frac{p(\tilde y, y|\theta )p(\theta)}{p(\theta, y)} =
\frac{p(\tilde y|\theta )p(y |\theta) p(\theta)}{p(y| \theta)p(\theta)} = p(\tilde y |\theta).$

The posterior predictive distribution of $ \tilde y$ is thus,

$ \displaystyle p(\tilde y|y ) = \int_\Theta p(\tilde y | \theta,y) p(\theta | y) d\theta = \int_\Theta p(\tilde y | \theta) p(\theta | y) d\theta $

where $ \Theta$ is the support of $ \theta$.
Now, how do we obtain the samples from $ p(\tilde y|y)$? The method you describe is sometimes called the method of composition, which works as follows:

for s = 1,2,...,S do 
draw $\theta^{(s)}$ from $ p(\theta|y)$
draw $\tilde y^{(s)}$ from $ p(\tilde y|\theta^{(s)})$

where, in most situations, we have already the draws from $ p(\theta|y)$, so that only the second step is required.
The reason why this works is quite simple: First note that $ p(\tilde y, \theta | y) = p(\tilde y| \theta, y)p(\theta | y)$. Thus, sampling a parameter vector $\theta^{(s)}$ from $ p(\theta|y)$ and, then, using this vector to sample $ \tilde y^{(s)}$ from $ p(\tilde y | \theta^{(s)}) = p(\tilde y | \theta^{(s)}, y)$ yields samples from the joint distribution $ p(\tilde y, \theta|y)$. It follows that, the sampled values $\tilde y^{(s)}, s=1,2,...,S$ are samples from the marginal distribution, $p(\tilde y|y)$. 
A: To address your first question: yes, the observations are not independent if you don't know the value of $\theta$. Say, you've observed that $\tilde{y}_1$ has rather extreme value. It may be an indication that the unknown value of $\theta$ itself is extreme, and, thus, you should expect other observations to be extreme as well.
