您所询问的分布称为Poisson Binomial 分布，具有相当复杂的 pmf（有关更广泛的描述，请参见 Wikipedia）

$$\Pr(X=x) = \sum\limits_{A\in F_x} \prod\limits_{i\in A} p_i \prod\limits_{j\in A^c} (1-p_j)$$

通常，问题是直接使用它进行大量试验会非常慢。还有其他计算 pmf 的方法，例如递归公式，但它们在数值上不稳定。解决这些问题的最简单方法是近似方法（例如由Hong, 2013描述）。如果我们定义

$$\mu = \sum_{i=1}^n p_i$$

$$\sigma = \sqrt{ \sum_{i=1}^n p_i(1-p_i) }$$

$$\gamma = \sigma^{-3} \sum_{i=1}^n p_i (1 - p_i) (1 - 2p_i)$$

然后我们可以通过小数定律或Le Cams 定理用泊松分布近似 pmf

$$\Pr(X = x) \approx \frac{\mu^x \exp(-\mu)}{x!}$$

但它发现通常二项式近似表现更好（Choi and Xia, 2002）

$$\Pr(X = x) \approx \mathrm{Binom} \left( n, \frac{\mu}{n} \right)$$

您可以使用正态近似

$$f(x) \approx \phi \left( \frac{x + 0.5 - \mu}{ \sigma} \right)$$

或 cdf 可以使用所谓的精细正态近似来近似（Volkova，1996）

$$F(x) \approx \max\left(0, \ g \left( \frac{x + 0.5 - \mu}{ \sigma} \right) \right)$$

其中。$g(x) = \Phi(x) + \gamma(1-x^2) \frac{\phi(x)}{6}$

另一种选择当然是蒙特卡罗模拟。

简单dpbinom的 R 函数将是

dpbinom <- function(x, prob, log = FALSE,
                    method = c("MC", "PA", "NA", "BA"),
                    nsim = 1e4) {

  stopifnot(all(prob >= 0 & prob <= 1))
  method <- match.arg(method)

  if (method == "PA") {
    # poisson
    dpois(x, sum(prob), log)
  } else if (method == "NA") {
    # normal
    dnorm(x, sum(prob), sqrt(sum(prob*(1-prob))), log)
  } else if (method == "BA") {
    # binomial
    dbinom(x, length(prob), mean(prob), log)
  } else {
    # monte carlo
    tmp <- table(colSums(replicate(nsim, rbinom(length(prob), 1, prob))))
    tmp <- tmp/sum(tmp)
    p <- as.numeric(tmp[as.character(x)])
    p[is.na(p)] <- 0
        
    if (log) log(p)
    else p 
  }
}

大多数方法（以及更多）也在 R poibin包中实现。

---

陈 LHY (1974)。关于泊松二项式到泊松分布的收敛性。概率年鉴，2（1），178-180。

Chen, SX 和 Liu, JS (1997)。泊松二项分布和条件伯努利分布的统计应用。中央统计 7, 875-892。

陈 SX (1993)。泊松二项分布、条件伯努利分布和最大熵。技术报告。哈佛大学统计系。

Chen, XH, Dempster, AP 和 Liu, JS (1994)。加权有限总体抽样以最大化熵。生物计量学 81, 457-469。

王永辉 (1993)。关于独立试验的成功次数。统计学 3(2): 295-312。

洪，Y.（2013）。关于计算泊松二项分布的分布函数。计算统计与数据分析，59、41-51。

沃尔科娃，AY（1996 年）。对独立随机指标之和的中心极限定理的改进。概率论及其应用 40, 791-794。

Choi, KP 和 Xia, A. (2002)。近似独立试验的成功次数：二项式与泊松。应用概率年鉴，14（4），1139-1148。

Le Cam, L. (1960)。泊松二项分布的近似定理。太平洋数学杂志 10（4），1181-1197。

一种方法是使用生成函数。您的问题的解决方案是多项式中$x^n$

$$\prod_{i=1}^{20}(p_ix + 1-p_i).$$

这是根据蒂姆的答案（这将是指数时间）在泊松二项式分布中求和的动态规划等价物（伯努利变量数量的二次时间）。

和的二次时间动态规划算法的 Python 代码：$n$$p$

import numpy as np

def calculated_probability(ps, n):
    total = np.zeros((ps.shape[0] + 1,))
    total[0] = 1.0
    for p in ps:
        total = p * np.roll(total, 1) + (1 - p) * total
    return total[n]
    
rng = np.random.default_rng(12345)
ps = rng.uniform(size=10000)
print(calculated_probability(ps, 5000))  # 0.008196669065619853

通过实施Kahan summation algorithm可以提高其数值精度，但可能没有什么好处，因为运行总计中的相邻条目（即加数）通常在幅度上没有太大差异。