# Marginal distribution is negative binomial under Poisson distribution with Gamma prior

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|\theta \sim \text{Poisson}(\theta)\\ & \text{Prior: } \hspace{1cm} \theta \sim \text{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}|\theta)f(\theta)d\theta\\ & =\int \left[\prod^n_{i=1}f(y_i|\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 NB(r,p)$$ then the pmf is $$f(Y=y)=\begin{pmatrix} y+r-1 \\ y \end{pmatrix}p^y(1-p)^r$$

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 }_=\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}$$.