Risk of extinction of Schrödinger's cats

Question

I am interested in how uncertainty can be accounted for when considering the risk of extinction of a species. Forgive me for extending a rather tired thought experiment, but at least it's familiar territory and I hope it illustrates what I am trying to understand.

Let's say that Schrödinger was not satisfied with killing and not killing only one cat, so he went out and collected the last 15 remaining Himalayan Snow-Cats. He put each one in a box with a vial of poison, hammer, and a trigger device for releasing the hammer. Each trigger device has a known probability of releasing the hammer within any given hour, and the poison will take 5 hours to kill a cat (things got a bit trippy when he used radioactive decay, so this time he steers clear of quantum mechanics). After one hour, Schrödinger receives a notice from the ethics board telling him he's nuts and ordering him to release the cats immediately. He pushes a button which opens a cat-flap at the back of each box, thus releasing them all back into the wild. By the time the humane society check the boxes, all the cats have bolted. Before anyone can restrain him, Schrödinger detonates the remaining triggers, so nobody knows how many of the cats were poisoned.

The probability that each of the cats is poisoned is as follows:

1. 0.17
2. 0.46
3. 0.62
4. 0.08
5. 0.40
6. 0.76
7. 0.03
8. 0.47
9. 0.53
10. 0.32
11. 0.21
12. 0.85
13. 0.31
14. 0.38
15. 0.69

After 5 hours have elapsed and all poisoned cats have died, what is the probability that there are:

• (A) More than 10 Himalayan Snow-Cats still alive
• (B) 10 or less still alive
• (C) 5 or less still alive
• (D) 2 or less still alive

Answer 1

Let $p_k, k=1, \ldots, K=15$ be the survival probabilities for individual snow cats.

1. Initialize the number of cats accounted for so far $k \leftarrow 0$, the vector of the binomial survival probabilities $(\pi_0^{(0)}, \pi_1^{(0)}, \ldots, \pi_K^{(0)}) \leftarrow (1, 0, \ldots, 0, 0)$
2. While there still are unaccounted cats, $k \le K$, repeat steps 3-5:
3. Increase $k \leftarrow k+1$
4. Update the 0 outcome (death): $\pi_j^{(k)} \leftarrow \pi_j^{(k-1)} (1-p_k), j=0, \ldots, K$
5. Update the 1 outcome (survival): $\pi_{j+1}^{(k)} \leftarrow \pi_{j+1}^{(k)} + \pi_j^{(k-1)} p_k, j=0, \ldots, K$

I would probably be paranoid about accounting for everything, and make sure that my probabilities still sum up to 1 at each iteration after step 5.

My results are:

4.091e-08  
2.647e-06  
.00006039  
.00069791  
.00479963  
.02141555  
.06519699  
.13945642  
.21276277  
.23238555  
.18045155  
.09775029  
.03565983  
.00823336  
.0010688  
.00005826  

The sum of the first three is the prob of being critically endangered, the last 5, the prob of least concern, etc.

Answer 2

It's unclear what "risk of extinction category" means, but it appears the question asks to compute the distribution of the sum of 15 independent binomial variates having the given expectations. This is a convolution and it's most efficiently done with the Fast Fourier Transform.

Here is an example in R, whose convolve function uses the FFT:

x <- c(0.17,0.46,0.62,0.08,0.40,0.76,0.03,0.47,0.53,0.32,0.21,0.85,0.31,0.38,0.69)
z <- 1
for (u in sort(x)) z <- convolve(z, c(u, 1-u), type="open")
z

[1] 5.826e-05 1.069e-03 8.233e-03 3.566e-02 9.775e-02 1.805e-01 2.324e-01 2.128e-01

[9] 1.395e-01 6.520e-02 2.142e-02 4.800e-03 6.979e-04 6.039e-05 2.647e-06 4.091e-08

As a check, in Mathematica the same results were obtained with

Product[1 - z + z t, {z, x}] // Expand

$4.091095\times 10^{-8} t^{15}+2.647052\times 10^{-6} t^{14}+\cdots+0.0010688 t+0.0000582614$

Now, for instance, the chance of $10$ or fewer poisonings is computed in R as

sum(z[1:11])

[1] 0.9944

Edit

Using R's convolve function is inefficient for larger problems, because it repeatedly performs an FFT and its inverse. It turns out that the direct algorithm for convolution--as described by StasK--is plenty fast enough. Here is an R implementation.

convolve.binomial <- function(p) {
  # p is a vector of probabilities of Bernoulli distributions.
  # The convolution of these distributions is returned as a vector
  # `z` where z[i] is the probability of i-1, i=1, 2, ..., length(p)+1.
  n <- length(p) + 1
  z <- c(1, rep(0, n-1))
  for (p in p) z <- (1-p)*z + p*c(0, z[-n])
  return(z)
}

(Thanks to an anonymous editor suggesting improvements that clarify the algorithm.)

This takes $O(n^2)$ time for $n$ distributions--the quadratic behavior is not good--but it's still quite fast. As an example, let's generate 10,000 random probabilities (instead of using 15 given ones) and form the convolution of the corresponding Bernoulli distributions:

x <- runif(10000)
system.time(y <- convolve.binomial(x))

This still takes less than 3 seconds.