Flipping a set of unfair coins [closed]

Question

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?

Example:

$P(H_1) = .38$

$P(H_2) = .18$

$P(H_3) = .71$

$P(H_4) = .66$

$P(H_5) = .29$

Answer 1

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

$$H_1H_2H_3$$

$$H_1H_2H_4$$

$$H_1H_2H_5$$

$$H_1H_3H_4$$

$$H_1H_3H_5$$

$$H_1H_4H_5$$

$$H_2H_3H_4$$

$$H_2H_3H_5$$

$$H_2H_4H_5$$

$$H_3H_4H_5$$

The respective probabilities for these are

$$P(H_1H_2H_3)=0.38\cdot0.18\cdot0.71\cdot0.34\cdot0.71$$

$$P(H_1H_2H_4)=0.38\cdot0.18\cdot0.29\cdot0.66\cdot0.71$$

$$P(H_1H_2H_5)=0.38\cdot0.18\cdot0.29\cdot0.34\cdot0.29$$

$$P(H_1H_3H_4)=0.38\cdot0.82\cdot0.71\cdot0.66\cdot0.71$$

$$P(H_1H_3H_5)=0.38\cdot0.82\cdot0.71\cdot0.34\cdot0.29$$

$$P(H_1H_4H_5)=0.38\cdot0.82\cdot0.29\cdot0.66\cdot0.29$$

$$P(H_2H_3H_4)=0.62\cdot0.18\cdot0.71\cdot0.66\cdot0.71$$

$$P(H_2H_3H_5)=0.62\cdot0.18\cdot0.71\cdot0.34\cdot0.29$$

$$P(H_2H_4H_5)=0.62\cdot0.18\cdot0.29\cdot0.66\cdot0.29$$

$$P(H_3H_4H_5)=0.62\cdot0.82\cdot0.71\cdot0.66\cdot0.29$$

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.

Answer 2

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.

Answer 3

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.

Answer 4

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...

Answer 5

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.

Answer 6

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 $$

Answer 7

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)