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 
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