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]

Roughly how big is N, in practice?

[etherial.livejournal.com]

No idea. That's part of the reason I told him to use a monte carlo method.

[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

[neuromancerzss.livejournal.com]

Haha, you weren't kidding when you said the problem got you two geeking.

I ended up using a recursive function that went through a list of PN's and built toward the target, Y, while keeping track of the probability of the path it was on (and summing them as the recursion returned). At each stage it branched into "heads" and "tails" (multiplying by the correct probabilities) and then continued on. When you hit the target, the probability (for that stage) was 1, or when you hit the end of the list it was 0. I think it's the same thing as your method except it doesn't differentiate between [1,1,1,1] and [1,1,1,-1] when Y = 3 (just bundles the N=4 cases together and says 1).

This code was actually for the case where Y was the sum of the value of all coins that landed heads and the value of each coin flip, VN=PN, but then I realized that with a few parameters you could set the VN and Y to be whatever you want and converted it to then handle my other more standard calculation (it was actually to combine a number of the previous experiments and judge how likely you would be to get Z successful experiments).

I'm not sure that it's mathematically sounds given my actual case (for me Y is a threshold target for an expected value function and VN is really PN*1), but it's an interesting problem.