Let's say I have $5$ unfair coins. Each with an independent, known, probability of landing on heads. I flip each coin once. How can I find the probability that I get $3$ or more heads?


$P(H_1) = .38$

$P(H_2) = .18$

$P(H_3) = .71$

$P(H_4) = .66$

$P(H_5) = .29$

  • 1
    $\begingroup$ Now you have only 5 different coins, should you have many more coins, say more than 16, and given that your probabilities for heads for each of those coins is not too small, then you could use a computer and approximate with a normal distribution. For details, see this extensive answer on the subject $\endgroup$
    – JJJ
    Mar 16, 2018 at 3:16
  • 1
    $\begingroup$ Chances are if you have 5 unfair coins, then these coins would be distinct as well, and breaking it down into cases would be your best bet. To find the case for 3 heads, for example, you would need to add the 10 different possible combinations that can make three heads appear $\endgroup$
    – WaveX
    Mar 16, 2018 at 3:16

7 Answers 7


Let $X$ be the number of heads we observe.

Suppose we wish to find $P(X=3)$

There are ${5 \choose 3}=10$ ways to select the coins that will show up as heads. In particular











The respective probabilities for these are











Summing these, we get

$$P(X=3)\approx 0.286$$

Similarly, for finding $P(X=4)$, there are ${5 \choose 4}=5$ ways to pick the four successes and for finding $P(X=5)$, there are ${5 \choose 5}=1$ ways to pick the five successes.

These computations aren't very fun so perhaps a computer program can be implemented.

  • $\begingroup$ ...The probabilities of heads for each coin is different... $\endgroup$
    – NL628
    Mar 16, 2018 at 3:10
  • 4
    $\begingroup$ @NL628 I misread the problem initially but the "..." at the beginning and end of your comment were needless. $\endgroup$
    – Remy
    Mar 16, 2018 at 8:29

I get approximately 0.3841, according to this Python script which performs the necessary evaluation of cases:

def all_outcomes(num_coins) :
    return [[(m / 2**k) % 2 for k in range(num_coins)[::-1]] for m in range(2**num_coins)]

def at_least_three_heads(coll) :
    # assume list of 0 (tails) and 1 (heads)
    return sum(coll) >= 3

def joint_probability(outcome, coins) :
    ret = 1
    for (is_heads, p) in zip(outcome, coins) :
        ret *= (p if is_heads else 1-p)
    print outcome, ret
    return ret

def compute_probability(coins, success_fn) :
    outcomes = all_outcomes(len(coins))
    return sum([joint_probability(outcome, coins) for outcome in outcomes if success_fn(outcome)])

print compute_probability([0.38, 0.18, 0.71, 0.66, 0.29], at_least_three_heads)

This number makes sense, because most of your coins come up tails most of the time. For comparison, a toss of five fair coins will have three or more heads in half of all cases—there's either three or more heads, or three or more tails. So the probability for five fair coins is 0.5.


Well, the easiest method and most clear one I suppose would be casework. Complementary counting would give the same number of cases.

Find the probability of each of the following:

Try (1, 2, 3) heads and (4, 5) tails, or (1, 2, 4) heads and (3, 5) tails, etc.

Then do (1, 2, 3, 4) heads and (5) tails, etc.

Then do (1, 2, 3, 4, 5) heads.

Then add up the probabilities of everything. To make it easier, I suppose you can write a simple program or code, but this problem is pretty hard as the numbers you gave are pretty random.

  • $\begingroup$ Small question: I'm new to stackexchange mathematics, so would I get flagged if I write code in my solution that could help the OP? $\endgroup$
    – NL628
    Mar 16, 2018 at 3:11

Using a computer might be the best bet for this type of problem. A library that keeps track of combinations, like the Python 3 library trotter (pip3 install trotter, shameless plug :), can keep things sane:

import trotter as tr

ps = [0.38, 0.18, 0.71, 0.66, 0.29]

# Combinations of indices.
combos_sets = [tr.Combinations(i, range(5)) for i in [3, 4, 5]]

prob_sum = 0

def prod(xs):
    p = 1
    for x in xs:
        p *= x
    return p

print("Coins heads: probability of outcome")
for combo_set in combos_sets:
    for combo in combo_set:
        p = prod([ps[i] if i in combo else 1 - ps[i] for i in range(5)])
        print("{}: {}".format(combo, p))
        prob_sum += p

print("Total probability:", prob_sum)

This gives a probability of 0.3841282976...


The distribution of heads is a Poisson binomial distribution. There is no closed form for the probability, but you can hand-code something as proposed by the other answerers. Alternatively, look to standard libraries for this distribution.

Almost exactly the same question was asked yesterday at CrossValidated: Probability mass function with variable probability? There are quite a few threads on the Poisson binomial at CV, e.g., Success of Bernoulli trials with different probabilities.

  • 1
    $\begingroup$ (+1) This is the correct answer. Additionally, there's no need to hand-code anything if you use R -- the R package poibin has tools for the Poisson binomial distribution $\endgroup$
    – duckmayr
    Mar 16, 2018 at 9:43

Encode using probability generating functions (pgfs). Let $X$ be the number of heads after performing the experiment. Then $X=\sum_{i=1}^5 X_i$ where $X_i$ is a Bernoulli trial for the $i$ th coin $H_i$ with probability of success $p_i$. Let $g_X$ be the pgf of $X$. Since the $X_i$ are independent, it follows that $$ \sum_{i=0}^5P(X=i)t^i=g_X(t)=Et^{X}=\prod_i Et^{X_i}=\prod_i [(1-p_i)+p_it] $$ Using the probabilities given in the question $$ \begin{align} g_X(t)&=(0.62+0.38t)(0.82+0.18t)(0.29+0.71t)(0.34+0.66t)(0.71+0.29t)\\ &=0.00929515 t^5 + 0.0888525 t^4 + 0.285981 t^3 + 0.379892 t^2 + 0.200389 t + 0.0355911 \end{align} $$ Thus $$ P(X\geq3)=0.38412865 $$

  • $\begingroup$ Sorry for my ignorance, but what are E and t in this expression? $\endgroup$
    – numbers_q
    Mar 16, 2018 at 19:26

Most convenient is to condition on the state of the "first" coin:

$Pr(3, p_1,p_2,p_3...) = p_1 Pr(2,p_2,p_3...) + (1-p_1) Pr(3, p_2,p_3,...)$

and work recursively, with obvious end conditions $Pr(0,...) = 1$ and $Pr(k,p_1,p_2,...,p_j) = 0 \text{ if } j \lt k$

Here is some Python code:

def pr(ps, n) :
  if n == 0 :
    return 1
  if len(ps) < n :
    return 0
  if len(ps) == n :
    return reduce(lambda x,y : x*y, ps)
  return ps[0] * pr(ps[1:], n-1) + (1-ps[0])*pr(ps[1:],n)

p = [0.38, 0.18, 0.71, 0.66, 0.29]
print pr(p,3)
  • $\begingroup$ I tried exactly this using R, but it runs out of heap space even with $n=5$... $\endgroup$ Mar 16, 2018 at 9:18
  • $\begingroup$ You should include your import statements, $\endgroup$ Mar 16, 2018 at 22:04
  • $\begingroup$ This Python code is self-contained and does not require any imports. $\endgroup$
    – pepster
    Mar 17, 2018 at 22:23

Not the answer you're looking for? Browse other questions tagged .