Suppose we have random variable $Y$ has the Poisson distribution with parameter $\theta$, and $\theta$ has a Gamma prior distribution, i.e.
$$\begin{aligned}
& \text{Data: }\hspace{1cm} y\mid\theta \sim \operatorname{Poisson}(\theta)\\
& \text{Prior: } \hspace{1cm} \theta \sim \operatorname{Gamma}(\alpha,\beta)
\end{aligned}$$
We want to show the marginal distribution of a random sample ${\bf y}=y_1,\cdots,y_n$ follows the Negative Binomial distribution.
My attempt:
$$\begin{aligned}
f({\bf y})& =\int f({\bf y},\theta)\,d \theta\\
& =\int f({\bf y}\mid\theta)f(\theta)\,d\theta\\
& =\int \left[\prod^n_{i=1}f(y_i\mid\theta)\right]f(\theta)\,d\theta\\
& \propto \int \left[\prod^n_{i=1}e^{-\theta} \theta^{y_i} \right] \theta^{\alpha-1}e^{-\beta\theta} \, d\theta\\
& =\int\theta^{\sum^n_{i=1}y_i+\alpha-1}e^{-(n+\beta) \theta} \, d\theta\\
& =\frac{\Gamma(\sum^n_{i=1}y_i+\alpha)}{(n+\beta)^{\sum^n_{i=1}y_i+\alpha}}
\end{aligned}$$
But I did not see why it is a Negative Binomial distribution, since if
$$y \sim \operatorname{NB}(r,p)$$
then the pmf is
$$f(Y=y)=\binom{y+r-1} y p^y(1-p)^r$$
Best Answer
Let $Y\sim \mathcal{P}(\lambda) \underset{\lambda}{\wedge} \Gamma(a,b)$ and $f$ denote the pdf of the gamma distribution with parameters $a$ and $b$. \begin{align*} P(Y=y) & = \int_0^{+\infty} e^{-\lambda} \frac{\lambda^y}{y!} f(\lambda) d\lambda \\ & = \frac{b^a}{\Gamma(a)y!} \underbrace{\int_0^{+\infty} \lambda^{y+a-1} e^{-(1+b)\lambda} d\lambda }_{\displaystyle =\frac{\Gamma(y+a)}{(1+b)^{y+a}}} \\ & = \frac{\Gamma(y+a)}{\Gamma(a)y!} \left(\frac{1}{1+b}\right)^y \left(\frac{b}{1+b}\right)^a \\ P(Y=y) & = \binom{y+a-1}{y} p^y (1-p)^a \end{align*} where $p:=\frac{1}{1+b}$.