# Probability distribution for different probabilities

If I wanted to get the probability of 9 successes in 16 trials with each trial having a probability of 0.6 I could use a binomial distribution. What could I use if each of the 16 trials has a different probability of success?

• @whuber In your explanation of Normal approximation, the calculations of mean and standard deviation is different with the description in Wikipedia. In Wiki, the mean is np, and the standard deviation is np(1-p). So, in this problem, for the Normal approximation of varying probability of success in binomial distribution, the mean is p1+p2+p3+p4+p5+...+pi, and the variance is p1(1-p1)+p2(1-p2)+...+pi(1-pi). Am I right? – David Nov 25 '13 at 10:34
• See Wikipedia on the Poisson binomial distribution. Also a search term that turns up a few hits here. – Glen_b Dec 2 '13 at 1:56
• @David When all the $p_i$ are equal to a common value $p$, then $p_1+p_2+\cdots+p_n = np$ and $p_1(1-p_1)+\cdots+p_n(1-p_n)=np(1-p)$, showing that the Wikipedia description you refer to is just a special case. – whuber Oct 24 '14 at 16:09
• – user83346 Aug 27 '15 at 8:16

This is the sum of 16 (presumably independent) Binomial trials. The assumption of independence allows us to multiply probabilities. Whence, after two trials with probabilities $p_1$ and $p_2$ of success the chance of success on both trials is $p_1 p_2$, the chance of no successes is $(1-p_1)(1-p_2)$, and the chance of one success is $p_1(1-p_2) + (1-p_1)p_2$. That last expression owes its validity to the fact that the two ways of getting exactly one success are mutually exclusive: at most one of them can actually happen. That means their probabilities add.

By means of these two rules--independent probabilities multiply and mutually exclusive ones add--you can work out the answers for, say, 16 trials with probabilities $p_1, \ldots, p_{16}$. To do so, you need to account for all the ways of obtaining each given number of successes (such as 9). There are $\binom{16}{9} = 11440$ ways to achieve 9 successes. One of them, for example, occurs when trials 1, 2, 4, 5, 6, 11, 12, 14, and 15 are successes and the others are failures. The successes had probabilities $p_1, p_2, p_4, p_5, p_6, p_{11}, p_{12}, p_{14},$ and $p_{15}$ and the failures had probabilities $1-p_3, 1-p_7, \ldots, 1-p_{13}, 1-p_{16}$. Multiplying these 16 numbers gives the chance of this particular sequence of outcomes. Summing this number along with the 11,439 remaining such numbers gives the answer.

Of course you would use a computer.

With many more than 16 trials, there is a need to approximate the distribution. Provided none of the probabilities $p_i$ and $1-p_i$ get too small, a Normal approximation tends to work well. With this method you note that the expectation of the sum of $n$ trials is $\mu = p_1 + p_2 + \cdots + p_n$ and (because the trials are independent) the variance is $\sigma^2 = p_1(1-p_1) + p_2(1-p_2) + \cdots + p_n(1-p_n)$. You then pretend the distribution of sums is Normal with mean $\mu$ and standard deviation $\sigma$. The answers tend to be good for computing probabilities corresponding to a proportion of successes that differs from $\mu$ by no more than a few multiples of $\sigma$. As $n$ grows large this approximation gets ever more accurate and works for even larger multiples of $\sigma$ away from $\mu$.

• Computer scientists call these "Poisson trials" to distinguish them from Bernoulli trials. In addition to Central Limit Theorem approximations, there are also good tail bounds available. Here's one. Google searches on "Chernoff bounds for Poisson trials" will turn up the results you might find in a typical CS treatment. – cardinal Apr 13 '11 at 17:10
• @Cardinal That nomenclature is interesting. It would be valid for very small $p_i$, but otherwise seems misleading, because the distribution otherwise is not well approximated by Poisson distributions. (There's another discussion on CV about this question, where "16" is replaced by 10,000 and we do examine tail probabilities, but I haven't been able to find it again.) – whuber Apr 13 '11 at 17:23
• yes, I agree on the name. I found it a bit odd when I first encountered it. I've given it here more as a useful term for searching. It appears computer scientists consider these probabilities often in dealing with certain algorithms. I'd be interested in reading that other question if you happen to find it. Is it this one perhaps? – cardinal Apr 13 '11 at 17:26
• @cardinal is right that we "CS folks" call them Poisson trials. in fact for this case, a standard Chernoff-Hoeffding bound will give you exactly the bound the OP is asking for. – Suresh Venkatasubramanian Apr 13 '11 at 22:50
• As per the comment by @David yesterday, there is something wrong with your statement of the Normal approximating mean as $$\mu = (p_1 + p_2 + \cdots + p_n)/n$$ We are summing 16 Bernoulli rvs, each of which can take value 0 or 1, so the sum will have domain of support from 0 to 16, not between 0 and 1. Worth checking your sd too. – wolfies Nov 26 '13 at 14:12

