How can I efficiently model the sum of Bernoulli random variables?

Question

I am modeling a random variable ($Y$) which is the sum of some ~15-40k independent Bernoulli random variables ($X_i$), each with a different success probability ($p_i$). Formally, $Y=\sum X_i$ where $\Pr(X_i=1)=p_i$ and $\Pr(X_i=0)=1-p_i$.

I am interested in quickly answering queries such as $\Pr(Y<=k)$ (where $k$ is given).

Currently, I use random simulations to answer such queries. I randomly draw each $X_i$ according to its $p_i$, then sum all $X_i$ values to get $Y'$. I repeat this process a few thousand times and return the fraction of times $\Pr(Y'\leq k)$.

Obviously, this is not totally accurate (although accuracy greatly increases as the number of simulations increases). Also, it seems I have enough data about the distribution to avoid the use simulations. Can you think of a reasonable way to get the exact probability $\Pr(Y\leq k)$?

p.s.

I use Perl & R.

EDIT

Following the responses I thought some clarifications might be needed. I will shortly describe the setting of my problem. Given is a circular genome with circumference c and a set of n ranges mapped to it. For example, c=3*10^9 and ranges={[100,200],[50,1000],[3*10^9-1,1000],...}. Note all ranges are closed (both ends are inclusive). Also note that we only deal with integers (whole units).

I am looking for regions on the circle that are undercovered by the given n mapped ranges. So to test whether a given a range of length x on the circle is undercovered, I test the hypothesis that the n ranges are mapped randomly. The probability a mapped range of length q>x will fully cover the given range of length x is (q-x)/c. This probability becomes quite small when c is large and/or q is small. What I'm interested is the number of ranges (out of n) which cover x. This is how Y is formed.

I test my null hypothesis vs. one sided alternative (undercoverage). Also note I am testing multiple hypothesis (different x lengths), and sure to correct for this.

Answer 1

If it often resembles a Poisson, have you tried approximating it by a Poisson with parameter $\lambda = \sum p_i$ ?

EDIT: I've found a theoretical result to justify this, as well as a name for the distribution of $Y$: it's called the Poisson binomial distribution. Le Cam's inequality tells you how closely its distribution is approximated by the distribution of a Poisson with parameter $\lambda = \sum p_i$. It tells you the quality of this approx is governed by the sum of the squares of the $p_i$s, to paraphrase Steele (1994). So if all your $p_i$s are reasonably small, as it now appears they are, it should be a pretty good approximation.

EDIT 2: How small is 'reasonably small'? Well, that depends how good you need the approximation to be! The Wikipedia article on Le Cam's theorem gives the precise form of the result I referred to above: the sum of the absolute differences between the probability mass function (pmf) of $Y$ and the pmf of the above Poisson distribution is no more than twice the sum of the squares of the $p_i$s. Another result from Le Cam (1960) may be easier to use: this sum is also no more than 18 times the largest $p_i$. There are quite a few more such results... see Serfling (1978) for one review.

Answer 2

@onestop provides good references. The Wikipedia article on the Poisson binomial distribution gives a recursive formula for computing the exact probability distribution; it requires $O(n^2)$ effort. Unfortunately, it's an alternating sum, so it will be numerically unstable: it's hopeless to do this computation with floating point arithmetic. Fortunately, when the $p_i$ are small, you only need to compute a small number of probabilities, so the effort is really proportional to $O(n \log(\sum_i{p_i}))$. The precision necessary to carry out the calculation with rational arithmetic (i.e., exactly, so that the numerical instability is not a problem) grows slowly enough that the overall timing may still be approximately $O(n^2)$. That's feasible.

As a test, I created an array of probabilities $p_i = 1/(i+1)$ for various values of $n$ up to $n = 2^{16}$, which is the size of this problem. For small values of $n$ (up to $n = 2^{12}$) the timing for the exact calculation of probabilities was in seconds and scaled quadratically, so I ventured a calculation for $n = 2^{16}$ out to three SDs above the mean (probabilities for 0, 1, ..., 22 successes). It took 80 minutes (with Mathematica 8), in line with the predicted time. (The resulting probabilities are fractions whose numerators and denominators have about 75,000 digits apiece!) This shows the calculation can be done.

