The Gambler's Ruin problem is a special case of a "random walk". In this notebook we use random numbers to do Monte Carlo simulations and compare the results with some mathematical expressions. The basic idea is that a gambler has a finite amount of money and continues betting until the money runs out. In the version used here there are two gamblers and the game ends when one gambler is ruined (and the other has won all the money).
Suppose two gamblers have $n_1$ and $n_2$ pennies to start. Each round of betting is won by player 1 with probability $p$ or by player 2 with probability $q=1-p$. The loser transfers a penny to the winner. The game ends when one player has no pennies left.
It can be shown that player 1 wins with probability $$ P_1 = \frac{1-(q/p)^{n_1}} {1-(q/p)^{n_1+n_2}} $$
and Player 2 wins with probability $$ P_2 = 1- P_1 = \frac{(q/p)^{n_1}-(q/p)^{n_1+n_2}} {1-(q/p)^{n_1+n_2}} $$ See http://en.wikipedia.org/wiki/Gambler's_ruin#Unfair_coin_flipping.
In [ ]:
%pylab inline
In [ ]:
p = 0.6
q = 1-p
n1 = 3
n2 = 3
P1 = (1-(q/p)**n1) / (1.-(q/p)**(n1+n2))
P2 = ((q/p)**n1 - (q/p)**(n1+n2)) / (1.-(q/p)**(n1+n2))
print P1,P2, P1+P2
In [ ]:
from numpy import random
random.seed(12345)
A function to take a single random walk:
In [ ]:
max_steps = 500 # maximum number of steps for each random walk
def walk(n1,n2,p,verbose=True):
n1path = [n1] # keep track of n1 each step of the walk
for nstep in range(1,max_steps+1):
r = random.random(1)
if r <= p:
n1 = n1+1
n2 = n2-1
else:
n1 = n1-1
n2 = n2+1
n1path.append(n1)
if verbose:
print "In step %i, r = %7.4f and n1 = %i, n2 = %i" % (nstep,r,n1,n2)
if min(n1,n2) == 0:
break
if verbose:
print "Stopped after %i steps with n1 = %i, n2 = %i" % (nstep,n1,n2)
if n1 == 0:
winner = 2
elif n2 == 0:
winner = 1
else:
winner = 0
return nstep, winner, n1path
Test this out by taking a single walk and plotting the results
In [ ]:
n1 = 4
n2 = 4
nstep, winner, n1path = walk(n1, n2,.6)
plot(n1path,'o-')
xlabel('Step in walk')
ylabel('Number of pennies player 1 has')
yt = yticks(range(n1+n2+1))
Try several walks and plot the results together:
In [ ]:
random.seed(1234)
for k in range(7):
nstep, winner, n1path = walk(5,5,.6,verbose=False)
plot(n1path,'o-')
xlabel('Step in walk')
ylabel('n1')
In [ ]:
p = 0.6
n1 = 3
n2 = 3
wins = [0,0,0]
kwalks = 1000
for k in range(kwalks):
nsteps, winner, n1path = walk(n1, n2, p, verbose=False)
wins[winner] += 1
frac1wins = wins[1] / float(kwalks)
frac2wins = wins[2] / float(kwalks)
if wins[0] > 0:
print "Warning: %i walks out of %i did not result in a win by either player" \
% (wins[0],kwalks)
print "Player 1 won %s fraction of the time, Player 2 won %s fraction of the time" \
% (frac1wins,frac2wins)
q = 1-p
P1 = (1-(q/p)**n1) / (1.-(q/p)**(n1+n2))
P2 = 1 - P1
print "True probabilities are P1 = %7.4f, P2 = %7.4f" % (P1,P2)
There's also a formula for the expected value of the number of steps in the walk, which for $p \neq 1/2$ is given by $$ \frac{n_1}{q-p} - \frac{n_1+n_2}{q-p} \frac{1-(q/p)^{n_1}}{1-(q/p)^{n_1+n_2}} $$ See for example http://www.dam.brown.edu/people/huiwang/classes/APMA1200/RWalk_07.pdf.
We can test this out:
In [ ]:
p = 0.8
q = 1 - p
n1 = 50
n2 = 50
wins = [0,0,0]
kwalks = 1000
nsteps_total = 0
for k in range(kwalks):
nsteps, winner, n1path = walk(n1, n2, p, verbose=False)
nsteps_total += nsteps
wins[winner] += 1
if wins[0] > 0:
print "Warning: %i walks out of %i did not result in a win by either player" \
% (wins[0],kwalks)
print "The average path length is %7.4f" % (nsteps_total/float(kwalks))
mean_length = (n1/(q-p)) - ((n1+n2)/(q-p))*(1-(q/p)**n1)/(1-(q/p)**(n1+n2))
print "True mean path length is %7.4f" % mean_length