Skip to content

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
R simulation of chutes and ladders
# 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)
# 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()
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()
@jameshunterbr

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.