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

## annoporci

/ April 15, 2013thanks for sharing! a few comments around the code would have made it easier to follow.

## annoporci

/ April 15, 2013same 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()