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$$ available?


5 Answers 5


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

  • $\begingroup$ 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?) $\endgroup$
    – May
    Dec 6, 2012 at 22:33
  • $\begingroup$ @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. $\endgroup$ Dec 7, 2012 at 20:26
  • 1
    $\begingroup$ @MichaelKuhn: so you used poisson binomial distribution function en.wikipedia.org/wiki/Poisson_binomial_distribution, unfortunately I have many samples and I need to use an approximation. $\endgroup$
    – May
    Dec 11, 2012 at 0:16
  • $\begingroup$ @May: These seems to be an R package for this distribution: cran.r-project.org/web/packages/poibin/poibin.pdf $\endgroup$ Dec 11, 2012 at 9:09
  • 1
    $\begingroup$ Page Not Found ... $\endgroup$
    – Dor
    Sep 13, 2015 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")

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)$.

  • $\begingroup$ Unfortunately, I cannot say anything about the Variance. In what direction would the normal approximation go? $\endgroup$
    – Lagerbaer
    Mar 30, 2011 at 22:31
  • 1
    $\begingroup$ 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$. $\endgroup$ Mar 30, 2011 at 23:29
  • $\begingroup$ I'd be interested in an estimate on the expected value. $\endgroup$
    – Lagerbaer
    Mar 31, 2011 at 3:03
  • $\begingroup$ 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? $\endgroup$ Apr 3, 2011 at 4:41
  • $\begingroup$ 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. $\endgroup$
    – Lagerbaer
    Apr 3, 2011 at 4:44

It is possible to get a Chernoff bound using the standard moment generating function method: $$ \begin{align} \Pr[S\ge s] &\le E[\exp[t \sum_i X_i]]\exp(-st) \\&= \exp\left(\sum_i 1 + (e^t-1) p_i\right) \exp(-st) \\&\le \exp\left(\sum_i \exp((e^t-1) p_i)-st\right) \\&= \exp\left(s-\sum_ip_i-s\log\frac{s}{\sum_i p_i}\right) \end{align}, $$ where we took $t=\log(s/\sum_ip_i)$. This is basically equal to the standard Chernoff bound for equal probabilities, just replaced with the sum (or average if you set $s=n s'$.)

Here we (surprisingly) used the inequality $1+x\le e^x$, but a slightly stronger bound may be possible without it. It'll just be much more messy.

Another way to look at the bound is that we bound each variable with a poisson distribution with the same mean.


If you are looking for a software package to approximate the sum of non-identical binomial random variables, you can try this one:





Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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