One alternative to @whuber's normal approximation is to use "mixing" probabilities, or a hierarchical model. This would apply when the $p_i$ are similar in some way, and you can model this by a probability distribution $p_i\sim Dist(\theta)$ with a density function of $g(p|\theta)$ indexed by some parameter $\theta$. you get a integral equation:

$$Pr(s=9|n=16,\theta)={16 \choose 9}\int_{0}^{1} p^{9}(1-p)^{7}g(p|\theta)dp$$

The binomial probability comes from setting $g(p|\theta)=\delta(p-\theta)$, the normal approximation comes from (I think) setting $g(p|\theta)=g(p|\mu,\sigma)=\frac{1}{\sigma}\phi\left(\frac{p-\mu}{\sigma}\right)$ (with $\mu$ and $\sigma$ as defined in @whuber's answer) and then noting the "tails" of this PDF fall off sharply around the peak.

You could also use a beta distribution, which would lead to a simple analytic form, and which need not suffer from the "small p" problem that the normal approximation does - as beta is quite flexible. Using a $beta(\alpha,\beta)$ distribution with $\alpha,\beta$ set by the solutions to the following equations (this is the "mimimum KL divergence" estimates):

$$\psi(\alpha)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[p_{i}]$$ $$\psi(\beta)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[1-p_{i}]$$

Where $\psi(.)$ is the digamma function - closely related to harmonic series.

We get the "beta-binomial" compound distribution:

$${16 \choose 9}\frac{1}{B(\alpha,\beta)}\int_{0}^{1} p^{9+\alpha-1}(1-p)^{7+\beta-1}dp ={16 \choose 9}\frac{B(\alpha+9,\beta+7)}{B(\alpha,\beta)}$$

This distribution converges towards a normal distribution in the case that @whuber points out - but should give reasonable answers for small $n$ and skewed $p_i$ - but not for multimodal $p_i$, as beta distribution only has one peak. But you can easily fix this, by simply using $M$ beta distributions for the $M$ modes. You break up the integral from $0<p<1$ into $M$ pieces so that each piece has a unique mode (and enough data to estimate parameters), and fit a beta distribution within each piece. then add up the results, noting that making the change of variables $p=\frac{x-L}{U-L}$ for $L<x<U$ the beta integral transforms to:

$$B(\alpha,\beta)=\int_{L}^{U}\frac{(x-L)^{\alpha-1}(U-x)^{\beta-1}}{(U-L)^{\alpha+\beta-1}}dx$$

• +1 This answer contains some interesting and clever suggestions. The last one looks particularly flexible and powerful. – whuber Nov 13 '13 at 21:37
• Just to take something very simple and concrete, suppose (i) $p_i = \frac{i}{17}$ and (ii) $p_i = \sqrt{i}/17$, for $i = 1$ to 16. What would be the solution to your $\alpha$ and $\beta$ estimates, and thus your estimates for $P(X=9)$ given $n= 16$, as per the OP's problem? – wolfies Nov 26 '13 at 14:06
• Great answer and proposal, especially the beta! It'd be cool to see this answer written in its general form with $n$ and $s$. – pglpm Oct 16 '19 at 9:45

Let $X_i$ ~ $Bernoulli(p_i)$ with probability generating function (pgf):

$$\text{pgf} = E[t^{X_i}] = 1 - p_i (1-t)$$

Let $S = \sum_{i=1}^n X_i$ denote the sum of $n$ such independent random variables. Then, the pgf for the sum $S$ of $n=16$ such variables is:

\begin{align*}\displaystyle \text{pgfS} &= E[t^S] \\&= E[t^{X_1}] E[t^{X_2}] \dots E[t^{X_{16}}] \text{ (... by independence)} \\ &= \prod _{i=1}^{16} \left(1-p_i(1-t) \right)\end{align*}

We seek $P(S=9)$, which is:

$$\frac{1}{9!}\frac{d^9 \text{pgfS}}{dt^9}|_{t=0}$$

ALL DONE. This produces the exact symbolic solution as a function of the $p_i$. The answer is rather long to print on screen, but it is entirely tractable, and takes less than $\frac{1}{100}$th of a second to evaluate using Mathematica on my computer.

Examples

If $p_i = \frac{i}{17}, i= 1 \text{ to } 16$, then: $P(S=9) = \frac{9647941854334808184}{48661191875666868481} = 0.198268 \dots$

If $p_i = \frac{\sqrt{i}}{17}, i= 1 \text{ to } 16$, then: $P(S=9) = 0.000228613 \dots$

More than 16 trials?

With more than 16 trials, there is no need to approximate the distribution. The above exact method works just as easily for examples with say $n = 50$ or $n = 100$. For instance, when $n = 50$, it takes less than $\frac{1}{10}$th of second to evaluate the entire pmf (i.e. at every value $s = 0, 1, \dots, 50$) using the code below.

Mathematica code

Given a vector of $p_i$ values, say:

n = 16;   pvals = Table[Subscript[p, i] -> i/(n+1), {i, n}];


... here is some Mathematica code to do everything required:

pgfS = Expand[ Product[1-(1-t)Subscript[p,i], {i, n}] /. pvals];
D[pgfS, {t, 9}]/9! /. t -> 0  // N


0.198268

To derive the entire pmf:

Table[D[pgfS, {t,s}]/s! /. t -> 0 // N, {s, 0, n}]


... or use the even neater and faster (thanks to a suggestion from Ray Koopman below):

CoefficientList[pgfS, t] // N


For an example with $n = 1000$, it takes just 1 second to calculate pgfS, and then 0.002 seconds to derive the entire pmf using CoefficientList, so it is extremely efficient.

• It can be even simpler. With[{p = Range@16/17}, N@Coefficient[Times@@(1-p+p*t),t,9]] gives the probability of 9 successes, and With[{p = Range@16/17}, N@CoefficientList[Times@@(1-p+p*t),t]] gives the probabilities of 0,...,16 successes. – Ray Koopman Dec 1 '13 at 20:24
• @RayKoopman That is cool. The Table for the $p$-values is intentional to allow for more general forms not suitable with Range. Your use of CoefficientList is very nice! I've added an Expand to the code above which speeds the direct approach up enormously. Even so, CoefficientList is even faster than a ParallelTable. It does not make much difference for $n$ under 50 (both approaches take just a tiny fraction of a second either way to generate the entire pmf), but your CoefficientList will also be a real practical advantage when n is really large. – wolfies Dec 2 '13 at 15:46

@wolfies comment, and my attempt at a response to it revealed an important problem with my other answer, which I will discuss later.

Specific Case (n=16)

There is a fairly efficient way to code up the full distribution by using the "trick" of using base 2 (binary) numbers in the calculation. It only requires 4 lines of R code to get the full distribution of $Y=\sum_{i=1}^{n} Z_i$ where $Pr(Z_i=1)=p_i$. Basically, there are a total of $2^n$ choices of the vector $z=(z_1,\dots,z_n)$ that the binary variables $Z_i$ could take. Now suppose we number each distinct choice from $1$ up to $2^n$. This on its own is nothing special, but now suppose that we represent the "choice number" using base 2 arithmetic. Now take $n=3$ so I can write down all the choices so there are $2^3=8$ choices. Then $1,2,3,4,5,6,7,8$ in "ordinary numbers" becomes $1,10,11,100,101,110,111,1000$ in "binary numbers". Now suppose we write these as four digit numbers, then we have $0001,0010,0011,0100,0101,0110,0111,1000$. Now look at the last $3$ digits of each number - $001$ can be thought of as $(Z_1=0,Z_2=0,Z_3=1)\implies Y=1$, etc. Counting in binary form provides an efficient way to organise the summation. Fortunately, there is an R function which can do this binary conversion for us, called intToBits(x) and we convert the raw binary form into a numeric via as.numeric(intToBits(x)), then we will get a vector with $32$ elements, each element being the digit of the base 2 version of our number (read from right to left, not left to right). Using this trick combined with some other R vectorisations, we can calculate the probability that $y=9$ in 4 lines of R code:

exact_calc <- function(y,p){
n       <- length(p)
z       <- t(matrix(as.numeric(intToBits(1:2^n)),ncol=2^n))[,1:n] #don't need columns n+1,...,32 as these are always 0
pz      <- z%*%log(p/(1-p))+sum(log(1-p))
ydist   <- rowsum(exp(pz),rowSums(z))
return(ydist[y+1])
}


Plugging in the uniform case $p_i^{(1)}=\frac{i}{17}$ and the sqrt root case $p_i^{(2)}=\frac{\sqrt{i}}{17}$ gives a full distribution for y as:

$$\begin{array}{c|c}y & Pr(Y=y|p_i=\frac{i}{17}) & Pr(Y=y|p_i=\frac{\sqrt{i}}{17})\\ \hline 0 & 0.0000 & 0.0558 \\ 1 & 0.0000 & 0.1784 \\ 2 & 0.0003 & 0.2652 \\ 3 & 0.0026 & 0.2430 \\ 4 & 0.0139 & 0.1536 \\ 5 & 0.0491 & 0.0710 \\ 6 & 0.1181 & 0.0248 \\ 7 & 0.1983 & 0.0067 \\ 8 & 0.2353 & 0.0014 \\ 9 & 0.1983 & 0.0002 \\ 10 & 0.1181 & 0.0000 \\ 11 & 0.0491 & 0.0000 \\ 12 & 0.0139 & 0.0000 \\ 13 & 0.0026 & 0.0000 \\ 14 & 0.0003 & 0.0000 \\ 15 & 0.0000 & 0.0000 \\ 16 & 0.0000 & 0.0000 \\ \end{array}$$

So for the specific problem of $y$ successes in $16$ trials, the exact calculations are straight-forward. This also works for a number of probabilities up to about $n=20$ - beyond that you are likely to start to run into memory problems, and different computing tricks are needed.

Note that by applying my suggested "beta distribution" we get parameter estimates of $\alpha=\beta=1.3206$ and this gives a probability estimate that is nearly uniform in $y$, giving an approximate value of $pr(y=9)=0.06799\approx\frac{1}{17}$. This seems strange given that a density of a beta distribution with $\alpha=\beta=1.3206$ closely approximates the histogram of the $p_i$ values. What went wrong?

General Case

I will now discuss the more general case, and why my simple beta approximation failed. Basically, by writing $(y|n,p)\sim Binom(n,p)$ and then mixing over $p$ with another distribution $p\sim f(\theta)$ is actually making an important assumption - that we can approximate the actual probability with a single binomial probability - the only problem that remains is which value of $p$ to use. One way to see this is to use the mixing density which is discrete uniform over the actual $p_i$. So we replace the beta distribution $p\sim Beta(a,b)$ with a discrete density of $p\sim \sum_{i=1}^{16}w_i\delta(p-p_i)$. Then using the mixing approximation can be expressed in words as choose a $p_i$ value with probability $w_i$, and assume all bernoulli trials have this probability. Clearly, for such an approximation to work well, most of the $p_i$ values should be similar to each other. This basically means that for @wolfies uniform distribution of values, $p_i=\frac{i}{17}$ results in a woefully bad approximation when using the beta mixing distribution. This also explains why the approximation is much better for $p_i=\frac{\sqrt{i}}{17}$ - they are less spread out.

The mixing then uses the observed $p_i$ to average over all possible choices of a single $p$. Now because "mixing" is like a weighted average, it cannot possibly do any better than using the single best $p$. So if the $p_i$ are sufficiently spread out, there can be no single $p$ that could provide a good approximation to all $p_i$.

One thing I did say in my other answer was that it may be better to use a mixture of beta distributions over a restricted range - but this still won't help here because this is still mixing over a single $p$. What makes more sense is split the interval $(0,1)$ up into pieces and have a binomial within each piece. For example, we could choose $(0,0.1,0.2,\dots,0.9,1)$ as our splits and fit nine binomials within each $0.1$ range of probability. Basically, within each split, we would fit a simple approximation, such as using a binomial with probability equal to the average of the $p_i$ in that range. If we make the intervals small enough, the approximation becomes arbitrarily good. But note that all this does is leave us with having to deal with a sum of indpendent binomial trials with different probabilities, instead of Bernoulli trials. However, the previous part to this answer showed that we can do the exact calculations provided that the number of binomials is sufficiently small, say 10-15 or so.

To extend the bernoulli-based answer to a binomial-based one, we simply "re-interpret" what the $Z_i$ variables are. We simply state that $Z_i=I(X_i>0)$ - this reduces to the original bernoulli-based $Z_i$ but now says which binomials the successes are coming from. So the case $(Z_1=0,Z_2=0,Z_3=1)$ now means that all the "successes" come from the third binomial, and none from the first two.

Note that this is still "exponential" in that the number of calculations is something like $k^g$ where $g$ is the number of binomials, and $k$ is the group size - so you have $Y\approx\sum_{j=1}^{g}X_j$ where $X_j\sim Bin(k,p_j)$. But this is better than the $2^{gk}$ that you'd be dealing with by using bernoulli random variables. For example, suppose we split the $n=16$ probabilities into $g=4$ groups with $k=4$ probabilities in each group. This gives $4^4=256$ calculations, compared to $2^{16}=65536$

By choosing $g=10$ groups, and noting that the limit was about $n=20$ which is about $10^7$ cells, we can effectively use this method to increase the maximum $n$ to $n=50$.

If we make a cruder approximation, by lowering $g$, we will increase the "feasible" size for $n$. $g=5$ means that you can have an effective $n$ of about $125$. Beyond this the normal approximation should be extremely accurate.

• @momo - I think this is ok, as my answers are two different ways to approach the problem. This answer is not an edited version of my first one - it is just a different answer – probabilityislogic Dec 1 '13 at 21:41
• For a solution in R that is extremely efficient and handles much, much larger values of $n$, please see stats.stackexchange.com/a/41263. For instance, it solved this problem for $n=10^4$, giving the full distribution, in under three seconds. (A comparable Mathematica 9 solution--see @wolfies' answer--also performs well for smaller $n$ but could not complete the execution with such a large value of $n$.) – whuber Dec 2 '13 at 17:09

The (in general intractable) pmf is $$\Pr(S=k) = \sum_{\substack{A\subset\{1,\dots,n\}\\ |A|=k}} \left( \prod_{i\in A} p_i \right)\left(\prod_{j\in \{1,\dots,n\}\setminus A} (1-p_j) \right) \, .$$ R code:

p <- seq(1, 16) / 17
cat(p, "\n")
n <- length(p)
k <- 9
S <- seq(1, n)
A <- combn(S, k)
pr <- 0
for (i in 1:choose(n, k)) {
pr <- pr + exp(sum(log(p[A[,i]])) + sum(log(1 - p[setdiff(S, A[,i])])))
}
cat("Pr(S = ", k, ") = ", pr, "\n", sep = "")


For the $p_i$'s used in wolfies answer, we have:

Pr(S = 9) = 0.1982677


When $n$ grows, use a convolution.

• Doing that with R code was really helpful. Some of us are more concrete-thinkers and it greatly helps to have an operational version of the generating function. – DWin Dec 1 '13 at 18:24
• @DWin I provide efficient R code in the solution to the same problem (with different values of the $p_i$) at stats.stackexchange.com/a/41263. The problem here is solved in 0.00012 seconds total computation time (estimated by solving it 1000 times) compared to 0.53 seconds (estimated by solving it once) for this R code and 0.00058 seconds using Wolfies' Mathematica code (estimated by solving it 1000 times). – whuber Dec 2 '13 at 16:59
• So $P(S=k)$ would follow a Poisson-Binomial Distribution. – fccoelho Oct 8 '14 at 12:12
• +1 Very useful post in my attempt at answering this question. I was wondering if using logs is more of a cool mathematical formulation than a real need. I am not too concerned about running times... – Antoni Parellada Oct 1 '16 at 5:46