R simulation of chutes and ladders

chutes_and_ladders.R

# function to simulate chutes and ladders game
chutesNladders <-
function(nplayers=1)
{
  transitions <- rbind(
    c(1, 38),
    c(4, 14),
    c(9, 31),
    c(16, 6),
    c(21, 42),
    c(28, 84),
    c(36, 44),
    c(48, 26),
    c(49, 11),
    c(51, 67),
    c(56, 53),
    c(62, 19),
    c(64, 60),
    c(71, 91),
    c(80, 100),
    c(87, 24),
    c(93, 73),
    c(95, 75),
    c(98, 78))

  transmat <- 1:106
  names(transmat) <- as.character(1:106)
  transmat[transitions[,1]] <- transitions[,2]

  lastpos <- 0
  curpos <- history <- NULL
  while(all(curpos < 100)) {
    curpos <- lastpos + sample(1:6, nplayers, repl=TRUE)
    curpos <- transmat[curpos]
    if(any(curpos > 100)) curpos[curpos > 100] <- lastpos[curpos > 100]

    lastpos <- curpos
    history <- rbind(history, curpos)
  }

  history
}


# for output from chutesNladders, which of the players won?
whowon <-
function(chutesOutput)
  min(which(chutesOutput[nrow(chutesOutput),] == 100))

# for output from chutesNladders, how many total spins were there?
nspins <-
function(chutesOutput)
  ncol(chutesOutput)*(nrow(chutesOutput)-1) +
    whowon(chutesOutput)

chutes_and_ladders_numerical.R

# Regarding Chutes & Ladders:
# numerical calculations (rather than simulations)

# calcTM: calculate transition matrix
#   (101 x 101), first is starting position; rest are board locations 1-100
calcTM <-
function()
{
  transitions <- rbind(
    c(1, 38),
    c(4, 14),
    c(9, 31),
    c(16, 6),
    c(21, 42),
    c(28, 84),
    c(36, 44),
    c(48, 26),
    c(49, 11),
    c(51, 67),
    c(56, 53),
    c(62, 19),
    c(64, 60),
    c(71, 91),
    c(80, 100),
    c(87, 24),
    c(93, 73),
    c(95, 75),
    c(98, 78))

  # simple transition matrix, not accounting for chutes and ladders
  tm <- matrix(0,ncol=107, nrow=107)
  dimnames(tm) <- list(as.character(0:106), as.character(0:106))
  tm[col(tm) > row(tm) & col(tm) < row(tm) + 7] <- 1/6

  # account for the chutes and ladders
  for(i in 1:nrow(transitions)) {
    tm[,transitions[i,2]+1] <- tm[,transitions[i,2]+1] + tm[,transitions[i,1]+1]
    tm[,transitions[i,1]+1] <- 0
  }

  # can't go past 100
  # (not efficient code, but who cares)
  for(i in 1:nrow(tm)) {
    tm[i,i] <- tm[i,i] + sum(tm[i,102:107])
    tm[i,102:107] <- 0
  }

  tm[1:101,1:101]
}

# calculate transition matrix
tm <- calcTM()
# start vector
start <- rbind(c(1, rep(0, 100)))
# vector to pull out end probability
end <- cbind(c(rep(0, 100), 1))

# calculation full distribution at each step
maxsteps <- 400
probs <- matrix(nrow=maxsteps, ncol=ncol(tm))
rownames(probs) <- as.character(1:maxsteps)
probs[1,] <- start %*% tm
for(i in 2:maxsteps)
  probs[i,] <- probs[i-1,] %*% tm

# Pr(game complete by round k, 1 player)
cdf <- probs %*% end

# pr(game complete by round k, n players)
cdf <- cbind(cdf)
for(n in 2:4)
  cdf <- cbind(cdf, 1-(1-cdf[,1])^n)
colnames(cdf) <- as.character(1:4)

# pr(complete game at round k, n players), by differences
stopprob <- apply(cdf, 2, function(a) diff(c(0, a)))

library(broman)
library(RColorBrewer)
col <- brewer.pal(3, "Dark2")
br <- seq(0.5, nrow(stopprob)+0.5, by=1)
png("chutes_and_ladders_rounds.png", height=900, width=1800, res=288, pointsize=12)
par(mar=c(3.1, 4.1, 1.1, 1.1))
grayplot(histlines(br, stopprob[,1]), type="l", lwd=2, xlim=c(0, 150), xaxs="i", yaxs="i",
         ylim=c(0, max(stopprob)*1.05), ylab="", xlab="No. rounds", yat=seq(0, 0.05, by=0.01),
         xat=sort(c(7, seq(0, 150, by=25))), vlines=seq(0, 150, by=25), mgp=c(1.6, 0.2, 0),
         hlines=seq(0, 0.05, by=0.01))
