16

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?

suomynonA
  • 6,895
Lagerbaer
  • 3,496

5 Answers5

11

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

NebulousReveal
  • 13,935
  • 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 06 '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 07 '12 at 20:26
  • 1
    @MichaelKuhn: so you used poisson binomial distribution function http://en.wikipedia.org/wiki/Poisson_binomial_distribution, unfortunately I have many samples and I need to use an approximation. – May Dec 11 '12 at 00:16
  • @May: These seems to be an R package for this distribution: http://cran.r-project.org/web/packages/poibin/poibin.pdf – Michael Kuhn Dec 11 '12 at 09:09
  • @MichaelKuhn: thank you. I also found this paper http://www.sciencedirect.com/science/article/pii/S0167947312003568# – May Dec 11 '12 at 15:40
  • @MichaelKuhn: Unfortunately I use MATLAB in my work, and I could not find MATLAB package for this distribution. – May Dec 11 '12 at 21:00
  • 1
    Page Not Found ... – Dor Sep 13 '15 at 00:30
  • It is not beautiful, but I implemented the algorithm from this paper in R for another project a few years ago at https://github.com/jvbraun/sumbin – jvbraun Sep 18 '16 at 18:33
6

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,]
    }

    result
}

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")
Michael Kuhn
  • 161
  • 1
  • 2
3

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

Douglas Zare
  • 2,963
  • Unfortunately, I cannot say anything about the Variance. In what direction would the normal approximation go? – Lagerbaer Mar 30 '11 at 22:31
  • 1
    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 03: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 03 '11 at 04: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 03 '11 at 04:44
2

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.

Thomas Ahle
  • 4,612
0

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

https://cran.r-project.org/web/packages/sinib/index.html

Paper:

https://journal.r-project.org/archive/2018/RJ-2018-011/RJ-2018-011.pdf

boxiangliu
  • 1