# Finite double summation - the upper limit of the inner summation depends on the outer summation's parameter

I am doing some work on finding the distribution of a sum of two independent random variables. In my actual work these two variables are independent, yet have a very different distribution.

I wanted to treat an example case, but I lack the knowledge of double summations. Before starting the example, let me say that I fully know how to obtain the cumulative distribution function in a number of ways, so this question is purely on how to compute double summations of the kind I find in the following example.

Let $X_1 \sim geom(p_1)$, $X_2 \sim geom(p_2)$, both with support $\{1,2,\ldots\}$ and $Y = X_1 + X_2$.

I want to compute the cumulative distribution function of $Y$, using the following approach,

\begin{align} \mathbb{P}(Y \leq k) &= \sum_{j = 1}^{k-1} \sum_{i=1}^{k-j} \mathbb{P}(X_1=i)\mathbb{P}(X_2=j) \\ &= \sum_{j = 1}^{k-1} \sum_{i=1}^{k-j} p_1(1-p_1)^{i-1} p_2(1-p_2)^{j-1}. \end{align}

Here is where I get stuck. How to compute the above double summation? I know, using Mathematica to evaluate the double summation, that the answer should be

\begin{align} \mathbb{P}(Y \leq k) &= \frac{p_1-p_1(1-p_2)^k+p_2((1-p_1)^k-1)}{p_1-p_2}, \end{align}

which is indeed correct when compared to e.g. a generalized negative binomial distribution with correct parameters.

I should state that the following steps depend on the formula for a geometric series:

$$\sum_{j=1}^{k} r^{j-1} = \frac{1-r^k}{1-r}$$

The sum in question is

$$p_1 p_2 \sum_{j=1}^{k-1} (1-p_2)^{j-1} \sum_{i=1}^{k-j} (1-p_1)^{i-1}$$

Evaluate the inner sum first:

$$\sum_{i=1}^{k-j} (1-p_1)^{i-1} = \frac{ 1 - (1-p_1)^{k-j} }{p_1}$$

Now the sum is

\begin{align} p_2 \sum_{j=1}^{k-1} (1-p_2)^{j-1} \left [ 1 - (1-p_1)^{k-j} \right ] &= 1 - (1-p_2)^{k-1} - p_2 \sum_{j=1}^{k-1} (1-p_2)^{j-1} (1-p_1)^{k-j} \\ \end{align}

The sum on the right-hand side may be evaluated as follows:

\begin{align}\sum_{j=1}^{k-1} (1-p_2)^{j-1} (1-p_1)^{k-j} &= (1-p_1)^{k-1} \sum_{j=1}^{k-1} \left ( \frac{1-p_2}{1-p_1} \right )^{j-1} \\ &= (1-p_1)^{k-1} \frac{ 1 - \left ( \frac{1-p_2}{1-p_1} \right )^{k-1}}{1 - \frac{1-p_2}{1-p_1}} \\ &= (1-p_1)\frac{(1-p_1)^{k-1} - (1-p_2)^{k-1}}{p_2-p_1} \\ \end{align}

Now we can put this all together:

\begin{align} 1 - (1-p_2)^{k-1} - p_2 (1-p_1)\frac{(1-p_1)^{k-1} - (1-p_2)^{k-1}}{p_2-p_1} \\ = 1 - \frac{(p_2-p_1)(1-p_2)^{k-1} + p_2 (1-p_1)[(1-p_1)^{k-1} - (1-p_2)^{k-1}]}{p_2-p_1} \\ \end{align}

which, after some cancellation and consolidation, produces the following result for the sum:

$$1 - \frac{p_1 (1-p_2)^k - p_2 (1-p_1)^k}{p_1-p_2}$$

and is equal to the result stated above.

• Thank you for you elaborate response, much appreciated. – user60307 Jan 31 '13 at 13:08
• You're welcome! To be honest, I like the form of the result I posted better than what you posted from Mathematica. You can see the symmetry in the behavior of the sum with respect to $p_1$ and $p_2$. Anyway, that was a lot of fun; feel free to post any more like this. – Ron Gordon Jan 31 '13 at 13:25
• I agree, it clearly shows that when $k$ increases, the term on the right goes to zero and the asymptote is 1. Thank you for your comments. – user60307 Jan 31 '13 at 13:30

If all you need is a decent approximation using a stratified sampling Monte Carlo approach:

1. Generate a table of values for the percentiles 0.5%, 1.5%,...,99.5% for each distribution (a finer mesh could be used for more accurate results)
2. Calculate all possible sums (or any other function of the two random variables)
3. Sort the values smallest to largest, count the n results (i.e. 10,000 for example in line 1), then assign percentiles 1/2n, 3/2n,..., (2n-1)/2n