Math Geekery

[etherial]

neuromancerzss asked me a quick math question and it got rosinavs and I geeking about probability. His question was basically, "If I have N coins, each with a different probability P_N of coming up heads, what is the probability that Y or more will come up heads?". This problem is trivial if they all have the same P_N, and moderately simple if N or Y is small (we're talking Y<=2 here, kids). But I thought "This should be no problem in matlab!"

1. Create an array of all the possibilities for flipping coins (1=heads and -1=tails). One coin looks like
[ 1;
-1]

- Two coins looks like -

[ 1 1;
1 -1;
-1 1;
-1 -1]

- Three looks like -

[ 1 1 1;
1 1 -1;
1 -1 1;
-1 1 1;
1 -1 -1;
-1 1 -1;
-1 -1 1;
-1 -1 -1]

Etc., etc. There are N_C_N=1 rows of all 1s at the top, N_C_(N-1)=N rows of all but one heads below it, N_C_(N-2) rows of all but two heads below that, etc. Each row is unique and has at least as many 1s as the row below it.

2. You then multiply each entry by (P_N - .5) and add .5 to each entry, resulting in an array with all the probabilities laid out.

3. Multiply the elements of each row to create a column vector giving the probability for that row.

4. Since these are all in order, sum N_C_K from k=Y to N and and take the sum of the top that-many rows, and Presto!, you have P(Y or more are heads). See? Trivial!

And man, did I forget how much I loved matlab. Tonight, I'm going to sit in bed and diagram out the first computer program I've written for myself in a long, long time - the Camelot Chores Generator.

Comments

[agthorr.livejournal.com]

There's also a dynamic programming solution based on the following recursion.

P(i,j,h) = Probability of getting exactly h heads using only coins i through j-1.
P(i,j,h) = 0 if i >= j
P(i,i+1,h) = P_i
P(i,j,h) = sum(P(i,(i+j)//2,k) * P((i+j)//2,j,h-k) for k = 0..h)

Then sum(P(0,n,k) for k=Y..n) to get the solution.

Once memoized, I think it's O(n log n).

[agthorr.livejournal.com]

P(i,j,h) = 0 if i >= j

Correction: P(i,j,h) = 0 if i==j or j-i < h