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