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.

Perhaps the simplest elementary yet rigorous proof employs the cumulant-generating functions.  If you insist, you can translate this into integrals involving the distribution functions.  Anything else would be so burdensome as to be of little or no interest from a statistical perspective, IMHO.
Recall that the cgf of a random variable $X$ with distribution $F$ is
$$\psi_F(t) = \log E\left[e^{itX}\right].$$
When $X$ has non-negative support (as is the case in this question), this can be expressed via integration by parts directly in terms of $F$ as
$$\psi_F(t) = \log\int e^{itx}\,\mathrm dF(x) = \log \int_0^\infty \frac{e^{itx} - 1}{it}(1 - F(x))\,\mathrm dx.$$
For the Negative Binomial distribution with parameters $p\in(0,1)$ and $r \gt 0,$ whose probability function is $f(k;p,r) = \binom{k+r-1}{k}(1-p)^k p^r$ and survival function is
$$1 - F(x;p,r) = I_p(x+1,r)$$
(the regularized Incomplete Beta function), the cgf is
$$\psi(t;p,r) = r\left(\log p - \log\left(1 - (1-p)e^{it}\right)\right).$$
Fix a number $\lambda\gt 0$ and let $p = r/(r+\lambda).$ Then
$$\psi(t;p(r,\lambda),r) = r\left(\log \left(1 - \frac{\lambda}{r+\lambda}\right) - \log\left(1 - \frac{\lambda}{r+\lambda}e^{it}\right)\right).$$
Its limiting value (if it exists) as $r\to\infty$ is the limit as $z = \lambda/(r+\lambda)$ shrinks to zero among positive values.  So, recognizing $r =\lambda/z - \lambda,$ rewrite the right hand side as
$$r\left(\log \left(1 - z\right) - \log\left(1 -ze^{it}\right)\right) = \left(\frac{\lambda}{z}-\lambda\right)(-z + ze^{it} + O(z^2)) = \lambda\left(e^{it}-1\right) + O(z).$$
Thus, the limit does exist and equals $\lambda\left(e^{it}-1\right),$ the cgf of the Poisson$(\lambda)$ distribution.  The Lévy Continuity Theorem asserts the limit of the Negative Binomial distribution function is the Poisson distribution function at every continuity point of the latter, which implies it converges everywhere (because all these functions share the same set of points of discontinuity; namely, the natural numbers.)

We are now in a better place to appreciate the content of the Lévy Continuity Theorem and to address the question about analyzing the limits of the distribution functions (or, equivalently, the survival functions).
We have seen what kind of relationship must hold among $p$ and $r$ asymptotically.  Continuing to fix $\lambda \gt 0$ and taking $p = r/(r+\lambda),$ let's analyze the survival function directly, the incomplete Beta function mentioned in the title:
$$1 - F(x;p,r) = I_p(x+1,r) = \frac{1}{B(x+1,r)}\int_0^{1-r/(r+\lambda)} t^x(1-t)^{r-1}\,\mathrm dt.$$
As $r$ grows large, the factor $(1-t)^{r-1}$ in the integrand dwindles except for $t$ extremely close to $0.$ This strongly suggests the classic approximation
$$(1-t)^{r-1} = \left[(1-t)^{1/t}\right]^{t(r-1)} = e^{-t(r-1)} + O(1/r).$$
Accordingly, change the variable of integration from $t$ to $z=t(r-1)$ (and do a little algebra along the way with the Beta function and the upper limit of integration) to find
$$1 - F(x;p,r) = \frac{\Gamma(x+r+1)}{\Gamma(x+1)\Gamma(r)(r-1)^{x+1}}\int_0^{\lambda(r-1)/(r+\lambda)} z^x\left(e^{-z}+O\left(\frac{1}{r}\right)\right)\,\mathrm dz.$$
The way forward is now so clear that it shouldn't be necessary to give the details.  Simply notice that for very large $r,$

*

*The basic relation $\Gamma(x+r+1) = \Gamma(r)(r)(r+1)\cdots(r+x)$ implies $$\frac{\Gamma(x+r+1)}{\Gamma(x+1)\Gamma(r)(r-1)^{x+1}} = \frac{1}{\Gamma(x+1)} + O\left(\frac{1}{r}\right).$$  (If you're struggling with this, take logarithms and expand each one to first order in its Maclaurin series.)


*The error in the integrand is bounded, $$\int_0^{\lambda(r-1)/(r+\lambda)} z^x\left(e^{-z}+O\left(\frac{1}{r}\right)\right)\,\mathrm dz = \int_0^{\lambda(r-1)/(r+\lambda)} z^xe^{-z}\,\mathrm dz + O\left(\frac{1}{r}\right).$$


*The region of integration is close to the simpler one bounded above by $\lambda = \lambda(r-1)/(r+\lambda) + O(1/r),$ $$\int_0^{\lambda(r-1)/(r+\lambda)} z^xe^{-z}\,\mathrm dz = \int_0^\lambda z^xe^{-z}\,\mathrm dz + O\left(\frac{1}{r}\right).$$
The Lévy Continuity Theorem has essentially taken care of all three parts of this analysis.
Putting these results together yields
$$\lim_{r\to \infty}1 - F(x;p(r,\lambda),r) = \frac{1}{\Gamma(x+1)}\int_0^\lambda z^x e^{-z}\,\mathrm dz.$$
This incomplete Gamma function (the other one mentioned in the title of this thread) is the Poisson$(\lambda)$ survival function for all integral $x\ge 0,$ QED.
