# Question about Banach's matchbox problem.

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.

-
The theoretical probabilities should add to 1 so you already have a problem. This is a programming question. –  Aaron Meyerowitz Dec 22 '10 at 22:55
Since the theoretical probabilities don't add up to 1, I would think there's something wrong with the analytical formula, and not with the program. –  Peter Shor Dec 22 '10 at 22:56
Don't forget that the stop time is when one box fails to get to -1. The 2nd box could be empty k=0 when the first is discovered to be –  Aaron Meyerowitz Dec 22 '10 at 23:05

Your mistake has to do with the definition of the problem. The man does not stop taking matches as soon as one of the boxes is empty (as in your code). He stops taking matches after he picks a box and finds it empty. This means that when Box A is empty, he doesn't immediately stop -- he continues until he picks Box A again (or until he empties Box B as well).

EDIT: In fact, it's trivial to fix your code -- you don't even have to rewrite your function! Just add the following:

def bar (n=3):
return foo(n+1) - 1


This just exploits the fact that picking the empty matchbox in the $n$-match problem is equivalent to emptying the box in the $(n+1)$-match problem.

-
Thanks! Big help! –  reckoner Dec 23 '10 at 0:34

What you're doing wrong is stopping when the first matchbox runs out. You have to allow $k=0$ to get the probabilities sum to 1, so you have to stop the simulation the first time an empty matchbox is pulled from a pocket (at which point the other pocket might have an empty matchbox in it as well).

-