Non-probabilistic proof for this infinite series identity related to the generalized negative binomial distribution? Using the alternative "Poisson-like" parameterization of the generalized negative binomial distribution in terms of its mean $\lambda$, it follows that the PMF can be written
$$ \frac{\Gamma(r+n)}{\Gamma(r)} \frac{r^r}{(r+\lambda)^{r+n}} \frac{\lambda^n}{n!} \,,$$
and therefore its mean $\lambda$ must equal
$$\sum_{n=0}^{\infty} n  \frac{\Gamma(r+n)}{\Gamma(r)} \frac{r^r}{(r+\lambda)^{r+n}} \frac{\lambda^n}{n!} = \lambda \,.$$
Pulling out all terms in the above sum which don't depend on $n$ leads to the identity:

$$\Gamma(r) \left(\frac{r+\lambda}{r}\right)^r \lambda = \sum_{n=1}^{\infty}  \frac{\Gamma(r+n)}{(n-1)!}  \left(\frac{\lambda}{r + \lambda}  \right)^n \,.  $$
Question: How could one prove the above identity "directly"/"non-probabilistically", presumably using properties of the $\Gamma$ function?

My guess is that because $\frac{\Gamma(r+n)}{\Gamma(r)} \approx r^n$ for $n$ large due to Stirling's approximation. Would someone just plug that into the first expression for the mean, and then try to apply the squeeze theorem? Can that really be good enough when the approximation is only good for large $n$ -- how could one possibly show that the error from smaller $n$ "disappears" as well? Specifically I tried doing this and got several unexpected additional factors.
 A: Let $a:=\lambda/(r+\lambda)$. By monotone convergence, swap integral and sum in the definition of the Gamma Function:
\begin{align*}
\sum_{n=1}^\infty \frac{\Gamma(r+n)}{(n-1)!}a^n &=
\int_0^{\infty}e^{-u}\sum_{n=1}^\infty \frac{1}{(n-1)!} u^{r+n-1} a^ndu\\
&=\int_0^\infty e^{-u}u^rae^{ua}du\\
&=a(1-a)^{-r-1}\Gamma(r+1)\\
&=\frac{\lambda}{\lambda+r}\left(\frac{r+\lambda}{r}\right)^{r+1}\Gamma(r+1)\\
&=\lambda \left(\frac{r+\lambda}{r}\right)^r\Gamma(r).
\end{align*}
A: Side note: I think you want the LHS to have $(\frac{r+\lambda}{r})^r$ instead of $(\frac{r}{r+\lambda})^r$. In any case, in the definition of the distribution, $r$ is a positive integer. Therefore, you get that $\Gamma(r) = (r-1)!$, and the identity can then be proven non-probabilistically by using the identity
$$\frac{d^k}{dx^k} \frac{1}{1-x} = \frac{k!}{(1-x)^{k+1}}$$
as well as noting that this also equals
$$\frac{d^k}{dx^k} \sum_{m=0}^\infty x^m = \sum_{m=k}^\infty m(m-1)...(m-k+1)x^{m-k}$$
$$=\sum_{m=0}^\infty (m+k)(m+k-1)...(m+1)x^{m}.$$
Of course you would be plugging in $x = (\frac{\lambda}{r+\lambda})$. Hope this helps!
