# Which test should I use?

I have a set of probabilities that an event will occur (yes/no). E.g. [0.87, 0.56, 0.97], and I need to know "What is that probability that at least X of these events occurs?". I've been looking into Binomial Probability Density and related tests, but they all seem to assume that the probabilities for each event are the same.

Thanks for your help! I hope this isn't too rudimentary for this forum.

• How large is your set? – whuber Feb 6 '13 at 20:02
• Are the events independent or dependent? Also, if this is homework, it should have the homework tag – Peter Flom - Reinstate Monica Feb 6 '13 at 20:02
• This isn't homework, I'm a web developer and don't have a stats background. The set will generally be around 12. They are independent. – robamaton Feb 6 '13 at 20:07
• Closely related: stats.stackexchange.com/questions/41251/…. – whuber Feb 6 '13 at 21:01

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

• Very good answer. What I wanted to propose are probability generating functions. In this case for $X_i$: $g_i(s) = E[s^{X_i}] = (1-p_i) + p_is$ and the generating function of the sum is $g(s) = \prod_{i=1}^n g_i(s)$. Taking the nth derivative wrt $s$ evaluated at $0$ gives all the point probabilites. Finally sum them up for the cdf. I just can not offer code for this, so your answer can be readily applied. – Ric Feb 7 '13 at 8:50
• @Richard The pgf approach works just fine--in fact, I used that method (in Mathematica, where it's a one-liner) to verify my example here. But under the hood, the product is carried out as a convolution, so it amounts to the same solution. BTW, there's no need to take a derivative: you only need to expand the product as a polynomial in $s$ and pick out the probabilities from the coefficients. – whuber Feb 7 '13 at 16:38
• yes, of course the derivatives are just there to extract the coefficients. And yes - in R, it has to be the way you presented it. I just thought I write down the pgf approach for its mathematical beauty ;) – Ric Feb 7 '13 at 17:05

Tools like the binomial assume equal probabilities because that greatly simplifys things, so you cannot simplify. If you need an exact answer then you will need to look at every possibly way that your event (at least X events occure), calculate the probability of each possibility then sum those up. There are tools that will generate every possible combination for you so that you can sum them up (the combn function in R is one) so this is doable.

If the set of probabilities is fixed, so you don't need to recompute every time and can just precompute the chances, and you can live with a good approximation then I would suggest simulation. Randomly generate an outcome with the given set of probabilities, repeat this a bunch of times (10,000 or even 100,000) and then just calculate how often you had X or more events. This would take only a few lines in R (and probably be similar in other tools as well).

There is a similar question here, see my example there (using R) for the more exact calculation.

• Thanks! That's very helpful. In practice, I'll need to quickly test to see what the new probability will be when a NEW event (with known probability) is added to the list. – robamaton Feb 6 '13 at 20:41
• If you know the probabilities for a given set of events and you add 1 more possible event, then you can work the new probability out as either X or more events occured in the original set AND the new event did not occur OR $X-1$ or more events occured in the original set AND the new event did occur (total of X events). If you have the probabilities pre-computed for the original set then with 2 multiplications and an addition you will have the new probabilities. – Greg Snow Feb 6 '13 at 20:58
• Perfect. I figured that part would be relatively straightforward, and I'm glad that the preprocessing can make this a quick calculation. – robamaton Feb 6 '13 at 21:00