Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
193 lines (171 sloc) 5.29 KB
# ***********************************
# 100 Prisoners, 100 lines of code
# Author: Benjamin Tovar
# Date: August 1, 2014
#
# post: http://tata-box-blog.blogspot.com/2014/08/update-100-prisoners-100-lines-of-code.html
#
# Based on the post: http://www.r-bloggers.com/100-prisoners-100-lines-of-code/
# ***********************************
# load the libraries
library(ggplot2)
library(gridExtra)
# set the variable values
n.prisoners=100
n.boxes=100
iters=1000
# ************************
# run the model
# this means, the simulation where
# the prisoners have a strategy
# (they will open their first box based on their
# own number, and if the box does not contain their
# number, they will move to the box based on the number
# inside the previous opened one)
# ************************
results <- run.model(n.prisoners,n.boxes,iters)
# ************************
# run the random model
# this means, the simulation where
# the prisoners do not have a strategy
# (they will open every box they can by random)
# ************************
results.random <- run.random.model(n.prisoners,n.boxes,iters)
# **********************
# PLOT THE HISTOGRAMS
# OF THE RESULTS
# **********************
# plot the distribution of the random model
p1 <- ggplot(data.frame(model=rep("random.model",iters),
values=results.random)) +
aes(x=values) +
geom_histogram(stat="bin",binwidth=1,
fill="red",colour="red",alpha=0.1) +
labs(title="RANDOM MODEL: Number of prisoners that succeeded the test",
x="Success event",
y="Frequency")
# plot the distribution of the model
p2 <- ggplot(data.frame(model=rep("model",iters),
values=results)) +
aes(x=values) +
geom_histogram(stat="bin",binwidth=1,
fill="lightblue",colour="blue",alpha=0.3) +
labs(title="MODEL: Number of prisoners that succeeded the test",
x="Success event",
y="Frequency")
# merge the plots
grid.arrange(p1, p2)
# **********************
# PLOT THE DENSITY OF THE RESULTS
# **********************
ggplot(data.frame(model=c(rep("model",iters),rep("random.model",iters)),
values=c(results,results.random))) +
aes(x=values) +
geom_density(aes(fill=as.factor(model),colour=as.factor(model)),alpha=0.3) +
labs(title="Density plot of number of prisoners that succeeded retrieving their number",
x="Success event",
y="Density")
# **********************
# PLOT THE HISTOGRAMS
# **********************
ggplot(data.frame(model=c(rep("model",iters),rep("random.model",iters)),
values=c(results,results.random))) +
aes(x=values) +
geom_histogram(position="dodge",stat="bin",binwidth=1,
aes(fill=as.factor(model),
colour=as.factor(model)),alpha=0.3) +
geom_freqpoly(binwidth=1,aes(fill=as.factor(model),colour=as.factor(model))) +
labs(title="Histograms of number of prisoners that succeeded retrieving their number",
x="Success event",
y="Frequency")
# **********************
# PLOT THE BOXPLOT
# **********************
dat <- data.frame(value=c(results,results.random),
strategy=c(rep("model",iters),rep("random.model",iters)))
ggplot(dat) +
aes(x=strategy,y=value) +
geom_boxplot(aes(fill=as.factor(strategy),
colour=as.factor(strategy)),alpha=0.5) +
geom_jitter(colour="darkgreen",alpha=0.15) +
labs(title="Boxplot of the models",
x="Model",
y="Success event")
###################
# FUNCTIONS
###################
run.model <- function(n.prisoners=100,n.boxes=100,iters=1000) {
results = rep(0,iters)
# Labels for our prisoners
prisoners = 1:n.prisoners
# run the model
for(i in 1:iters) {
# A random permutation:
boxes = sample(1:n.boxes,n.boxes)
# Track how many "winners" we have
foundIt = 0
# Main loop over the prisoners
for(prisoner in prisoners) {
# Track the prisoners path
path = c(prisoner)
tries = 1
# Look first in the box that matches your own number
inBox = boxes[prisoner]
while(tries < 50) {
path = c(path, inBox)
if(inBox == prisoner) {
foundIt = foundIt + 1
break;
} else {
# Follow that number to the next box
inBox = boxes[inBox]
}
tries = tries+1
}
}
# How many prisoners found their numbers?
results[i] = foundIt
}
return(results)
}
run.random.model <- function(n.prisoners=100,n.boxes=100,iters=1000) {
results = rep(0,iters)
# Labels for our prisoners
prisoners = 1:n.prisoners
# run the model
for(i in 1:iters) {
# A random permutation:
boxes = sample(1:n.boxes,n.boxes)
# Track how many "winners" we have
foundIt = 0
# Main loop over the prisoners
for(prisoner in prisoners) {
# create a copy of the boxes object to
# keep an index of the opened boxes for
# each prisoner
boxes.temp <- boxes
# Track the prisoners path
path = c(prisoner)
tries = 1
# Look first in the box that matches your own number
inBox = sample(boxes,1)
while(tries < 50) {
path = c(path, inBox)
if(inBox == prisoner) {
foundIt = foundIt + 1
break;
} else {
# choose randomly among any of the remaining boxes
# delete that box from the box pool
boxes.temp <- boxes.temp[boxes.temp!=inBox]
# choose among the remaining boxes
inBox = sample(boxes.temp,1)
}
tries = tries+1
}
}
# How many prisoners found their numbers?
results[i] = foundIt
}
return(results)
}