# How to numerically evaluate the CDF of this random variable?

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.

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.
• Thanks for the answer. I understand the normal approximation for the binomial distribution, but am stuck on your first sentence. Writing out the sum $Y+Z$ as a convolution, I obtain the following: $$\sum_{t=0}^x\binom{n}{t}(p+q-pq)^t(1-p-q+pq)^{n-t}\binom{m-n}{x-t}q^{x-t}(1-q)^{m-n-x+t}.$$ Since $1-p-q+pq=(1-p)(1-q)$, I can simplify the convolution as follows: $$\sum_{t=0}^x\binom{n}{t}(p+q-pq)^t(1-p)^{n-t}\binom{m-n}{x-t}q^{x-t}(1-q)^{m-x}$$ but I can't get the expression in (1). Could please elaborate how you get there? – M.B.M. Mar 26 '14 at 23:00
• See Edit. The key is to compute every $P(X=x)$ in one stroke. – Did Mar 27 '14 at 7:15
• Thanks for the edit! I still have questions. First, what is $A_x$? It is summation over $x$ from $x=0$ to $x=m$? I summed it that way (using Mathematica) and it spat out an expression involving ${}_2F_1$... – M.B.M. Mar 27 '14 at 7:37
• From $x=t$ to $x=m$. – Did Mar 27 '14 at 7:42
• Ahh, I see -- there is $x-t$ in the binomial. I must say, this is a fiendishly clever answer. I never though of using generating functions for this (and, as the matter of fact, until now I didn't know what generating functions were good for other than exercises in a graduate probability class...) Thanks! – M.B.M. Mar 27 '14 at 7:53