The independence assumption implies the cumulative probabilities (chances that $X$ or fewer occur) can be computed by means of convolution. Convolutions are fast and efficient to compute; far faster than enumerating all the possibilities. (Although with just a dozen variables or so there are just $2^{12} \approx 2000$ possibilities, the number grows exponentially: beyond $30$ or so variables it would take too long to do the calculations by enumeration.)
Here is an R
implementation of the convolution.
convolution <- function(y) {
w <- z <- y[, 1]
apply(y[, -1, drop=FALSE], 2, function(x) w <<- convolve(w, rev(x), type="open"))
f <- cumsum(w)
names(f) <- 0:(length(f)-1)
return(f)
}
The input is supposed to be a matrix of probability distributions, one per column. For a binomial outcome as in this question, there will be just two rows: the first row contains the chance that the outcome is zero and the second row, the chance that it is one. For example, suppose three events have chances $4/5$, $2/3$, and $1/2$ of occurring. The input columns would be $(1/5, 4/5)$, $(1/3, 2/3)$, and $(1/2, 1/2)$ in any order. Thus:
y <- c(4/5, 2/3, 1/2)
convolution(rbind(1-y, y))
will specify the second row of the matrix, compute the first row (by subtracting the second row from unit), and perform the convolution. The output is
0 1 2 3
0.03333333 0.26666667 0.73333333 1.00000000
The values of $X$ occur on the top row (the names) and the cumulative probabilities on the second row. For instance,
0.03333333 is the chance that $X \le 0$, whence $0.96777777 = 1 - 0.03333333$ is the chance that at least $1$ event occurs.
0.26666667 is the chance that $X \le 1$, whence $0.7333333 = 1 - 0.26666667$ is the chance that at least $2$ events occur.
0.73333333 is the chance that $X \le 2$, whence $0.2666667 = 1 - 0.73333333$ is the chance that at least $3$ events occur.
1.00000000 is the chance that $X \le 3$, whence $0 = 1 - 1$ is the chance that at least $4$ events occur.
(You can ask R
to compute these chances directly via 1 - convolution(rbind(1-y, y))
.)
This approach comes to the fore with larger datasets. For instance, let's consider $1000$ independent binomial distributions. (I need to use that many because the calculation is too fast to be timed with smaller numbers!) I'll generate their parameters at random and time the computation of the convolution:
p <- rbeta(10^3,2,1)
system.time(z <- 1 - convolution(rbind(1-p, p)))
user system elapsed
0.34 0.02 0.36
(Those are seconds of computing time.)