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.
Thanks in advance.

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.
 A: 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.
A: If all you need is a decent approximation using a stratified sampling Monte Carlo approach:


*

*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)

*Calculate all possible sums (or any other function of the two random variables)

*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

