Intuition behind posterior predictive distribution I've recently encountered the "posterior predictive distribution" $$p(\bar{x}|X)=E_\theta[p(\bar{x}|\theta)]=\int_\theta p(\bar{x}|\theta)\hspace{0.5mm}p(\theta|X)d\theta$$
where $\bar{x}$ is a new point, $\theta$ is the (vector or 1-D) parameters of the distribution and $X$ is the already observed sample.
However I'm not sure I understand where this formula comes from. I know that we don't know the true value of $\theta$. I think that since we don't know $\theta$, we say "why settle for one estimate of it, and not scan all it's possible values?" Is that correct? Even if it is, I'm not certain that I understand the expected value part. 
Sorry if my thoughts are a bit jumbled. Any insight would be appreciated, thanks.
 A: Let $X$ denotes the observations and $\theta \in \Theta$ the parameter. In a Bayesian approach, both are considered random quantities. 
The first step of modeling is to define a statistical model, i.e. the distribution of $X$ given $\theta$, which can be written as $X \mid \theta \sim p(\cdot \mid \theta)$. This is mainly done by expliciting a likelihood function.
Thus our statistical model describe the conditional distribution of $X$ given $\theta$.
From a Bayesian perspective, we also define a prior distribution for $\theta$ on $\Theta$: $\theta \sim \pi(\theta)$.
The prior predictive distribution
Before observing any data, what we have is simply the chosen model, $p(x \mid \theta)$, and the prior distribution of $\theta$, $\pi(\theta)$. 
One can then ask to see what is the marginal distribution of $X$, that is, the distribution of $X \mid \theta $ averaged over all possible values of $\theta$.
This can be simply written using expectation:
\begin{align*}
p(x) &= \mathbb{E}_\theta \Big [ p(x \mid \theta)  \Big ] \\
&= \int_\Theta p(x \mid \theta) \pi(\theta) d\theta.
\end{align*}
The posterior predictive distribution
The interpretation is the same than for the prior predictive distribution, is it the marginal distribution of $X \mid \theta$ averaged over all values of $\theta$.
But this time the "weighting" function to be used is not $\pi(\theta)$ but our updated knowledge about $\theta$ after observing data $X^*$: $\pi(\theta \mid X^*)$.
Using the known known Bayes theorem we have:
$$
\pi(\theta \mid X^*) = \frac{p(X^* \mid \theta) \pi(\theta)}{p(X^*)}
$$
And thus, the marginal distribution of $X \mid (X^*,\theta)$ averaged over $\Theta$ is:
$$
p(x \mid X^*) = \int_\Theta p(x \mid \theta) \pi(\theta \mid X^*)d\theta
$$
Example: Gamma-Poisson mixture.
Suppose our observations are made of counts, $X$, and we define a Poisson model that is: $X \mid \lambda \sim \mathcal{P}(\lambda)$.
From a Bayesian perspective, we also define a prior distribution for $\lambda$.
For mathematical reasons, it is appealing to use a Gamma distribution, $\lambda \sim \mathcal{G}(a,b)$.    
The prior predictive distribution
One particulariy of this Gamma-Poisson mixture, is that the marginal distribution will be distribution as a Negative-Binomial random variable.
That is, if  $X \mid \lambda \sim \mathcal{P}(\lambda)$ and $\lambda \sim \mathcal{G}(a,b)$ then, $X \sim \mathcal{NB}\big (a,\frac{b}{b+1} \big )$.
Thus the prior predictive distribution of $X$ is a Negative Binomial distribution $\mathcal{NB}\big (a,\frac{b}{b+1} \big )$.
The posterior predictive distribution
Now, say we have observed $n$ counts $X =(X_1,\dots,X_n)$.
First, thanks to the choice of a Gamma prior for $\lambda$, the posterior distribution of $\lambda$ can be easily derived as being also a Gamma distribution:
$$
\lambda \mid X \sim \mathcal{G} \bigg ( a  + \sum_{i=1}^n X_i , b+n \bigg)
$$
From what we saw for the prior predictive distribution, the posterior predictive distribution of $X$ will also be a Negative-Binomial:
$$
\mathcal{NB} \bigg ( a  + \sum_{i=1}^n X_i, \frac{b+n}{b+1+n} \bigg )
$$
Here is an example where $a=100$, $b=2$ and we observe the vector of counts $X=(85,80,70,65,71,92)$ :

Here is the R code to produce the plot :
### Gamma-Poisson mixture: prior and posterior predictive distributions :
require(ggplot2)
# Parameters of the prior distribution of lambda
a<-100
b<-2
x<-0:150
y1<-dnbinom(x,a,b/(b+1))


# Vector of observations and posterior predictive distribution
X<-c(85,80,70,65,71,92)
n<-length(X)
XS<-sum(X)
x<-0:150
au<-a+XS
bu<-b+n
y2<-dnbinom(x,size=au,prob=bu/(bu+1))


plot1<-ggplot() + aes(x=x,y=y1,colour="Prior") + geom_line(size=1)+
       geom_line(aes(x=x,y=y2,colour="Post"))+
       scale_colour_manual(breaks=c("Prior",    "Post"),
       values=c("#cd7118","#1874cd"),labels=c("Prior Predictive",
       "Postererior    Predictive"))+
       ggtitle("Prior and posterior predictive distributions for a=100 and b=2")+
       labs(x="x", y="Density") +
       theme(
           panel.background = element_rect(fill = "white",
                                colour = "white",
                                size = 0.5, linetype = "solid"),

           axis.line = element_line(size = 0.2, linetype = 'solid',
                         colour = "black"), 
           axis.text = element_text(size=10),
           axis.title = element_text(size=10),
           legend.title = element_blank(),
           legend.background = element_blank(),
           legend.key = element_blank(),
           legend.position = c(.7,.5)
         )
 plot1

