public
Last active

R simulation of chutes and ladders

  • Download Gist
chutes_and_ladders.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
# 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
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
# 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
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
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.

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.