An alternative is to run a long simulation (a million trials ought to do). It only has to be done once, because the $p_i$ don't change.

Answer 3

The sum of $n$ Bernoulli random variables is the binomial distribution. So in R to obtain $Pr(Y < k)$, just use the pbinom function:

##k=1000, n=2000 bernoulli, p=0.5
pbinom(1000, 2000, 0.5)
0.5089195

However, your scenario isn't that straightforward as $p$ isn't constant. Since your Bernoulli random variables are independent, I would replace $p$ with $\bar p$. Since $n$ is very large, the error associated with your estimate of $\bar p$ would be very small.

---

In your case you might even consider using the Normal approximation to the binomial distribution:

pnorm(1000, 2000*0.5, 2000*0.25)
0.5

Answer 4

Well, based on your description and the discussion in the comments it is clear that $Y$ has mean $\sum_i p_i$ and variance $\sum_i p_{i}(1-p_{i})$. The shape of $Y$'s distribution will ultimately depend on the behavior of $p_i$. For suitably "nice" $p_i$ (in the sense that not too many of them are really close to zero), the distribution of $Y$ will be approximately normal (centered right at $\sum p_i$). But as $\sum_i p_i$ starts heading toward zero the distribution will be shifted to the left and when it crowds up against the $y$-axis it will start looking a lot less normal and a lot more Poisson, as @whuber and @onestop have mentioned.

From your comment "the distribution looks Poisson" I suspect that this latter case is what's happening, but can't really be sure without some sort of visual display or summary statistics about the $p$'s. Note however, as @whuber did, that with sufficiently pathological behavior of the $p$'s you can have all sorts of spooky things happen, like limits that are mixture distributions. I doubt that is the case here, but again, it really depends on what your $p$'s are doing.

As to the original question of "how to efficiently model", I was going to suggest a hierarchical model for you but it isn't really appropriate if the $p$'s are fixed constants. In short, take a look at a histogram of the $p$'s and make a first guess based on what you see. I would recommend the answer by @mpiktas (and by extension @csgillespie) if your $p$'s aren't too crowded to the left, and I would recommend the answer by @onestop if they are crowded left-ly.

By the way, here is the R code I used while playing around with this problem: the code isn't really appropriate if your $p$'s are too small, but it should be easy to plug in different models for $p$ (including spooky-crazy ones) to see what happens to the ultimate distribution of $Y$.

set.seed(1)
M <- 5000
N <- 15000
p <- rbeta(N, shape1 = 1, shape2 = 10)
Y <- replicate(M, sum(rbinom(N, size = 1, prob = p)))

Now take a look at the results.

hist(Y)
mean(Y)
sum(p)
var(Y)
sum(p*(1 - p))

Have fun; I sure did.

Answer 5

With different $p_i$ your best bet I think is normal approximation. Let $B_n=\sum_{i=1}^np_i(1-p_i)$. Then

\begin{align*} B_n^{-1/2}\left(\sum_{i=1}^nX_i-\sum_{i=1}^np_i\right)\to N(0,1), \end{align*} as $n\to\infty$, provided that for each $\varepsilon>0$

\begin{align*} B_n^{-1}\sum_{i=1}^nE\left((X_i-p_i)^2\mathbf{1}\{|X_i-p_i|>\varepsilon B_n^{1/2}\}\right)\to 0, \end{align*} as $n\to\infty$, which for Bernouli variables will hold if $B_n\to\infty$. This is so called Lindeberg condition which is sufficient and necessary for convergence to the standard normal.

Update: The approximation error can be calculated from following inequality:

