Simulating the Gambler’s Ruin

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.

Number of Time Until Annihilation

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 \frac{1-\left(\frac{q}{p}\right)^i}{1-\left(\frac{q}{p}\right)^N}. In this example it works out to be \frac{1-\left(\frac{.4}{.6}\right)^5}{1-\left(\frac{.4}{.6}\right)^{10}} \approx 0.8703.

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

Gambler's Ruin

m4s0n501
Leave a comment

4 Comments

  1. annoporci

     /  April 15, 2013

    thanks for sharing! a few comments around the code would have made it easier to follow.

    Reply
  2. annoporci

     /  April 15, 2013

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

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *


− one = 2

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

%d bloggers like this: