# What is the probability distribution of this random sum of non-iid Bernoulli variables?

I'm trying to find the probability distribution of a sum of a random number of variables which aren't identically distributed. Here's an example:

John works at a customer service call center. He receives calls with problems and tries to solve them. The ones he can't solve, he forwards them to his superior. Let's assume that the number of calls he gets in a day follows a Poisson distribution with mean $\mu$. The difficulty of each problem varies from pretty simple stuff (which he can definitely deal with) to very specialized questions which he won't know how to solve. Assume that the probability $p_i$ he'll be able to solve the i-th problem follows a Beta distribution with parameters $\alpha$ and $\beta$ and is independent of the previous problems. What is the distribution of the number of calls he solves in a day?

More formally, I have:

$Y = I(N > 0)\sum_{i = 0}^{N} X_i$ for $i = 0, 1, 2, ..., N$

where $N \sim \mathrm{Poisson}(\mu)$ , $(X_i | p_i) \sim \mathrm{Bernoulli}(p_i)$ and $p_i \sim \mathrm{Beta}(\alpha, \beta)$

Note that, for now, I'm happy to assume that the $X_i$'s are independent. I'd also accept that the parameters $\mu, \alpha$ and $\beta$ do not affect each other although in a real-life example of this when $\mu$ is large, the parameters $\alpha$ and $\beta$ are such so that the Beta distribution has more mass on low success rates $p$. But let's ignore that for now.

I can calculate $P(Y = 0)$ but that's about it. I can also simulate values to get an idea of what the distribution of $Y$ looks like (it looks like Poisson but I don't know if that's down to the numbers of $\mu, \alpha$ and $\beta$ I tried or whether it generalizes, and how it might change for different parameter values). Any idea of what this distribution is or how I could go about deriving it?

Please note that I have also posted this question on TalkStats Forum but I thought that it might get more attention here. Apologies for cross-posting and many thanks in advance for your time.

EDIT: As it turns out (see the very helpful answers below - and thanks for those!), it is indeed a $\mathrm{Poisson}(\frac{\mu\alpha}{\alpha + \beta})$ distribution, something which I was guessing based on my intuition and some simulations, but was not able to prove. What I now find surprising though, is that the Poisson distribution only depends on the mean of the $\mathrm{Beta}$ distribution but is not affected by its variance.

As an example, the following two Beta distributions have the same mean but different variance. For clarity, the blue pdf represents a $\mathrm{Beta}(2, 2)$ and the red one $\mathrm{Beta}(0.75, 0.75)$.

However, they would both result in the same $\mathrm{Poisson}(0.5\mu)$ distribution which, to me, seems slightly counter-intuitive. (Not saying that the result is wrong, just surprising!)

• For fixed $N$ there is Poisson-binomial distribution but your problem is more complicated then this.
– Tim
May 11, 2017 at 9:02
• Thanks, I know of the Poisson-binomial distribution but $N$ is random here. May 11, 2017 at 9:05
• You might take a look at the compound Poisson, but you might need to do some work with the 0s to make that useful May 12, 2017 at 6:10

The calls (that is, the $X_i$) arrive according to a Poisson process. The total number of calls $N$ follows a Poisson distribution. Divide the calls into two types, e.g. whether $X_i = 1$ or $X_i = 0$. The goal is to determine the process that generates the $1$s. This is trivial if $X_i = 1$ with a fixed probability $p$: by the superposition principle of Poisson processes, the full process thinned to just the $1$s would also be a Poisson process, with rate $p\mu$. In fact this is the case, we just require an additional step to get there.

Marginalize over $p_i$, so that $$\mathrm{Pr}(X_i|\alpha, \beta) = \int_0^1 p_i^{X_i} (1-p_i)^{1-X_i} \frac{p_i^{\alpha-1} (1-p_i)^{\beta-1}}{\mathcal{B}(\alpha, \beta)} dp_i = \frac{\mathcal{B}(X_i + \alpha, 1 - X_i + \beta)}{\mathcal{B}(\alpha, \beta)}$$

Where $\mathcal{B}(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a + b)}$ is the beta function. Using the fact that $\Gamma(x+1) = x\Gamma(x)$, the above simplifies to;

$$\mathrm{Pr}(X_i = 1|\alpha, \beta) = \frac{\Gamma(1+\alpha)\Gamma(\beta)}{\Gamma(1+\alpha+\beta)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} = \frac{\alpha}{\alpha+\beta}$$ In other words, $X_i \sim \mathrm{Bernoulli}(\frac{\alpha}{\alpha+\beta})$. By the superposition property, $Y$ is Poisson distributed with rate $\frac{\alpha \mu}{\alpha+\beta}$.

A numerical example (with R) ... in the figure, the vertical lines are from simulation and red points are the pmf derived above:

draw <- function(alpha, beta, mu)
{ N <- rpois(1, mu); p = rbeta(N, alpha, beta); sum(rbinom(N, size=1, prob=p)) }

pmf <- function(y, alpha, beta, mu)
dpois(y, alpha*mu/(alpha+beta))

y <- replicate(30000,draw(4,5,10))
tb <- table(y)

# simulated pmf
plot(tb/sum(tb), type="h", xlab="Y", ylab="Probability")
# analytic pmf
points(0:max(y), pmf(0:max(y), 4, 5, 10), col="red")


1. Since $p_i$ is a random variable with a $\operatorname{Beta}(\alpha,\beta)$ you have $\mathbb{E}[p_i]= \dfrac{\alpha}{\alpha+\beta}$ and this is in fact the probability that John actually solves the $i$th problem, independently of all the others.

2. Since the total number of problems in a day has a Poisson distribution with parameter $\mu$ and each will be solved with probability $\dfrac{\alpha}{\alpha+\beta}$, the number John solves each day has a Poisson distribution with parameter $\dfrac{\mu\alpha}{\alpha+\beta}$

3. Your calculation of the probability he does not solve any problems should be $\mathbb{P}(Y=0) = e^{-{\mu\alpha}/({\alpha+\beta})}$