Distribution of $10$ unique coins

Question

You have $10$ coins each with a different probability ($p$) of landing heads. $p = [.05,.1,.15,.2,.25,...,.45,.5].$

What is the likelihood that at least $5$ coins land heads? Is there any more efficient way of calculating than doing the grunt work calculation of all the probabilities of the different outcomes?

Answer 1

Typically conditional probability is the most efficient computationally. Here is an example python code for your case:

 def pgtk(ps, k) :
  if k == 0 :
    return 1

  if len(ps) < k:
    return 0

  return ps[0] * pgtk(ps[1:], k-1) + (1-ps[0])*pgtk(ps[1:],k)


 >>>>  pgtk([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], 5) 
       0.0974653128125

Answer 2

An efficient way of doing the grunt work would be to use probability generating functions.

Given a discrete random variable $X$ which takes values in $\{0,1,2,\ldots\}$ the probability generating function $G_X$ is defined as, $$ G_X(s) = \mathbb{E}[s^X] = \sum_{k=0}^{\infty}\mathbb{P}(X=k) s^k $$

Probability generating functions have some useful properties. Specifically, the coefficient of the $s^k$ term gives the probability that $X=k$. Therefore, given the generating function $G_X$ of $X$, we can compute, $$ \mathbb{P}(X=k) = \left[ \frac{d^k}{s^k} G_X(s) \right]_{s=0} $$

Given independent random variables $X$ and $Y$ with probability generating functions $G_X$ and $G_Y$, the sum of the variables has generating function, $$G_{X+Y}(s) = G_X(s)G_Y(s)$$

Suppose $X$ is heads with probability $p$. Then $X$ has probability generating function, $$ G_X(s) = (1-p) + ps $$

In your case, you can compute the probability generating function of each of the coins and use this to compute the generating function of their sum. You can then use the formula above to compute $$\mathbb{P}(sum \geq 5) = \sum_{k=5}^{\infty} \mathbb{P}(sum=k) $$

If you want an exact answer, I would suggest using something like Mathematica or Sympy to help with this computation. Otherwise you could do it numerically in python

For instance,

p = {1/20, 1/10, 3/20, 1/5, 1/4, 3/10, 7/20, 2/5, 9/20, 1/2}
k = 5
Total[CoefficientList[Fold[Times, Map[(1 - #) + # s &, p]], s][[k + 1 ;;]]]

This gives, $$ \mathbb{P}(sum \geq 5) = \frac{311889001}{3200000000} \approx 0.0974653$$

If we think about what is going on, we are really just using polynomials in the variable $s$ to keep track of all of the "grunt work" which we would have to do if we wanted to write down all the combinations of ways to get a certain number of heads.

Answer 3

Tyler Chen has the right approach, you can use generating functions:

For a Bernoulli law you have (where $p$ is head probability): $$ G_X(s)=(1−p)+ps $$ The trick of generating functions is to encode all your possible outputs of your discrete random variable {0,1,2,...} in a polynomial where {0,1,2,...} are distinct powers of your polynomial variable s. For the Bernoulli example, you only have two possible outputs 0 with probability 1-p and 1 with probability p, hence the $(1−p)s^0+ps^1$ polynomial. One nice property is that the sum X+Y of two independent variables is the product of the associated generating functions $G_X(s)*G_Y(s)$.

You throw your $10$ coins, hence $X_1+X_2+...+X_{10}$ is translated into: $G_{X_1}(s)*...*G_{X_{10}}(s)$. You will get a polynomial. For each head a coefficient $s$ is involved, hence to get to probability of having at least $5$ heads you simply have to sum all the polynomial coefficients associated to $s^k$ with $k\ge5$.

Here is the Mathematica code detailing each step:

G[p_] := (1 - p) + p*s    
pi = Map[Rationalize, {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5}]
P = Expand[Apply[Times, Map[G, pi]]]
Apply[Plus, CoefficientList[P, s][[6 ;; -1]]] 
(* 6 because Mathematica indices start at 1 and not 0 *)

Your probability is $\frac{311889001}{3200000000}\approx 0.09746534 $

---

Below you will find the complete polynomial: $P=G_{X_1}(s)*...*G_{X_{10}}(s)$

$$ \frac{567 s^{10}}{1600000000}+\frac{55089 s^9}{3200000000}+\frac{523683 s^8}{1600000000}+\frac{5262769 s^7}{1600000000}+\frac{6297587 s^6}{320000000}+\frac{59321001 s^5}{800000000}+\frac{287756563 s^4}{1600000000}+\frac{446890087 s^3}{1600000000}+\frac{213926463 s^2}{800000000}+\frac{91671039 s}{640000000}+\frac{26189163}{800000000} $$