How can I calculate the probability for multiple trials with different probabilities?

So for a bit of recreational spreadsheet creation I've been pondering probabilities.

Let's say I have two types of trial, A and B, with a 40% and 25% chance of success respectively. How can I calculate the probability, for instance, that I get at least two successes if I run trial A three times and B twice?

As I understand it, I can use a binomial distribution to calculate each one individually: 35.2% for at least two successes from A, or 6.25% from B. In addition, I would reckon the probability of getting at least one success in the combined set to be 87.85%, as I calculate a 12.15% chance of getting exactly zero successes on all five trials.

I'd greatly appreciate a little help working this out, because it's been bothering me for a couple of days.

• you need a convolution of Binomials, I doubt that analytical expression exists. I'm afraid you have to approximate or simulate, see further en.wikipedia.org/wiki/Convolution_of_probability_distributions Commented May 20, 2020 at 14:52
• @Aksakal The answer can be expressed as a sum of hypergeometric functions (which are polynomials in the odds ratio).
– whuber
Commented May 20, 2020 at 15:23
• Seems that I've accidentally stumbled across something a lot more complex than it appeared on the surface... Commented May 20, 2020 at 15:29
• @whuber Could you explain in bit more detail as to how the answer to this problem can be expressed as a sum of hypergeometrics? I am quite a novice to this area of study. Commented May 20, 2020 at 16:15

Suppose you conduct the first trial $$n$$ times and (independently) second one $$m$$ times, and that the chances of success are $$p$$ and $$q$$ respectively. Let $$A$$ be the total number of successes in the first instance, $$B$$ the total in the second, and $$X=A+B$$ be the total number of successes. Obviously $$X$$ is an integer between $$0$$ and $$m+n$$ (inclusive). For any such integer $$x,$$ let's find an expression for the chance that $$X=x.$$

One such expression exploits the probability axiom that says the chance of an event is the sum of the chances of mutually disjoint events of which it is comprised. Here, the event $$X=x$$ is comprised of the events $$A=a, B=x-a$$ where $$a$$ ranges over all possible counts (of successes of $$A$$).

The independence of $$A$$ and $$B$$ implies the chance an event $$A=a,B=x-a$$ is the product of the component chances. Since $$A$$ and $$B$$ have Binomial distributions, we have immediately

$$\Pr(A=a,B=x-a) = \left(\binom{n}{a} p^a(1-p)^{n-a}\right)\left(\binom{m}{x-a} q^{x-a}(1-q)^{m-(x-a)}\right).$$

Summing these over $$a$$ and doing a little algebraic simplification yields

\eqalign{\Pr(X=x) &= (1-p)^n q^x(1-q)^{m-x}\,\sum_{a=0}^x \binom{n}{a}\binom{m}{x-a} \left(\frac{p(1-q)}{(1-p)q}\right)^a \\ &= \phi^x (1-p)^n (1-q)^{m}\,\sum_{a=0}^x \binom{n}{a}\binom{m}{x-a} t^a \\ &= \phi^x (1-p)^n (1-q)^{m}\binom{m}{x}\,\sum_{a=0}^\infty \frac{(-n)^a (-x)^a }{a! (m-x+1)^a} (-1)^{2a}t^a \\ &= \phi^x (1-p)^n(1-q)^{m}\binom{m}{x}\,_2F_1(-n,-x;m-x+1;t) }

where

$$\phi = \frac{q}{1-q}$$

is the odds for $$B,$$

$$t = \frac{p}{1-p}\,/\,\frac{q}{1-q}$$

is the odds ratio for $$A$$ relative to $$B,$$

$$(z)^s = z(z+1)\cdots(z+s-1)$$

is the "rising factorial" (or Pochhammer symbol), and $$\,_2F_1$$ is the Riemann hypergeometric function (which, in this case, obviously reduces to a polynomial in $$t$$ of degree no greater than $$x$$).

Find the chance of the event $$X\ge x$$ (as in the question) by summing over the individual possibilities of $$x$$ or, when $$x$$ is small, by computing the chance of its complement,

