duration.py

#! /usr/bin/env python3
## NAME: duration.py
## USAGE: From shell prompt: python3 duration.py
##        or within interactive python3 environment, import filename
## REQUIRED ARGUMENTS: none
## OPTIONS: none
## DESCRIPTION: Experiment to measure duration 
##              of random walk of n steps done k times each,
##              starting from each value from ruin to victory,
##              continuing until reaching either ruin or victory
##              equivalent to absorbing random walk, or gambler's ruin
## DIAGNOSTICS: none
## CONFIGURATION AND ENVIRONMENT: requires wheel matplotlib numpy scipy
## DEPENDENCIES: requires wheel matplotlib numpy scipy
## INCOMPATIBILITIES: none known
## PROVENANCE:  Steven R. Dunbar Wed Aug  8, 2012  5:30 AM
## BUGS AND LIMITATIONS: None known
## FEATURES AND POTENTIAL IMPROVEMENTS:  
## AUTHOR:  Steve Dunbar
## VERSION: Version 1.0 as of Wed Aug  8, 2012  5:33 AM
##          Version 1.1 as of Tue Mar  1, 2016  7:31 AM, upgrade to python3
##          Version 1.2 as of Fri Mar 11, 2016  7:40 AM  add matplotlib graphics
##          Version 1.3 as of Sun May 11, 2020  6:34 AM  transition to numpy
## KEYWORDS: Coin flips, Gambler's Ruin, Absorbing Random Walk

import numpy as np
import matplotlib.pyplot as plt

p = 0.495
n = 1000
k = 500
victory = 100;
ruin = 0;
interval = victory - ruin + 1

winLose = 2*( np.random.random((n,k,interval)) <= p ) - 1
totals = np.cumsum(winLose, axis = 0)

start = np.multiply.outer( np.ones((n+1,k), dtype=int), np.arange(ruin, victory+1, dtype=int))
paths = np.zeros((n+1,k,interval), dtype=int)
paths[ 1:n+1, :,:] = totals
paths = paths + start

def match(a,b,nomatch=None):
    return  b.index(a) if a in b else nomatch
  
# arguments: a is a scalar, b is a python list, value of nomatch is scalar
# returns the position of first match of its first argument in its second argument
# but if a is not there, returns the value nomatch
# modeled on the R function "match", but with less generality

hitVictory = np.apply_along_axis(lambda x:( match(victory,x.tolist(),nomatch=n+2)), 0, paths)
hitRuin = np.apply_along_axis(lambda x:match(ruin,x.tolist(),nomatch=n+2), 0, paths)

# If no ruin or victory on a walk, nomatch=n+2 sets the hitting
# time to be two more than the number of steps, one more than
# the column length.

durationUntilRuinOrVictory = np.minimum(hitVictory, hitRuin)
durationMasked = np.ma.masked_greater(durationUntilRuinOrVictory, n)

startValues = np.arange(ruin, victory+1, dtype=int)
meanDuration = np.mean(durationUntilRuinOrVictory, axis = 0)
durationFunction = np.polyfit( startValues, meanDuration, 2)

print(( "Duration function is: ", durationFunction[2], "+", durationFunction[1], "x+", durationFunction[0], "x^2")) 

# should return coefficients to (x-ruin)*(victory - x), descending degree order

plt.scatter(startValues, meanDuration)
X = np.linspace(ruin, victory, 256, endpoint=True)
Y = durationFunction[2] + durationFunction[1]*X + durationFunction[0]*X**2
plt.plot(X,Y)
plt.show()