How to show that the negative binomial CDF converges to the Poisson CDF? (Incomplete beta vs incomplete gamma functions) Question: Is there a straightforward proof of the following relationship between the (lower, non-regularized) incomplete beta function $\mathcal{B}(x; a ,b)$ and the (upper, non-regularized) incomplete gamma function $\Gamma(s,t)$:
$$\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right) \sim \frac{\Gamma(\alpha) \Gamma(k+1, \mu)}{\Gamma(\alpha + k + 1)}  \quad \text{as } \alpha \to \infty ,$$
where $f(\alpha) \sim g(\alpha)$ as $\alpha \to \infty$ means $\displaystyle\lim_{\alpha \to \infty} \frac{f(\alpha)}{g(\alpha)}=1$?
Notice in particular that this would be analogous to the relationship $\mathcal{B}(\alpha, k+1) = \frac{\Gamma(\alpha) \Gamma(k+1)}{\Gamma(\alpha + k + 1)}$ that one has for the (complete) beta function and the (complete) gamma function. I believe I was able to show that the above asymptotic relationship would be equivalent to showing the convergence of the negative binomial CDF to the Poisson CDF.
Background+definitions: The negative binomial distribution with "number of successes" $\alpha > 0$, $\alpha \in \mathbb{R}$, expectation/mean $\mu > 0$, $\mu \in \mathbb{R}$, and "probability of success" $\frac{\alpha}{\alpha + \mu}$ has the following CDF for $k \ge 0$, $ k \in \mathbb{N}$:
$$\frac{\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)}{\mathcal{B}(\alpha, k+1)}, \quad \text{where} \quad \mathcal{B}(x; a, b) := \int_{0}^x t^{a-1} (1-t)^{b-1} \mathrm{d}t . $$
The Poisson distribution with expectation/mean $\mu > 0$, $\mu \in \mathbb{R}$ has the following CDF for $k \ge 0$, $k \in \mathbb{N}$:
$$\frac{\Gamma(k+1, \mu)}{\Gamma(k+1)}, \quad \text{where} \quad \Gamma(s,t) := \int_{t}^{\infty} e^{-u} u^{s-1} \mathrm{d}u . $$
Note: One can also show the convergence in distribution using either the PMF or MGF, and that it is easier than using the CDF. However, I want to understand the proof specifically using the CDFs. (As a sanity check that the formulas given for the CDFs are at least plausible.)
Here are some related, but distinct, questions on this website: (1) (2) (3) (4) (5) (6) (7)
As a community wiki answer I typed up an outline of what I tried so far, which either does not work or which I was not able to complete. Either way the calculations involved are not straightforward.

The general approach of using a change of variables does work, with a slight modification of the strategy:

*

*Engineer the change of variables so that the lower bound of the integral approaches $\mu$ as $\alpha \to \infty$.

Specifically, consider the change of variables $u = \alpha (- \ln t)$, e.g. $\mathrm{d}t = -\frac{1}{\alpha} e^{-u/\alpha} \mathrm{d}u$, so that
$$\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right) = 
 \int^{\infty}_{\alpha \ln \left( \frac{\alpha + \mu}{\alpha} \right)}
e^{-u}  \frac{(1 - e^{-\frac{u}{\alpha}})^k }{\alpha}  \mathrm{d}u , \tag{$\ddagger$}$$
where the right-hand side has been simplified, in particular using $\int_b^a -f(x) \mathrm{d}x = \int_a^b f(x) \mathrm{d}x$.
Notice that $e^{-u}$ is already part of the integrand for the incomplete gamma function, and that
$$\lim_{\alpha \to \infty} \alpha \ln \left( \frac{\alpha+\mu}{\alpha} \right) = \lim_{\alpha \to \infty} \ln\left(\left(\frac{\alpha+\mu}{\alpha} \right)^{\alpha} \right) = \ln \left(\lim_{\alpha \to \infty} \left( \frac{\alpha+\mu}{\alpha}\right)^{\alpha}   \right) = \ln(e^{\mu})  = \mu$$
because $\lim_{\alpha \to \infty} \left(\frac{\alpha+\mu}{\alpha}\right)^{\alpha} = e^{\mu}$, so the lower bound of $(\ddagger)$ really does approach $\mu$ as $\alpha \to \infty$.
Pushing the Pochhammer symbol $\frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)}$ under the integral sign, i.e. writing
$$\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right) \frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)} = 
 \int^{\infty}_{\alpha \ln \left( \frac{\alpha + \mu}{\alpha} \right)}
e^{-u}  \frac{(1 - e^{-\frac{u}{\alpha}})^k }{\alpha} \frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)}  \mathrm{d}u , \tag{$\ddagger\ddagger$} $$
the above observations mean that we are done as soon as we show that
$$\frac{(1 - e^{-\frac{u}{\alpha}})^k }{\alpha} \frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)} \sim u^k \quad \text{as }\alpha \to \infty . \tag{$\ddagger\ddagger\ddagger$} $$
The desired asymptotic relationship $(\ddagger\ddagger\ddagger)$ follows from these two asymptotic relationships:

*

*$\frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)} \sim \alpha^{k+1}$ as $\alpha \to \infty$ which was mentioned in the original question,

*$e^{-\frac{u}{\alpha}} \sim \frac{\alpha}{\alpha + u} = \left(\left( \frac{\alpha}{\alpha + u} \right)^{\alpha} \right)^{\frac{1}{\alpha}} $ as $\alpha \to \infty$, related to $\lim_{\alpha \to \infty} \left( \frac{\alpha}{\alpha + u} \right)^{\alpha} = e^{-u}$.

Substituting them into the left-hand side of $(\ddagger \ddagger \ddagger)$ we get
$$\frac{(1 - e^{-\frac{u}{\alpha}})^k }{\alpha} \frac{\Gamma(\alpha+k+1)}{\Gamma(\alpha)} \sim \frac{(1 - \frac{\alpha}{\alpha + u})^k}{\alpha} \alpha^{k+1}  =  \left(\frac{\alpha}{\alpha + u}\right)^k u^k \sim u^k \quad \text{as }\alpha \to \infty . $$
Note: The idea for a change of variables that turns the incomplete beta integral into something more closely resembling an upper incomplete gamma integral, $v = - \ln t$, comes from equation (3.3) of Uniform Asymptotic Expansions of the Incomplete Gamma Functions and the Incomplete Beta Function by N.M. Temme, Mathematics of Computation, vol. 29, no. 132,
October 1975, pp. 1109-1114.
