Hi,
I've been struggling with this for awhile ( http://en.wikipedia.org/wiki/Banach%27s_matchbox_problem)
and I put together this little bit of Python code
def foo(n=3):
from numpy import rand
x = [1,]*n # one box of matches with elements 1
y = [0,]*n # the other box with elements 0
c=[]
while x and y:
if rand()<0.5: # 1/2 probability of picking either box
c.append(x.pop())
else:
c.append(y.pop())
if x: return len(x) # num matches in non-empty box
if y: return len(y)
which simulates the act of picking a match from one box or the other ($1$ or $0$) for this code. The problem is that I can't get the analytical solution for the probability of having $k$ matches left in the remaining box:
$ 1/2^\{(2 n -k )\} Binomial(2 n -k, n) $
For example, if I run 1000 cases for the case where $n=3$ as in
>>> y = [ foo() for i in range ( 1000)]
I get 1 match left in the remaining box 382 times, 2 matches left 375 times, and 3 matches left 243 times. Thus, my estimate of the probability of getting 1 match left in the remaining box is 382/1000, but the analytical solution is
$ 1/2^\{(2 * 3 -1 )\} Binomial(2 *3 -1, 3)= 0.3125$ vs. 0.382
$ 1/2^\{(2 *3 -2 )\} Binomial(2 *3 -2, 3) = 0.25$ vs. 0.375
$ 1/2^\{(2 *3 -3 )\} Binomial(2 *3 -3, 3) = 0.125$ vs. 0.243
which is nowhere near my estimates, and I have tried this with 10e6 trials without getting any closer.
Any help appreciated.