If 20 independent Bernoulli trials are carried out each with a different probability of success and therefore failure. What is the probability that exactly n of the 20 trials was successful?
Is there a better way of calculating these probabilities rather than simply summing together the combinations of success and failure probabilities?
The distribution you are asking about is called the Poisson Binomial distribution, with rather complicated pmf (see Wikipedia for broader description)
$$ \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) $$
Generally, the problem is that using it directly for a larger number of trials would be very slow. There are also other methods of calculating the pmf, e.g. recursive formulas, but they are numerically unstable. The easiest way around those problems are approximation methods (described e.g. by Hong, 2013). If we define
$$ \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)$$
then we can approximate pmf with Poisson distribution via law of small numbers or Le Cams theorem
$$ \Pr(X = x) \approx \frac{\mu^x \exp(-\mu)}{x!} $$
but it sees that generally Binomial approximation behaves better (Choi and Xia, 2002)
$$ \Pr(X = x) \approx \mathrm{Binom} \left( n, \frac{\mu}{n} \right)$$
you can use Normal approximation
$$ f(x) \approx \phi \left( \frac{x + 0.5 - \mu}{ \sigma} \right) $$
or cdf can be approximated using so-called refined Normal approximation (Volkova, 1996)
$$ F(x) \approx \max\left(0, \ g \left( \frac{x + 0.5 - \mu}{ \sigma} \right) \right)$$
where $g(x) = \Phi(x) + \gamma(1-x^2) \frac{\phi(x)}{6}$.
Another alternative is of course a Monte Carlo simulation.
Simple dpbinom R function would be
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
}
}
Most of the methods (and more) are also implemented in R poibin package.
Chen, L.H.Y. (1974). On the Convergence of Poisson Binomial to Poisson Distributions. The Annals of Probability, 2(1), 178-180.
Chen, S.X. and Liu, J.S. (1997). Statistical applications of the Poisson-Binomial and conditional Bernoulli distributions. Statistica Sinica 7, 875-892.
Chen, S.X. (1993). Poisson-Binomial distribution, conditional Bernoulli distribution and maximum entropy. Technical Report. Department of Statistics, Harvard University.
Chen, X.H., Dempster, A.P. and Liu, J.S. (1994). Weighted finite population sampling to maximize entropy. Biometrika 81, 457-469.
Wang, Y.H. (1993). On the number of successes in independent trials. Statistica Sinica 3(2): 295-312.
Hong, Y. (2013). On computing the distribution function for the Poisson binomial distribution. Computational Statistics & Data Analysis, 59, 41-51.
Volkova, A. Y. (1996). A refinement of the central limit theorem for sums of independent random indicators. Theory of Probability and its Applications 40, 791-794.
Choi, K.P. and Xia, A. (2002). Approximating the number of successes in independent trials: Binomial versus Poisson. The Annals of Applied Probability, 14(4), 1139-1148.
Le Cam, L. (1960). An Approximation Theorem for the Poisson Binomial Distribution. Pacific Journal of Mathematics 10(4), 1181–1197.
One approach is to use generating functions. The solution to your problem is the coefficient $x^n$ in the polynomial
$$\prod_{i=1}^{20}(p_ix + 1-p_i).$$
This is the dynamic programming equivalent (quadratic time in the number of Bernoulli variables) of doing the summation in the Poisson Binomial distribution from Tim's answer (which would be exponential time).
Here's the Python code of the quadratic-time dynamic programming algorithm for a given $n$ and $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
Its numerical precision could be improved by implementing the Kahan summation alagorithm, but there's probably very little benefit since adjacent entries in the running totals (which are the addends) are usually not that different in magnitude.
