The gambler’s ruin problem is one where a player has a probability p of winning and probability q of losing. For example let’s take a skill game where the player x can beat player y with probability 0.6 by getting closer to a target. The game play begins with player x being allotted 5 points and player y allotted 10 points. After each round a player’s points either decrease by one or increase by one we can determine the probability that player x will annihilate player y. The player that reaches 15 wins and the player that reach zero is annihilated. There is a wide range of application for this type of problem that goes being gambling.
This is actually a fairly simple problem to solve on pencil and paper and to determine an exact probability. Without going into too much detail we can determine the probability of annihilation by . In this example it works out to be
.
But this is a relatively boring approach and coding up an R script makes everything that much better. So here is a simulation of this same problem estimating that same probability plus it provides additional information on the distribution of how many times this game would have to be played.
gen.ruin = function(n, x.cnt, y.cnt, x.p){
x.cnt.c = x.cnt
y.cnt.c = y.cnt
x.rnd = rbinom(n, 1, p=x.p)
x.rnd[x.rnd==0] = -1
y.rnd = x.rnd*-1
x.cum.sum = cumsum(x.rnd)+x.cnt
y.cum.sum = cumsum(y.rnd)+y.cnt
ruin.data = cumsum(x.rnd)+x.cnt
if( any( which(ruin.data>=x.cnt+y.cnt) ) | any( which(ruin.data< =0) ) ){ cut.data = 1+min( which(ruin.data>=x.cnt+y.cnt), which(ruin.data< =0) ) ruin.data[cut.data:length(ruin.data)] = 0 } return(ruin.data) } n.reps = 10000 ruin.sim = replicate(n.reps, gen.ruin(n=1000, x.cnt=5, y.cnt=10, x.p=.6)) ruin.sim[ruin.sim==0] = NA hist( apply(ruin.sim==15 | is.na(ruin.sim), 2, which.max) , nclass=100, col='8', main="Distribution of Number of Turns", xlab="Turn Number") abline(v=mean(apply(ruin.sim==15 | is.na(ruin.sim), 2, which.max)), lwd=3, col='red') abline(v=median(apply(ruin.sim==15 | is.na(ruin.sim), 2, which.max)), lwd=3, col='green') x.annihilation = apply(ruin.sim==15, 2, which.max) ( prob.x.annilate = length(x.annihilation[x.annihilation!=1]) / n.reps ) state.cnt = ruin.sim state.cnt[state.cnt!=15] = 0 state.cnt[state.cnt==15] = 1 mean.state = apply(ruin.sim, 1, mean, na.rm=T) plot(mean.state, xlim=c(0,which.max(mean.state)), ylim=c(0,20), ylab="Points", xlab="Number of Plays", pch=16, cex=.5, col='green') lines(mean.state, col='green') points(15-mean.state, pch=16, cex=.5, col='blue') lines(15-mean.state, col='blue') [/code]
thanks for sharing! a few comments around the code would have made it easier to follow.
same with some comments I added to help me understand:
gen.ruin = function(x.cnt, y.cnt)
{
# n = number of observations of the binomial distribution
# s = number of trials
# p = probability of success on each trial
x.rnd = rbinom(n = 1000, s = 1, p = 0.6)
x.rnd[x.rnd==0] = -1
y.rnd = x.rnd*-1
x.cumsum = cumsum(x.rnd)+x.cnt
y.cumsum = cumsum(y.rnd)+y.cnt
ruin.data = x.cumsum # same as x.cumsum
if( any( which( ruin.data >= x.cnt+y.cnt ) ) | any( which( ruin.data = x.cnt+y.cnt ), which( ruin.data <= 0) )
ruin.data[cut.data:length(ruin.data)] = 0
}
return(ruin.data)
}
# set number of replications
n.reps = 10000
# set initial counters
x.cnt.init = 5
y.cnt.init = 10
# replicate n.reps simulations of ruin.sim
ruin.sim = replicate(n.reps, gen.ruin(x.cnt = x.cnt.init, y.cnt = y.cnt.init))
ruin.sim[ruin.sim==0] = NA
# plot the histogram
# pdf("histogram.pdf")
hist( apply( ruin.sim==15 | is.na(ruin.sim), 2, which.max )
, nclass = 100
, col = '8'
, main = "Distribution of Number of Turns"
, xlab = "Turn Number"
)
# add a vertical line at the mean value
abline( v = mean(apply(ruin.sim==15 | is.na(ruin.sim), 2, which.max) )
, lwd = 3
, col = 'red'
)
# add a vertical line at the median value
abline( v = median(apply(ruin.sim==15 | is.na(ruin.sim), 2, which.max) )
, lwd = 3
, col = 'green'
)
# dev.off()
# _apply_ applies a function to margins of an array or matrix
# second arg: for a matrix 1 indicates rows, 2 indicates columns
# set maximum number of wins
win.max = 15
# conditions for annihilation: return the maximizer
# mean over columns of mean.state
x.annihilation = apply(ruin.sim==win.max, 2, which.max)
# probability of annihilation for n.reps replications
prob.x.annilate = length(x.annihilation[x.annihilation!=1]) / n.reps
print(prob.x.annilate)
# copy ruin.sim as state.cnt
state.cnt = ruin.sim
# look inside state.cnt: cases not equal to win.max
# print(state.cnt!=win.max)
state.cnt[state.cnt!=win.max] = 0
# look inside state.cnt: cases equaul to win.max
# print(state.cnt==win.max)
state.cnt[state.cnt==win.max] = 1
# mean over rows of mean.state
mean.state = apply( ruin.sim, 1, mean, na.rm = TRUE )
# plot of connected points for score by number of plays
# pdf("score.pdf")
plot( mean.state
, xlim = c(0,which.max(mean.state))
, ylim = c(0,20)
, ylab = "Score"
, xlab = "Number of Plays"
, pch = 16
, cex = .5
, col = 'green'
)
lines(mean.state, col = 'green')
points(win.max-mean.state, pch = 16, cex = .5, col = 'blue')
lines(win.max-mean.state, col = 'blue')
# dev.off()
Simulating the Gambler’s Ruin in R http://t.co/jMKqK7a02f
RT @ds_ldn: Simulating the Gambler’s Ruin in R http://t.co/jMKqK7a02f