#
#
# Please participate in the R&SS Client Feedback Survey.
# https://unt.az1.qualtrics.com/SE/?SID=SV_diLLVU9iuJZN8C9
#
###### Random Jumping Fleas -- how many empty boxes after 50 bells.
#
# So, here's the setup:
# Taken from: https://www.r-bloggers.com/flea-circus/
#
# An old riddle found on X validated asking for Monte Carlo resolution but originally given
# on Project Euler:
# A 30×30 grid of squares contains 30² fleas, initially one flea per square. When a bell is
# rung, each flea jumps to an adjacent square at random. What is the expected number of
# unoccupied squares after 50 bell rings, up to six decimal places?
# The debate on X validated is whether or not a Monte Carlo resolution is feasible. Up to six
# decimals, certainly not. But with some lower precision, certainly.
#
# It's important to note, the original problem above specifies (thankfully) fleas can only 'jump'
# to one of 4 adjacent boxes / cells (i.e. no diagonal jumps; which makes things much easier).
#
# So, we're going to start with a 5 X 5 grid of squares...scaling up will presumably be easy later.
grid <- data.frame(matrix(rep(1, 49), ncol = 7))
grid <- matrix(rep(1, 49), ncol = 7)
grid[1,] <- NA
grid[7,] <- NA
grid[,1] <- NA
grid[,7] <- NA
fleas.df <- data.frame(id=seq(1:25), x=rep(seq(2,6),5),
y=c(rep(2,5), rep(3,5), rep(4,5), rep(5,5), rep(6,5)),
new.x=rep(0,25), new.y=rep(0,25))
grid
fleas.df
jumps <- c("n","s","e","w")
bells <- 10
while(bells > 0){
j <- 1
while(j < nrow(fleas.df)){
move <- sample(jumps, 1)
if(move == "e"){fleas.df[j,4] <- fleas.df[j,2] + 1
fleas.df[j,5] <- fleas.df[j,3]}
if(move == "w"){fleas.df[j,4] <- fleas.df[j,2] - 1
fleas.df[j,5] <- fleas.df[j,3]}
if(move == "n"){fleas.df[j,5] <- fleas.df[j,3] + 1
fleas.df[j,4] <- fleas.df[j,2]}
if(move == "s"){fleas.df[j,5] <- fleas.df[j,3] - 1
fleas.df[j,4] <- fleas.df[j,2]}
if(is.na(grid[fleas.df[j,4], fleas.df[j,5]]) == FALSE){
grid[fleas.df[j,2],fleas.df[j,3]] <- grid[fleas.df[j,2],fleas.df[j,3]] - 1
grid[fleas.df[j,4],fleas.df[j,5]] <- grid[fleas.df[j,4],fleas.df[j,5]] + 1
fleas.df[j,2] <- fleas.df[j,4]
fleas.df[j,3] <- fleas.df[j,5]
fleas.df[j,4] <- 0
fleas.df[j,5] <- 0
j <- j + 1}
}
bells <- bells - 1
}
fleas.df
grid
length(which(as.vector(na.omit(as.vector(grid))) == 0))
# Clean up.
ls()
rm(bells, fleas.df, grid, j, jumps, move)
ls()
###################################################################################################
#
# Now attempt to put it all into a function....
flea.circus <- function(bells = 50, boxes.rc = 30){
jumps <- c("n","s","e","w")
N <- boxes.rc*boxes.rc
n <- boxes.rc + 2
grid <- matrix(rep(1,n*n), ncol = n)
grid[1,] <- rep(NA,n)
grid[n,] <- rep(NA,n)
grid[,1] <- rep(NA,n)
grid[,n] <- rep(NA,n)
n1 <- 2
n2 <- boxes.rc + 1
fleas.df <- data.frame(matrix(rep(0,N*4), ncol = 4))
fleas.df[,1] <- rep(seq(n1,n2), n2-1)
fleas.df[,2] <- as.vector(t(matrix(rep(seq(n1,n2), n2-1), ncol = n2-1)))
while(bells > 0){
j <- 1
while(j < nrow(fleas.df)){
move <- sample(jumps, 1)
if(move == "e"){fleas.df[j,3] <- fleas.df[j,1] + 1
fleas.df[j,4] <- fleas.df[j,2]}
if(move == "w"){fleas.df[j,3] <- fleas.df[j,1] - 1
fleas.df[j,4] <- fleas.df[j,2]}
if(move == "n"){fleas.df[j,4] <- fleas.df[j,2] + 1
fleas.df[j,3] <- fleas.df[j,1]}
if(move == "s"){fleas.df[j,4] <- fleas.df[j,2] - 1
fleas.df[j,3] <- fleas.df[j,1]}
if(is.na(grid[fleas.df[j,3], fleas.df[j,4]]) == FALSE){
grid[fleas.df[j,1],fleas.df[j,2]] <- grid[fleas.df[j,1],fleas.df[j,2]] - 1
grid[fleas.df[j,3],fleas.df[j,4]] <- grid[fleas.df[j,3],fleas.df[j,4]] + 1
fleas.df[j,1] <- fleas.df[j,3]
fleas.df[j,2] <- fleas.df[j,4]
fleas.df[j,3] <- 0
fleas.df[j,4] <- 0
j <- j + 1}
}
bells <- bells - 1
}
output <- length(which(as.vector(na.omit(as.vector(grid))) == 0))
return(output)
}
system.time(flea.circus(bells = 50, boxes.rc = 30))
# That's roughly 11 seconds per run; so if we were to do a monte carlo simulation to build
# a distribution of the number of empty boxes, say n = 5000, that would take roughly
# 15 hours to complete (obviously, 60 seconds in a minute and 60 minutes in an hour).
(11*5000) / (60*60)
# Therefore, I don't recommend doing it with the usually suggested n = 5000. But if you've got
# a day to wait and perhaps an old or leftover machine you don't use on a regular basis...
n <- 5000
mc.results <- rep(0,n)
for (k in 1:n){
mc.results[k] <- flea.circus(bells = 50, boxes.rc = 30)
}; rm(n,k)
summary(mc.results)
hist(mc.resutls)
# Instead, let's try it with only n = 250 which should take about 45 minutes....while
# printing the iteration number (k); and then after 2 iterations also printing the
# summary and a histogram of the results.
mc.results <- as.vector(0)
for (k in 1:250){
mc.results <- c(mc.results, flea.circus(bells = 50, boxes.rc = 30))
if(k == 1){mc.results <- mc.results[-1]}
print(k)
if(k > 2){
print(summary(mc.results))
hist(mc.results)
}
}; rm(k)
summary(mc.results)
hist(mc.resutls)
# End script.