How to numerically evaluate the CDF of this random variable?

Question

I have a discrete random variable $X=0,1,2,\ldots$ with the following probability mass function:

$$P(X=x)=\sum_{t=0}^x\binom{n}{t}p^t(1-p)^{n-t}\binom{m-t}{x-t}q^{x-t}(1-q)^{m-x}\tag{1}$$

where $0<p,q<1$, $n\gg0$ and $m=an$ with $a>1$. The values I am using are $p=2\times10^{-3}$, $q=10^{-4}$, $n=10^5$ and $a=30$ (so $m=3\times10^6$). I am interested in making $p$ even smaller (as low as $2\times 10^6$ and maybe lower), and increasing $n$ to as much as $10^6$ (with corresponding increase in $m$).

Essentially, $X$ is a "binomial mixture" of binomial random variables.

I would like to calculate its cumulative distribution $F_X(s)=\sum_{x=0}^{s}P(X=x)$ for relatively small values of $s$ (where most of this distribution should be concentrated). I would happy with tight upper and lower bounds or approximations of the CDF of $X$.

Right now, I am having a hard time even evaluating (1) at the values above. Does anyone have ideas? What am I doing wrong? What else can I try?

WHAT I'VE DONE

1. Tried plugging the values above into Wolfram Mathematica, and just evaluating the sum as is. Result: overflow (due to binomials).

2. Tried sum $\sum_{x=0}^s\hat{p}(x)$ over the following Gaussian approximation of the p.m.f.: $$\hat{p}(x)=\int_{0}^{x}\frac{\exp\left[-\frac{(t-np)^2}{2np(1-p)}\right]}{\sqrt{2\pi np(1-p)}}\frac{\exp\left[-\frac{(x-t-(m-t)q)^2}{2(m-t)q(1-q)}\right]}{\sqrt{2\pi (m-t)q(1-q)}}dt$$

   The approximation is due to the local limit theorem for binomial random variables, however, I am not sure if it applies here. I got $\approx2.85\times 10^{-94}$ evaluating the sum for $s=40$ and the parameters above in Wolfram Mathematica, which isn't right.

Answer 1

Let $r=p+q-pq$, then $X=Y+Z$ where $Y$ is binomial $(n,r)$, $Z$ is binomial $(m-n,q)$, and $Y$ and $Z$ are independent.

Thus, the gaussian approximations are $$Y=nr+\sqrt{nr(1-r)}\,U_n, $$ and $$Z=(m-n)q+\sqrt{(m-n)q(1-q)}\,V_{m-n},$$ where $U_n\to U$ and $V_{m-n}\to V$ in distribution, for some independent standard normal random variables $U$ and $V$. Thus, $$X=np(1-q)+mq+\varrho W_{m,n},$$ where $W_{m,n}\to W$ in distribution, for some standard normal random variable $W$, with $$ \varrho^2=nr(1-r)+(m-n)q(1-q). $$

Edit: To prove the first assertion, one can identify the generating function of some random variable $X$ with the distribution in the question. This is $$ E(s^X)=\sum_xP(X=x)s^x=\sum_t\binom{n}{t}(sp)^t(1-p)^{n-t}A_x, $$ where $$ A_x=\sum_x\binom{m-t}{x-t}(sq)^{x-t}(1-q)^{m-x}=(sq+1-q)^{m-t}, $$ hence $$ E(s^X)=(sq+1-q)^{m-n}\sum_t\binom{n}{t}(sp)^t(1-p)^{n-t}(sq+1-q)^{n-t}, $$ and the sum on the RHS is $$ (sp+(1-p)(sq+1-q))^n=(1-r+rs)^n, $$ hence $E(s^X)=E(s^Y)E(s^Z)$, as desired.