# 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...

• Do you think that {$X_i | p_i \sim \mathrm{Bernoulli}(p_i)$ and $p_i \sim \mathrm{Beta}(\alpha, \beta)$} is the same as { $X_i \sim \mathrm{Bernoulli}(\alpha/(\alpha+\beta)$}? – user158565 May 29 '17 at 17:37
• I thought that this is what was shown in the other thread (stats.stackexchange.com/questions/278897/…). Is it wrong? – Constantinos May 29 '17 at 17:44
• I think the marginal variance of $X_i$ is larger when beta distribution is added than no beta interfering. What do you think? – user158565 May 29 '17 at 17:51
• What exactly do you mean by "marginal variance"? The way I see it, $X_i$ can take values of 0 or 1. If the probability of taking the value 1 is $\frac{\alpha}{\alpha + \beta}$ then it must be a Bernoulli variable with that mean. I also ran a small simulation: iSim <- 500000 ; a <- 4 ; b <- 9 ; set.seed(123) ; p <- rbeta(iSim, shape1 = a, shape2 = b) ; X <- sum(runif(iSim) < p)/iSim The empirical and theoretical mean are the same. That's not a proof of course, but it points to that direction, no? – Constantinos May 29 '17 at 18:51
• I verified that variances are equal ($\alpha \beta / (\alpha + \beta)^2$) – user158565 May 29 '17 at 19:12

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.