Which test should I use?

Question

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.

Answer 1

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

Answer 2

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.