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.
 A: 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.
A: 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.
A: (community wiki) Failed attempt: What we want to show is
$$ \lim_{\alpha \to \infty} \frac{\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)}{\mathcal{B}(\alpha, k+1)} = \frac{\Gamma(k+1, \mu)}{\Gamma(k+1)} . \tag{$\dagger$}$$
Using $\mathcal{B}(\alpha, k+1) = \frac{\Gamma(\alpha) \Gamma(k+1)}{\Gamma(\alpha + k + 1)}$, this is the same as showing
$$\lim_{\alpha \to \infty} \frac{\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)\Gamma(\alpha + k + 1)}{\Gamma(\alpha) \Gamma(k+1)} = \frac{\Gamma(k+1, \mu)}{\Gamma(k+1)}  .$$
Obviously we can can cancel out the $\Gamma(k+1)$'s in the denominators to get
$$\lim_{\alpha \to \infty} \frac{\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)\Gamma(\alpha + k + 1)}{\Gamma(\alpha)} = \Gamma(k+1, \mu) , \tag{**}$$
which means that proving $(\dagger)$ is equivalent to showing
$$\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 , \tag{*}$$
as I had claimed at the beginning of the question.
Based on how analogous $(*)$ is to the identity related the (complete) beta and gamma functions, I have to imagine that the slickest possible proof of $(\dagger)$ would prove $(*)$ directly, using some facts about special functions that I don't know about.
However the closest I've come so far to attempting to prove $(\dagger)$ is by attacking $(**)$ via brute force and not using any facts about special functions (because I don't know many).
The first observation is that the ratio of the $\Gamma$ functions in $(**)$ is the Pochhammer symbol a.k.a. rising factorial :
$$\frac{\Gamma(\alpha + k + 1)}{\Gamma(\alpha)}  = \prod_{i=0}^k (\alpha + i) .$$
This follows using induction from the functional equation that shows that the gamma function interpolates the factorials, $\Gamma(x+1)=x \Gamma(x)$. This is also related to how $\frac{\Gamma(\alpha + k + 1)}{\Gamma(\alpha)} \sim \alpha^{k+1}$ as $\alpha \to \infty$. Therefore we can rewrite $(**)$ as
$$ \lim_{\alpha \to \infty} \left[\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)  \left( \prod_{i=0}^k (\alpha + i) \right) \right] = \Gamma(k+1, \mu) \tag{***} .$$
Then I did the following brute-force change of variables on the integral defining $\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)$ to force its bounds of integration to match those of the integral defining $\Gamma(k+1, \mu)$:
$$ u = \frac{\alpha}{\left( \frac{\alpha}{\alpha + \mu} \right) -t}  - \alpha , \quad \mathrm{d}t = \frac{\alpha}{(\alpha + u)^2} \mathrm{d}u , \quad t = \frac{\alpha(u - \mu)}{(\alpha + \mu)(\alpha + u)}, \quad 1 - t = \frac{\alpha^2 + 2\alpha\mu + \mu u}{(\alpha+\mu)(\alpha + u)} .$$
So in other words I think I showed that
$$\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right) = \int_{\mu}^{\infty} \left( \frac{\alpha(u - \mu)}{(\alpha + \mu)(\alpha + u)} \right)^{\alpha - 1} \left( \frac{\alpha^2 + 2\alpha\mu + \mu u}{(\alpha+\mu)(\alpha + u)} \right)^k \frac{\alpha}{(\alpha + u)^2} \mathrm{d}u . $$
Then, using that $\frac{\Gamma(\alpha + k + 1)}{\Gamma(\alpha)} = \prod_{i=0}^k (\alpha + i) \sim \alpha^{k+1}$ as $\alpha \to \infty$, pulling $\frac{\Gamma(\alpha + k + 1)}{\Gamma(\alpha)} \sim \alpha^{k+1}$ inside of the integral, and then simplifying and grouping terms, I think I was able to show that
$$\lim_{\alpha \to \infty} \left[\mathcal{B}\left(\frac{\alpha}{\alpha + \mu} ; \alpha, k + 1 \right)  \left( \prod_{i=0}^k (\alpha + i) \right) \right] = \lim_{\alpha \to \infty} \int_{\mu}^{\infty} \left( \frac{\alpha}{\alpha + u}  \right)^{\alpha + k + 1} \cdot \frac{(\alpha^2 + 2 \alpha \mu + \mu u)^k (u-\mu)^{\alpha - 1}}{(\alpha + \mu)^{\alpha + k -1}} \mathrm{d}u . $$
Now it's straightforward enough to show that
$$ \lim_{\alpha \to \infty} \left( \frac{\alpha}{\alpha + u}  \right)^{\alpha + k + 1} = \lim_{\alpha \to \infty} \left( \frac{\alpha + u}{\alpha}  \right)^{-\alpha} = \left(\lim_{\alpha \to \infty}\left( 1 + \frac{u}{\alpha} \right)^{\alpha} \right)^{-1} = e^{-u} ,$$
so we would actually be done now if it is possible to at least show that
$$ \frac{(\alpha^2 + 2 \alpha \mu + \mu u)^k (u-\mu)^{\alpha - 1}}{(\alpha + \mu)^{\alpha + k -1}} \sim u^k \quad \text{as }\alpha \to \infty \tag{****} .$$
However, that seems to be more wishful thinking on my part than anything else, and I probably made an algebra or other error on the way there, so the result $(****)$ probably is not true. And even if it is, the calculations involved to get there are so tedious and not straightforward and not intuitive that it would not be a very straightforward proof.
On the other hand, if there were some fact about special functions that would allow one to show $(*)$ in one fell swoop, that in my opinion would be a much "slicker" argument. So that's why I phrased the question as asking about $(*)$.
