Probability of Outcomes Algorithm

Question

I have a probability problem, which I need to simulate in a reasonable amount of time. In simplified form, I have 30 unfair coins each with a different known probability. I then want to ask things like "what is the probability that exactly 12 will be heads?", or "what is the probability that AT LEAST 5 will be tails?".

I know basic probability theory, so I know I can enumerate all (30 choose x) possibilities, but that's not particularly scalable. The worst case (30 choose 15) has over 150 million combinations. Is there a better way to approach this problem from a computational standpoint?

Any help is greatly appreciated, thanks! :-)

Answer 1

You can use a dynamic programming approach.

For example, to calculate the probability of 12 heads out of 30 coins, let P(n, k) be the probability that there's k heads from the first n coins.

Then P(n, k) = p_n * P(n - 1, k - 1) + (1 - p_n) * P(n - 1, k)

(here p_i is the probability the i'th coin is heads).

You can now use this relation in a dynamic programming algorithm. Have a vector of 13 probabilities (that represent P(n - 1, i) for i in 0..12). Build a new vector of 13 for P(n, i) using the above recurrence relation. Repeat until n = 30. Of course, you start with the vector (1, 0, 0, 0, ...) for n=0 (since with no coins, you're sure to get no heads).

The worst case using this algorithm is O(n^2) rather than exponential.

Answer 2

This is actually an interesting problem. I was inspired to write a blog post about it covering in detail fair vs unfair coin tosses all the way to the OP's situation of having a different probability for each coin. You need a technique called dynamic programming to solve this problem in polynomial time.

General Problem: Given C, a series of n coins p₁ to p_n where p_i represents the probability of the i-th coin coming up heads, what is the probability of k heads coming up from tossing all the coins?

This means solving the following recurrence relation:

P(*n*,k,*C*,i) = p_i x P(*n*-1,k-1,C,*i*+1) + (1-p_i) x P(*n*,k,*C*,i+1)

A Java code snippet that does this is:

private static void runDynamic() {
  long start = System.nanoTime();
  double[] probs = dynamic(0.2, 0.3, 0.4);
  long end = System.nanoTime();
  int total = 0;
  for (int i = 0; i < probs.length; i++) {
    System.out.printf("%d : %,.4f%n", i, probs[i]);
  }
  System.out.printf("%nDynamic ran for %d coinsin %,.3f ms%n%n",
      coins.length, (end - start) / 1000000d);
}

private static double[] dynamic(double... coins) {
  double[][] table = new double[coins.length + 2][];
  for (int i = 0; i < table.length; i++) {
    table[i] = new double[coins.length + 1];
  }
  table[1][coins.length] = 1.0d; // everything else is 0.0
  for (int i = 0; i <= coins.length; i++) {
    for (int j = coins.length - 1; j >= 0; j--) {
      table[i + 1][j] = coins[j] * table[i][j + 1] +
          (1 - coins[j]) * table[i + 1][j + 1];
    }
  }
  double[] ret = new double[coins.length + 1];
  for (int i = 0; i < ret.length; i++) {
    ret[i] = table[i + 1][0];
  }
  return ret;
}

What this is doing is constructing a table that shows the probability that a sequence of coins from p_i to p_n contain k heads.

For a deeper introduction to binomial probability and a discussion on how to apply dynamic programming take a look at Coin Tosses, Binomials and Dynamic Programming.

Answer 3

Just wondering - clearly DP gives an exact solution (nice posts BTW - two thumbs up) but could we approximate the solution with a binomial distribution where p(binomial) = function of (p1,p2,p3....) of the individual coins. Once you implement the DP solution, you might compare the result to one using binomial probabilities and see if a quick approximation is good enough.