I observe independent, Poisson-distributed data $ D = \{x_1, ... x_n \} $ with mean parameter $ \mu $, i.e., $$x_i\stackrel{\text{iid}}{\sim}\mathcal{P}(\mu)$$ Over $ \mu $ I assume $ Gamma(\alpha_0, \beta_0) $ as a prior (where $ \alpha_0 $ and $ \beta_0 $ are some given hyperparameters). The posterior distribution is $ \mu|D \sim Gamma(a_0 + \sum\limits_{n=1}^N x_n), \beta_0 + N) $.$$\mu|D \sim Gamma\left(a_0 + \sum\limits_{n=1}^N x_n, \beta_0 + N\right) $$ I am interested in the mean and variance of the predictive distribution $ x_{N+1} | D $$p(x_{N+1} | D)$.
I know it's possible to show that the predictive distribution is a negative binomial, but, since I am only interested in mean and variance there should be a much quicker way by using identities, the law of total expectation and the law of total variance. However, I get confused when trying to apply them here.
The posterior predictive distribution for $ x_{N+1} $ is (leaving out conditioning on hyperparameters):
$ p(x_{N+1}|D) = \int_{\mu} p(x_{N+1}|\mu) \, p(\mu|D) \operatorname{d}\!\mu $ $$ p(x_{N+1}|D) = \int_{\mu} p(x_{N+1}|\mu) \, p(\mu|D) \operatorname{d}\!\mu $$
That is equivalent to the expectation:
$ p(x_{N+1}|D) = E_{\mu|D}\Big[p(x_{N+1}|\mu)\Big] $$$ p(x_{N+1}|D) = \mathbb{E}_{\mu|D}\Big[p(x_{N+1}|\mu)\Big] $$
The mean of the posterior distribution then is:
$ E_{x_{N+1}|D} \Big[E_{\mu|D}\Big[p(x_{N+1}|\mu)\Big]\Big]$$$ \mathbb{E}_{x_{N+1}|D} \Big[\mathbb{E}_{\mu|D}\Big[p(x_{N+1}|\mu)\Big]\Big]$$
Whichwhich is an iterated expectation to which I would like to apply the law of total expectation.
The law states: $ E_X[X] = E_Y[E_{X|Y}[X|Y]] $.$$ \mathbb{E}_X[X] = \mathbb{E}_Y[\mathbb{E}_{X|Y}[X|Y]] $$
Looking at the indices however, I don't know how to proceed from here: That is I cannot map my indices properly to the ones of the law of total expectation.