title(ylab="Probability", mgp=c(2.4, 0, 0))
for(i in 2:4)
  lines(histlines(br, stopprob[,i]), col=col[i-1], lwd=2)
for(i in 1:4) {
  dens <- histlines(br, stopprob[,i])
  ymx <- max(dens[,2])
  xmx <- max(dens[dens[,2]==max(dens[,2]),1])
  text(xmx+c(10, 9, 6, 6)[i], ymx*0.95, paste0(i, " player", c("", "s", "s","s")[i]),
       adj=c(0, 0.5), col=c("black", col)[i])
}
dev.off()

# pr(player j wins at round k)
#   cdf = Pr(end by round k for a single player)
prob.player.wins <-
function(cdf, n.players=4)
{
  if(is.na(match(n.players, 2:4)))
    stop("n.players should be in {2,3,4}")

  stopprob <- diff(c(0, cdf))

  result <- cbind(stopprob * (stopprob + 1-cdf)^(n.players-1))

  if(n.players == 2)
    return(cbind(result, (1-cdf)*stopprob))

  result <- cbind(result, (1-cdf)*stopprob*(stopprob + 1-cdf)^(n.players-2))

  if(n.players == 3)
    return(cbind(result, (1-cdf)^2*stopprob))

  cbind(result, (1-cdf)^2*stopprob*(stopprob + 1-cdf),
        (1-cdf)^3*stopprob)
}

p2 <- prob.player.wins(cdf[,1], 2)
p3 <- prob.player.wins(cdf[,1], 3)
p4 <- prob.player.wins(cdf[,1], 4)


png("advantage_to_first_player.png", height=900, width=1800, res=288, pointsize=12)
par(mar=c(3.1, 4.1, 1.1, 1.1))
grayplot(7:101, (p2[,1]/rowSums(p2))[7:101]*2, type="l", lwd=2, xlim=c(6.5, 100.5), xaxs="i", yaxs="i",
         ylim=c(0.998, 1.068), ylab="", xlab="No. rounds to complete game", yat=seq(1, 1.06, by=0.01),
         xat=c(7, seq(25, 100, by=25)), vlines=c(7,seq(25, 100, by=25)), mgp=c(1.6, 0.2, 0),
         hlines=seq(1, 1.06, by=0.01), col=col[1])
title(ylab="Relative advantage for player 1", mgp=c(2.4, 0, 0))
lines(7:101, (p3[,1]/rowSums(p3))[7:101]*3, col=col[2], lwd=2)
lines(7:101, (p4[,1]/rowSums(p4))[7:101]*4, col=col[3], lwd=2)
mx <- c(max(p2[,1]/rowSums(p2)*2, na.rm=TRUE),
        max(p3[,1]/rowSums(p3)*3, na.rm=TRUE),
        max(p4[,1]/rowSums(p4)*4, na.rm=TRUE))
text(rep(99, 3), mx+0.0015, paste(2:4, "players"), adj=c(1, 0), col=col)
dev.off()


# pr(finish after k spins)
#   cdf = Pr(end by round k for a single player)
probBySpins <-
function(cdf, n.players=4)
{
  if(is.na(match(n.players, 2:4)))
    stop("n.players should be in {2,3,4}")

  # prob player k wins at round n
  p <- prob.player.wins(cdf, n.players)
  p <- p/rowSums(p)

  cdf <- 1 - (1-cdf)^n.players
  stopprob <- diff(c(0, cdf))

  p <- as.numeric(t(p * stopprob))
  p[is.na(p)] <- 0
  p
}

stopprobBySpins <-
  cbind(stopprob[,1],
        probBySpins(cdf[,1], 2)[1:nrow(stopprob)],
        probBySpins(cdf[,1], 3)[1:nrow(stopprob)],
        probBySpins(cdf[,1], 4)[1:nrow(stopprob)])

breaks <- seq(0.5, nrow(stopprobBySpins)+0.5, by=1)
dens <- matrix(ncol=ncol(stopprobBySpins)+1, nrow=length(breaks)*2)
for(i in 1:ncol(stopprobBySpins)) {
  tmp <- histlines(breaks, stopprobBySpins[,i], use="density")
  if(i==1) dens[,1] <- tmp[,1]
  dens[,i+1] <- tmp[,2]
}

