# coding: utf-8 # Math 640 Project - Gambler's Ruin # ================================= # # Jeff Ames, Kafung Mok # # 2014-05-06 # Abstract # -------- # # We investigate the Gambler's Ruin problem, extending it to two # dimensions, and look into the following questions: # # * On an $N \times N$ grid, what is the expected duration of a game? # # * What is the probability of winning? # # * How does a grid size of $M \times N$ affect these results? # # * How do different definitions of winning (e.g., reaching edges # vs. reaching corners) affect these results? ### Initialization # In[42]: # Import required libraries import itertools import matplotlib as ml import matplotlib.pyplot as plt import numpy as np import sympy as sp # In[43]: # Set plots to render inline in the notebook get_ipython().magic('matplotlib inline') ### Expected duration # In[44]: # ExpD2 calculates the expected duration of a random walk in the # plane. It outputs the two-dimensional vector [f(1,1),f(1,2),..., # f(N-1,N-1)] of expected durations with prob. 1/4 of moving in each # of four directionsfff, where f(i,j) is the expected duration of the # walk starting at i,j: # # f(i,j) = p*f[i-1,j] + p*f[i+1,j] + p*f[i,j-1] + p*f[i,j+1] + 1 # # f(0,k) = 0 # f(N,k) = 0 # f(k,0) = 0 # f(k,M) = 0 # def ExpD2(M, N): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list p = 1/4 boundaries = [(f[0][0], 0), (f[M][0], 0), (f[0][N], 0), (f[M][N], 0)] for k in range(1, M): boundaries += [(f[k][0], 0), (f[k][N], 0)] for k in range(1, N): boundaries += [(f[0][k], 0), (f[M][k], 0)] eq = [] for j in range(1, N): for i in range(1, M): eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[45]: exp = ExpD2(12, 12) exp.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types exp = [[float(exp[i][j]) for j in range(len(exp[0]))] for i in range(len(exp))] # Plot expected duration in 2D H = np.array(exp) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() ### Probability of winning #### In one dimension # In[46]: # ProbW1D inputs a positive integer N, a number p between 0 and 1, and # outputs a list of length N-1 whose i-th entry tells you the # probability of exiting the casino a winner, if you start out with i # dollars, and must leave as soon as you are either broke or have N # dollars. def ProbW1D(N, p): f = [sp.Symbol("f[{0}]".format(i)) for i in range(N+1)] boundaries = [(f[0], 0), (f[N], 1)] eq = [(p*f[i+1] + (1-p)*f[i-1] - f[i]).subs(boundaries) for i in range(1, N)] soln = sp.solve(eq, f) return [var.subs(soln).subs(boundaries) for var in f] # In[56]: # Parameters: n is the amount you must have to win the game. probs is # an array of probabilities. The second item in the probs array is a # line style for plotting the results afterwards. n = 50 probs = [(0.45, 'g--'), (0.475, 'c:'), (0.5, 'b-'), (0.525, 'm.'), (0.55, 'r-.')] p = [ProbW1D(n, prob[0]) for prob in probs] # Convert from sympy types to native python types p = [[float(pr) for pr in p[i]] for i in range(len(p))] # Plot probabilities of winning args = [(range(n+1), p[i], probs[i][1]) for i in range(len(p))] args_list = itertools.chain(*args) # flatten to a one-dimensional list fig = plt.figure(figsize=(10, 6)) plt.plot(*args_list) plt.legend(["{:.2f}".format(p[0]) for p in probs], loc='center left', bbox_to_anchor=(1, 0.5)); #### In two dimensions # In[48]: # ProbW2D inputs positive integers M, N, a number p between 0 and 1, and # outputs a 2D list whose i,j-th entry tells you the # probability of exiting the casino a winner, if you start out with # i,j dollars, and must leave as soon as you are either broke or have # M or N dollars. def ProbW2D(M, N, p): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list boundaries = [(f[0][0], 0), (f[M][0], 1), (f[0][N], 1), (f[M][N], 1)] for k in range(1, M): boundaries += [(f[k][0], 0), (f[k][N], 1)] for k in range(1, N): boundaries += [(f[0][k], 0), (f[M][k], 1)] eq = [] for j in range(1, N): for i in range(1, M): eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[49]: prob = ProbW2D(12, 12, 1/4) prob.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types prob = [[float(prob[i][j]) for j in range(len(prob[0]))] for i in range(len(prob))] # Plot expected duration in 2D H = np.array(prob) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() ### Non-square grids # Suppose that the size of the grid is $M \times N$, where $M \ne N$. Then the expected duration is the following: # In[50]: exp = ExpD2(10, 20) exp.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types exp = [[float(exp[i][j]) for j in range(len(exp[0]))] for i in range(len(exp))] # Plot expected duration in 2D H = np.array(exp) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() # And the probability of winning is the following: # In[51]: prob = ProbW2D(10, 20, 1/4) prob.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types prob = [[float(prob[i][j]) for j in range(len(prob[0]))] for i in range(len(prob))] # Plot expected duration in 2D H = np.array(prob) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() ### Corners vs. edges # In this case, instead of winning when you reach the top or right edge # and losing when you reach the bottom or left edge, you must reach one # of the corners. We denote the bottom left corner as losing (both # currencies are at 0) and the other three currencies as winning. # # If you are at an interior square, you have equal probability of moving # in each of the four directions ($p = 1/4$). However, if you reach an # edge, we redistribute the probabilities since you can no longer go in # one direction, so the new probability $p_{\text{edge}}$ for each of the three directions # is $1/3$. # The expected duration is the following: # In[52]: def ExpD2corner(M, N): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list p = 1/4 boundaries = [(f[0][0], 0), (f[M][0], 0), (f[0][N], 0), (f[M][N], 0)] eq = [] for j in range(1, N): for i in range(1, M): eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] for j in range(1, N): p = 1/3 i = 0 eq += [(p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] i = M eq += [(p*f[i-1][j] + p*f[i][j-1] + p*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] for i in range(1, M): p = 1/3 j = 0 eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] j = N eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] + 1 - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[53]: exp = ExpD2corner(12, 12) exp.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types exp = [[float(exp[i][j]) for j in range(len(exp[0]))] for i in range(len(exp))] # Plot expected duration in 2D H = np.array(exp) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() # The probability of winning is the following: # In[54]: # ProbW2Dcorner inputs positive integers M, N, a number p between 0 and 1, and # outputs a list of length N-1 whose i-th entry tells you the # probability of exiting the casino a winner, if you start out with # i,j dollars, and must leave as soon as you are either broke or have # M or N dollars. def ProbW2Dcorner(M, N, p): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list boundaries = [(f[0][0], 0), (f[M][0], 1), (f[0][N], 1), (f[M][N], 1)] eq = [] for j in range(1, N): for i in range(1, M): eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] - f[i][j]).subs(boundaries)] for j in range(1, N): p = 1/3 i = 0 eq += [(p*f[i+1][j] + p*f[i][j-1] + p*f[i][j+1] - f[i][j]).subs(boundaries)] i = M eq += [(p*f[i-1][j] + p*f[i][j-1] + p*f[i][j+1] - f[i][j]).subs(boundaries)] for i in range(1, M): p = 1/3 j = 0 eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j+1] - f[i][j]).subs(boundaries)] j = N eq += [(p*f[i-1][j] + p*f[i+1][j] + p*f[i][j-1] - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[55]: prob = ProbW2Dcorner(12, 12, 1/4) prob.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types prob = [[float(prob[i][j]) for j in range(len(prob[0]))] for i in range(len(prob))] # Plot expected duration in 2D H = np.array(prob) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() ### Biased coins # Suppose instead that there are independent probabilities for each currency, with probability $p_i$ of winning a dollar and probability $1-p_i$ of losing a dollar. # Then the expected duration is of the following form: # In[68]: def ExpD2bias(M, N, p): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list boundaries = [(f[0][0], 0), (f[M][0], 0), (f[0][N], 0), (f[M][N], 0)] for k in range(1, M): boundaries += [(f[k][0], 0), (f[k][N], 0)] for k in range(1, N): boundaries += [(f[0][k], 0), (f[M][k], 0)] # Normalize so all probabilities sum to 1 norm = 2 # p[0] + (1-p[0]) + p[1] + (1-p[1]) p_left = (1-p[0]) / norm p_right = p[0] / norm p_up = p[1] / norm p_down = (1-p[1]) / norm eq = [] for j in range(1, N): for i in range(1, M): eq += [(p_left*f[i-1][j] + p_right*f[i+1][j] + p_down*f[i][j-1] + p_up*f[i][j+1] + 1 - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[80]: exp = ExpD2bias(12, 12, [2/3, 1/4]) exp.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types exp = [[float(exp[i][j]) for j in range(len(exp[0]))] for i in range(len(exp))] # Plot expected duration in 2D H = np.array(exp) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show() # And the probability of winning is of the following form: # In[77]: def ProbW2Dbias(M, N, p): f = [[sp.Symbol("f[{}][{}]".format(i,j)) for j in range(N+1)] for i in range(M+1)] f_list = itertools.chain(*f) # flatten to a one-dimensional list boundaries = [(f[0][0], 0), (f[M][0], 1), (f[0][N], 1), (f[M][N], 1)] for k in range(1, M): boundaries += [(f[k][0], 0), (f[k][N], 1)] for k in range(1, N): boundaries += [(f[0][k], 0), (f[M][k], 1)] norm = 2 # p[0] + (1-p[0]) + p[1] + (1-p[1]) p_left = (1-p[0]) / norm p_right = p[0] / norm p_up = p[1] / norm p_down = (1-p[1]) / norm eq = [] for j in range(1, N): for i in range(1, M): eq += [(p_left*f[i-1][j] + p_right*f[i+1][j] + p_down*f[i][j-1] + p_up*f[i][j+1] - f[i][j]).subs(boundaries)] soln = sp.solve(eq, f_list) return [[f[i][j].subs(soln).subs(boundaries) for j in range(N+1)] for i in range(M+1)] # In[79]: prob = ProbW2Dbias(12, 12, [2/3, 1/4]) prob.reverse() # so that (0,0) is in the bottom left # Convert from sympy types to native python types prob = [[float(prob[i][j]) for j in range(len(prob[0]))] for i in range(len(prob))] # Plot expected duration in 2D H = np.array(prob) fig = plt.figure(figsize=(10, 6)) plt.imshow(H) plt.colorbar(orientation='vertical') plt.show()