\begin{align*} \sup_x|F_n(x)-\Phi(x)|\le AL_n, \end{align*} where \begin{align*} L_n=B_n^{-3/2}\sum_{i=1}^nE|X_i-p_i|^3 \end{align*} and $F_n$ is the cdf of the scaled and centered sum of $X_i$.

As whuber pointed out, the convergence can be slow for bad behaved $p_i$. For $p_i=\frac{1}{1+i}$ we have $B_n\approx \ln n$ and $L_n\approx (\ln n)^{-1/2}$. Then taking $n=2^{16}$ we get that maximum deviation from the standard normal cdf is a whopping 0.3.

Answer 6

I came across your question while searching for a solution to this very problem. I wasn't terrifically satisfied with the answers here, but I think there's a pretty simple solution that gives you the exact distribution, and is pretty tractable.

The distribution of the sum of two discrete random variables is the convolution of their densities. So if you have $Z = X + Y$ where you know $P(X)$ and $P(Y)$ then you can compute:

$$P(Z=z) = \sum_{k=-\infty}^{\infty} P(X=k) \; P(Y=z-k)$$

(Of course for Bernoulli random variables you don't need to go quite to infinity.)

You can use this to find the exact distribution of the sum of your RVs. First sum two of the RVs together by convolving their PDFs (e.g. [0.3, 0.7] * [0.6, 0.4] = [0.18, 0.54, 0.28]). Then convolve that new distribution with your next Bernoulli PDF (e.g. [0.18, 0.54, 0.28] * [0.5, 0.5] = [0.09, 0.36, 0.41, 0.14]). Keep repeating this until all RVs have been added. And voila, the resulting vector is the exact PDF of the sum of all your variables.

I've verified with simulation that this produces the correct results. It doesn't rely on any asymptotic assumptions, and has no requirements that the Bernoulli probs are small.

There may also be some way to do this more efficiently than repeated convolution, but I haven't thought about it very deeply. I hope this is helpful to somebody!

Answer 7

(Because this is approach is independent of the other solutions posted, including one that I have posted, I'm offering it as a separate response).

You can compute the exact distribution in seconds (or less) provided the sum of the p's is small.

We have already seen suggestions that the distribution might approximately be Gaussian (under some scenarios) or Poisson (under other scenarios). Either way, we know its mean $\mu$ is the sum of the $p_i$ and its variance $\sigma^2$ is the sum of $p_i(1-p_i)$. Therefore the distribution will be concentrated within a few standard deviations of its mean, say $z$ SDs with $z$ between 4 and 6 or thereabouts. Therefore we need only compute the probability that the sum $X$ equals (an integer) $k$ for $k = \mu - z \sigma$ through $k = \mu + z \sigma$. When most of the $p_i$ are small, $\sigma^2$ is approximately equal to (but slightly less than) $\mu$, so to be conservative we can do the computation for $k$ in the interval $[\mu - z \sqrt{\mu}, \mu + z \sqrt{\mu}]$. For example, when the sum of the $p_i$ equals $9$ and choosing $z = 6$ in order to cover the tails well, we would need the computation to cover $k$ in $[9 - 6 \sqrt{9}, 9 + 6 \sqrt{9}]$ = $[0, 27]$, which is just 28 values.

The distribution is computed recursively. Let $f_i$ be the distribution of the sum of the first $i$ of these Bernoulli variables. For any $j$ from $0$ through $i+1$, the sum of the first $i+1$ variables can equal $j$ in two mutually exclusive ways: the sum of the first $i$ variables equals $j$ and the $i+1^\text{st}$ is $0$ or else the sum of the first $i$ variables equals $j-1$ and the $i+1^\text{st}$ is $1$. Therefore

$$f_{i+1}(j) = f_i(j)(1 - p_{i+1}) + f_i(j-1) p_{i+1}.$$

We only need to carry out this computation for integral $j$ in the interval from $\max(0, \mu - z \sqrt{\mu})$ to $\mu + z \sqrt{\mu}.$

When most of the $p_i$ are tiny (but the $1 - p_i$ are still distinguishable from $1$ with reasonable precision), this approach is not plagued with the huge accumulation of floating point roundoff errors used in the solution I previously posted. Therefore, extended-precision computation is not required. For example, a double-precision calculation for an array of $2^{16}$ probabilities $p_i = 1/(i+1)$ ($\mu = 10.6676$, requiring calculations for probabilities of sums between $0$ and $31$) took 0.1 seconds with Mathematica 8 and 1-2 seconds with Excel 2002 (both obtained the same answers). Repeating it with quadruple precision (in Mathematica) took about 2 seconds but did not change any answer by more than $3 \times 10^{-15}$. Terminating the distribution at $z = 6$ SDs into the upper tail lost only $3.6 \times 10^{-8}$ of the total probability.

Another calculation for an array of 40,000 double precision random values between 0 and 0.001 ($\mu = 19.9093$) took 0.08 seconds with Mathematica.

This algorithm is parallelizable. Just break the set of $p_i$ into disjoint subsets of approximately equal size, one per processor. Compute the distribution for each subset, then convolve the results (using FFT if you like, although this speedup is probably unnecessary) to obtain the full answer. This makes it practical to use even when $\mu$ gets large, when you need to look far out into the tails ($z$ large), and/or $n$ is large.

The timing for an array of $n$ variables with $m$ processors scales as $O(n(\mu + z \sqrt{\mu})/m)$. Mathematica's speed is on the order of a million per second. For example, with $m = 1$ processor, $n = 20000$ variates, a total probability of $\mu = 100$, and going out to $z = 6$ standard deviations into the upper tail, $n(\mu + z \sqrt{\mu})/m = 3.2$ million: figure a couple seconds of computing time. If you compile this you might speed up the performance two orders of magnitude.

Incidentally, in these test cases, graphs of the distribution clearly showed some positive skewness: they aren't normal.

For the record, here is a Mathematica solution:

pb[p_, z_] := Module[
  {\[Mu] = Total[p]},
  Fold[#1 - #2 Differences[Prepend[#1, 0]] &, 
   Prepend[ConstantArray[0, Ceiling[\[Mu] + Sqrt[\[Mu]] z]], 1], p]
  ]

(NB The color coding applied by this site is meaningless for Mathematica code. In particular, the gray stuff is not comments: it's where all the work is done!)

An example of its use is

pb[RandomReal[{0, 0.001}, 40000], 8]

---

Edit

An R solution is ten times slower than Mathematica in this test case--perhaps I have not coded it optimally--but it still executes quickly (about one second):

pb <- function(p, z) {
  mu <- sum(p)
  x <- c(1, rep(0, ceiling(mu + sqrt(mu) * z)))
  f <- function(v) {x <<- x - v * diff(c(0, x));}
  sapply(p, f); x  
}
y <- pb(runif(40000, 0, 0.001), 8)
plot(y)

Answer 8

I think other answers are great, but I didn't see any Bayesian ways of estimating your probability. The answer doesn't have an explicit form, but the probability can be simulated using R.

Here is the attempt:

$$ X_i | p_i \sim Ber(p_i)$$

$$ p_i \sim Beta(\alpha, \beta) $$

Using wikipedia we can get estimates of $\hat{\alpha}$ and $\hat{\beta}$ (see parameter estimation section).

Now you can generate draws for the $i^{th}$ step, generate $p_i$ from $Beta(\hat{\alpha},\hat{\beta})$ and then generate $ X_i$ from $Ber(p_i)$. After you have done this $N$ times you can get $Y = \sum X_i$. This is a single cycle for generation of Y, do this $M$(large) number of times and the histogram for $M$ Ys will be the estimate of density of Y.

$$Prob[Y \leq y] = \frac {\#Y \leq y} {M}$$

This analysis is valid only when $p_i$ are not fixed. This is not the case here. But I will leave it here, in case someone has a similar question.