Marginal distribution is negative binomial under Poisson distribution with Gamma prior

Question

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$$

Answer 1

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}$.