Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

suppose I have independent random variables $X_i$ which are distributed binomially via $$X_i \sim \mathrm{Bin}(n_i, p_i)$$.

Are there relatively simple formulae or at least bounds for the distribution $$S = \sum_i X_i$$ avaialble?

share|cite|improve this question
up vote 8 down vote accepted

See this paper (The Distribution of a Sum of Binomial Random Variables by Ken Butler and Michael Stephens).

share|cite|improve this answer
I have the same question and i read the paper (The Distribution of a Sum of Binomial Random Variables by Ken Butler and Michael Stephens). unfortunately the approximations are not clear to me ( for example how are the probabilities in Table 2 calculated?) – May Dec 6 '12 at 22:33
@May: I had the same problem these days and I ended up using the explicit formula given in the linked paper. Should be fine if you don't have too many samples. – Michael Kuhn Dec 7 '12 at 20:26
@MichaelKuhn: so you used poisson binomial distribution function, unfortunately I have many samples and I need to use an approximation. – May Dec 11 '12 at 0:16
@May: These seems to be an R package for this distribution: – Michael Kuhn Dec 11 '12 at 9:09
Page Not Found ... – Dor Sep 13 '15 at 0:30

This answer provides an R implementation of the explicit formula from the paper linked in the accepted answer (The Distribution of a Sum of Binomial Random Variables by Ken Butler and Michael Stephens). (This code can in fact be used to combine any two independent probability distributions):

# explicitly combine two probability distributions, expecting a vector of 
# probabilities (first element = count 0)
combine.distributions <- function(a, b) {

    # because of the following computation, make a matrix with more columns than rows
    if (length(a) < length(b)) {
        t <- a
        a <- b
        b <- t

    # explicitly multiply the probability distributions
    m <- a %*% t(b)

    # initialized the final result, element 1 = count 0
    result <- rep(0, length(a)+length(b)-1)

    # add the probabilities, always adding to the next subsequent slice
    # of the result vector
    for (i in 1:nrow(m)) {
        result[i:(ncol(m)+i-1)] <- result[i:(ncol(m)+i-1)] + m[i,]


a <- dbinom(0:1000, 1000, 0.5)
b <- dbinom(0:2000, 2000, 0.9)

ab <- combine.distributions(a, b)
ab.df <- data.frame( N = 0:(length(ab)-1), p = ab)

plot(ab.df$N, ab.df$p, type="l")
share|cite|improve this answer

One short answer is that a normal approximation still works well as long as the variance $\sigma^2 = \sum n_i p_i(1-p_i)$ is not too small. Compute the average $\mu = \sum n_i p_i$ and the variance, and approximate $S$ by $N(\mu,\sigma)$.

share|cite|improve this answer
Unfortunately, I cannot say anything about the Variance. In what direction would the normal approximation go? – Lagerbaer Mar 30 '11 at 22:31
Do you mean, "is the normal approximation an overestimate or an underestimate?" That depends on the range of values you are considering. Both distributions have total mass $1$. – Douglas Zare Mar 30 '11 at 23:29
I'd be interested in an estimate on the expected value. – Lagerbaer Mar 31 '11 at 3:03
To use a normal approximation, you have to know the mean and variance. (To use the more complicated approximations in the paper PEV cited, you need more information, such as the first 4 moments.) If you don't know the expected value, then what do you know about these binomial summands? – Douglas Zare Apr 3 '11 at 4:41
I know $n_i$ and $p_i$ of each of the summands, and hence I know the expected value and variance of each of the summands. – Lagerbaer Apr 3 '11 at 4:44

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.