$$\Pr(X \ge x) = 1 - \Pr(X \lt x) = 1 - \sum_{y=0}^{x-1} \Pr(X = y).$$

For tiny values of $$x$$ this won't be too bad; for larger values, you will want a good software library for computing values of the hypergeometric function.

Remarks

Convolution of the two binomial distributions (using the Fast Fourier Transform) is an attractive option for precise calculation.

When both of $$np+mq$$ and $$n(1-p)+m(1-q)$$ are not small (exceeding $$5$$ is often considered ok), the Normal approximation to the Binomial distributions will give a good approximation. Specifically, the approximating Normal distribution will have mean

$$\mu= np + mq,$$

variance

$$\sigma^2 = np(1-p) + mq(1-q),$$

and the chance is therefore approximated (using a continuity correction) by

$$\Pr(X \ge x) \approx \Phi\left(\frac{\mu - x + 1/2}{\sigma}\right)$$

where $$\Phi$$ is the CDF of the standard Normal distribution. If you're brave, you can also approximate the individual probabilities as

\eqalign{ \Pr(X = x) &= \Pr(X \ge x) - \Pr(X \ge x+1) \\ &\approx \Phi\left(\frac{\mu-x+1/2}{\sigma}\right) - \Phi\left(\frac{\mu-x-1/2}{\sigma}\right).}

As an example, with $$n=6,$$ $$p=0.40,$$ $$m=10,$$ and $$q=0.25$$ (the chances in the question, with the minimal numbers of trials for the approximation to hold), a simulation of 100,000 values of $$X$$ (shown by the line heights) is pretty well reproduced by the approximation (shown by the dots):

This R code produced the figure.

n <- 6
m <- 10
p <- 0.4
q <- 0.25
#
# Simulate X.
#
n.sim <- 1e5
A <- rbinom(n.sim, n, p)
B <- rbinom(n.sim, m, q)
X <- A+B
#
# Plot the simulation.
#
plot(0:(n+m), tabulate(X+1, n+m+1)/n.sim, type="h", ylab="Relative frequency", xlab="x")
#
# Plot the Normal approximation.
#
f <- function(x, n, p, m, q) {
mu <- n * p + m * q
sigma <- sqrt(n * p * (1-p) + m * q * (1-q))
pnorm((x + 1/2 - mu) / sigma) - pnorm((x-1 + 1/2 - mu) / sigma)
}
points(0:(n+m), f(0:(n+m), n, p, m, q), pch=21, bg="#e0000080")

• (+1) I was expecting to see your answer :) Commented May 20, 2020 at 21:44

The sum of independent non-identically distributed Bernoulli trials is called a Poisson-Binomial distribution. The probabilities can be easily calculated with an R package called poibin.

For the example in the OP's description (3 trials with 40% chance of success and 2 trials with 25% chance of success), the following code will find the pmf:

library(poibin)

p_success <- c(rep(0.4,3),rep(0.25,2))
kk <- 0:5

pmf <- dpoibin(pp=p_success,kk=kk)

pmf


Here is the output, which matches what the OP found for the probability of zero successes and gives the probabilities for 1-5 successes as well:

[1] 0.1215 0.3240 0.3375 0.1710 0.0420 0.0040


For at least two successes, we can calculate the probability of no success and only one success, sum them and subtract from $$1$$. You've already done the first. For one success, the cases are

• 1 success for A, 0 success for B $$\rightarrow {3\choose 1}p_a(1-p_a)^2(1-p_b)^2$$
• 0 success for A, 1 success for B $$\rightarrow (1-p_a)^3{2\choose 1}p_b(1-p_b)$$

This was a simple situation and the general case is harder as pointed out int the comments.

• Thank you. Unfortunately I was hoping the general case wouldn't be so problematic, but it looks like I might have to rethink how I'm structuring my data to avoid it, else leap into the deep end. Commented May 20, 2020 at 15:33