png("chutes_and_ladders_spins_exact.png", height=900, width=1800, res=288, pointsize=12)
par(mar=c(3.1, 4.1, 1.1, 1.1))
col <- brewer.pal(3, "Dark2")
grayplot(dens[,1], dens[,2], xlab="No. spins", ylab="", col="black", lwd=2, type="l",
         yat = seq(0, 0.025, by=0.005), hlines=seq(0, 0.025, by=0.005),
         xlim=c(0, 200), xaxs="i", yaxs="i", ylim=c(0, max(dens[,2])*1.03),
         xat=seq(0, 200, by=50), vlines=seq(0, 200, by=50), mgp=c(1.8, 0.3, 0))
title(ylab="Probability", mgp=c(2.7, 0, 0))
for(i in 3:ncol(dens))
  lines(dens[,1], dens[,i], lwd=2, col=col[i-2])
ymx <- apply(dens[,-1], 2, max)
xmx <- apply(dens[,-1], 2, function(a,b) max(b[a==max(a)]), dens[,1])
text(xmx+c(8, 12, 10, 20), ymx*0.95, paste0(1:4, " player", c("", "s", "s", "s")),
     adj=c(0, 0.5), col=c("black",col))
dev.off()

sim_chutesNladders.R

source("chutes_and_ladders.R")

n.players <- 2:4
n.sim <- 400000

# simulate different numbers of players and summarize results
simSum <- function(junk, n.players) 
{
  result <- rep(NA, length(n.players)*2)
  for(i in seq(along=n.players)) {
    out <- chutesNladders(n.players[i])
    result[i] <- nspins(out)
    result[i+length(n.players)] <- whowon(out)
  }
  result
}  
  
# simulation in parallel
library(parallel)
RNGkind("L'Ecuyer-CMRG")
mc.cores <- detectCores()
result <- mclapply(1:n.sim, simSum, n.players=n.players, mc.cores=mc.cores)
result <- matrix(unlist(result), nrow=n.sim, byrow=TRUE)
result.nspins <- result[,1:3]
result.whowon <- result[,4:6]
colnames(result.nspins) <- colnames(result.whowon) <- n.players

# save result
save(result.nspins, result.whowon, file="chutes_results.RData")

# average no. spins
colMeans(result.nspins)

# 90th percentile
apply(result.nspins, 2, quantile, 0.9)

# chance that the first player wins
colMeans(result.whowon == 1)

# histograms of the no. spins by no. players
library(RColorBrewer)
library(broman) # get it from github: http://github.com/kbroman/broman

breaks <- seq(-0.5, max(result.nspins)+0.5, by=1)
dens <- matrix(ncol=ncol(result.nspins)+1, nrow=length(breaks)*2)
for(i in 1:ncol(result.nspins)) {
  tmp <- histlines(result.nspins[,i], breaks=breaks, use="density")
  if(i==1) dens[,1] <- tmp[,1]
  dens[,i+1] <- tmp[,2]
}

png("chutes_and_ladders_spins.png", height=900, width=1800, res=288, pointsize=12)
par(mar=c(3.1, 4.1, 1.1, 1.1))
col <- brewer.pal(length(n.players), "Dark2")
grayplot(dens[,1], dens[,2], xlab="No. spins", ylab="", col=col[1], lwd=2, type="l",
         yat = seq(0, 0.02, by=0.005), hlines=seq(0, 0.02, by=0.005),
     xlim=c(0, 200), xaxs="i", yaxs="i", ylim=c(0, max(dens[,2])*1.05),
         xat=seq(0, 200, by=50), vlines=seq(0, 200, by=50), mgp=c(1.8, 0.3, 0))
title(ylab="Probability", mgp=c(2.7, 0, 0))
for(i in 2:length(n.players))
  lines(dens[,1], dens[,i+1], lwd=2, col=col[i])
ymx <- apply(dens[,-1], 2, max)
xmx <- apply(dens[,-1], 2, function(a,b) max(b[a==max(a)]), dens[,1])
text(xmx+c(12, 10, 20), ymx*0.95, paste(n.players, "players"), adj=c(0, 0.5), col=col)
dev.off()

Very nice. Students of R and simulation can learn a lot from this. Thanks.