Distribution of the Negative Binomial sum of Bernoulli random variables I'm trying to find the probability distribution of the random sum $Y$ of $N$ random $X_i$ variables where:
$$ Y = I(N > 0) \sum^N_{i=1} X_i $$
$$ N \sim \mathrm{NegBin}(r, \theta)$$
$$ X_i | p_i \sim \mathrm{Bernoulli}(p_i)$$
$$ p_i \sim \mathrm{Beta}(\alpha, \beta)$$
The setting is similar to a question I had asked recently (available here) but this time, instead of having a Poisson sum, I now have a Negative Binomial sum. 
Using results from the other question (e.g. that each $X_i \sim \mathrm{Bernoulli}(\frac{\alpha}{\alpha + \beta})$ but also simulations (comparing theoretical with simulated means and variances as well as PDFs), I believe that the answer must be:
$$ Y \sim \mathrm{NegBin} \bigg( r, \frac{\alpha \theta}{\alpha + \beta (1 - \theta)} \bigg)$$
However, I'm trying to derive this result but I'm stuck. Here's where (I think) I managed to get to, so far:
$$\mathbb{P}(Y = k) = \sum^\infty_{n=k} \mathbb{P}(Y=k \mid N=n) \mathbb{P}(N=n) = $$
$$ = \sum^\infty_{n=k}{{n}\choose{k}} 
\bigg( \frac{\alpha}{\alpha + \beta} \bigg)^k 
\bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg)^{n-k}
{{n + r - 1}\choose{n}}
(1 - \theta)^r \theta^n = $$
$$ = \frac{(1 - \theta)^r}{k! (r-1)!} \bigg( \frac{\alpha}{\beta} \bigg)^k
\sum^\infty_{n=k}\frac{(n + r - 1)!}{(n - k)!}
\bigg[ \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{n}$$
I'm not sure how to proceed from here. My gut feeling is that I need to make some sort of substitution but most of the stuff I've been taught during my University days has long deserted me! Maybe my initial belief about the result is worng, I may have made a mistake somewhere along the way, or perhaps I'm on the right track but lack the knowledge to convert that into the Negative Binomial PDF. Eitherway, I'd appreciate any suggestions on how to solve this.
UPDATE:
I'm not sure if it helps but we can substitute $j = n-k$ in the last summation. I managed to get:
$$\mathbb{P}(Y = k) = \frac{(1 - \theta)^r}{k! (r-1)!} \bigg[ \frac{\alpha \theta}{\alpha + \beta} \bigg]^k 
\sum^\infty_{j=0}\frac{(j + k + r - 1)!}{j!}
\bigg[ \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{j}$$
I still can't see how to solve the summation. Any help will be much appreciated...
 A: With some help from an offline friend on the summation, I've managed to prove it. For the sum:
$$\sum^\infty_{j=0}\frac{(j + k + r - 1)!}{j!}
\bigg[ \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{j}$$
let $x = \theta ( 1 - \frac{\alpha}{\alpha + \beta} )$ and $m = k + r - 1$. To solve it, and since $|x| < 1$, consider the power series:
$$\sum^\infty_{j = -m} x^{j + m} = \frac{1}{1 - x} = (1 - x)^{-1}$$
Differentiating both sides with respect to $x$ gives:
$$\sum^\infty_{j = -m} (j + m) \space  x^{j + m - 1} = (1 - x)^{-2}$$
The first summation term for $j = -m$ is zero so:
$$\sum^\infty_{j = -m} (j + m) \space  x^{j + m - 1} = 
\sum^\infty_{j = - m + 1} (j + m) \space  x^{j + m - 1} = (1 - x)^{-2}$$
Differentiating again and ignoring the first element which is zero leads to:
$$\sum^\infty_{j = -m + 1} (j + m)(j + m - 1) \space  x^{j + m - 2} = 
\sum^\infty_{j = - m + 2} (j + m) (j + m - 1) \space x^{j + m - 2} = 2(1 - x)^{-3}$$
Repeating the process will eventually lead to:
$$\sum^\infty_{j = 0} (j + m) (j + m - 1) ... (j + 2)(j + 1) \space  x^{j} = m! \space (1 - x)^{-(m+1)}$$
Substituting $x$ and $m$ in the result:
$$\sum^\infty_{j=0}\frac{(j + k + r - 1)!}{j!}
\bigg[ \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{j} = 
(k + r - 1)! \bigg[ 1 - \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{-(k + r)}$$
Using this result in the final (updated) expression of the original question, we get:
$$\mathbb{P}(Y = k) = \frac{(1 - \theta)^r}{k! (r-1)!} \bigg[ \frac{\alpha \theta}{\alpha + \beta} \bigg]^k (k + r - 1)! \bigg[ 1 - \theta \bigg( 1 - \frac{\alpha}{\alpha + \beta} \bigg) \bigg]^{-(k + r)}$$
Grouping together the factorial terms and the terms raised to the power $k$ and $r$ we get:
$$\mathbb{P}(Y = k) = {{k + r - 1}\choose{k}} 
\bigg[\frac{\alpha \theta}{\alpha + \beta (1 - \theta)}\bigg]^k 
\bigg[1 - \frac{\alpha \theta}{\alpha + \beta (1 - \theta)} \bigg]^r,  \space \space k = 0, 1, 2, ...$$
therefore the $\mathrm{Negative Binomial}$ random sum of $\mathrm{Bernoulli}$ random variables each with a success probability following a $\mathrm{Beta}$ distribution follows
$$ Y \sim \mathrm{NegBin} \bigg( r, \frac{\alpha \theta}{\alpha + \beta (1 - \theta)} \bigg)$$
as proposed. 
