目前尚不清楚“灭绝风险类别”的含义，但似乎问题要求计算具有给定期望的 15 个独立二项式变量之和的分布。这是一个卷积，使用快速傅里叶变换最有效。

这是 中的一个示例R，其convolve函数使用 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

作为检查，在Mathematica中获得了相同的结果

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$

现在，例如，机会$10$或更少的中毒计算R为

sum(z[1:11])

[1] 0.9944

编辑

对于较大的问题，使用R'convolve函数效率低下，因为它反复执行 FFT 及其逆运算。事实证明，直接卷积算法——正如 StasK 所描述的——足够快。这是一个R实现。

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

（感谢一位匿名编辑提出的改进以阐明算法。）

这需要$O(n^2)$的时间$n$分布——二次行为不好——但它仍然相当快。例如，让我们生成 10,000 个随机概率（而不是使用 15 个给定的概率）并形成相应伯努利分布的卷积：

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

这仍然需要不到 3 秒的时间。

让$p_k, k=1, \ldots, K=15$是单个雪猫的生存概率。

1. 初始化到目前为止占猫的数量$k \leftarrow 0$, 二项式生存概率的向量$(\pi_0^{(0)}, \pi_1^{(0)}, \ldots, \pi_K^{(0)}) \leftarrow (1, 0, \ldots, 0, 0)$
2. 虽然仍有下落不明的猫，$k \le K$，重复步骤 3-5：
3. 增加$k \leftarrow k+1$
4. 更新 0 结果（死亡）：$\pi_j^{(k)} \leftarrow \pi_j^{(k-1)} (1-p_k), j=0, \ldots, K$
5. 更新 1 结果（生存）：$\pi_{j+1}^{(k)} \leftarrow \pi_{j+1}^{(k)} + \pi_j^{(k-1)} p_k, j=0, \ldots, K$

我可能会偏执于考虑所有事情，并确保在第 5 步之后的每次迭代中我的概率总和仍为 1。

我的结果是：

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

前三个之和是极度濒危的概率，后五个之和是最不关心的概率，等等。