Random Numbers and Game of Life: Exercises and Solutions in R
The following exercises are from the book “Introduction to scientific programming and simulation using R” by O. Jones, R. Maillardet and A. Robinson (2009).
Solution file can be downloaded here.
Chapter 5, section 7, exercise 3, parts a, b, and c
In this question we simulate the rolling of a die. To do this we use the function runif(1), which returns a ‘random’ number in the range (0,1). To get a random integer in the range {1, 2, 3, 4, 5, 6}, we use ceiling(6*runif(1)), or if you prefer, sample(1:6,size=1) will do the same job.
- Suppose that you are playing the gambling game of the Chevalier de Mere. That is, you are betting that you get at least one six in 4 throws of a die. Write a program that simulates one round of this game and prints out whether you win or lose. Check that your program can produce a different result each time you run it.
win <- FALSE
for (i in 1:4) {
if (sample(1:6, size = 1) == 6) {
win <- TRUE
}
}
if (win) {
print("win")
} else {
print("lose")
}## [1] "win"- Turn the program that you wrote in part (a) into a function
sixes, which returnsTRUEif you obtain at least one six in n rolls of a fair die, and returnsFALSEotherwise. That is, the argument is the number of rolls n, and the value returned isTRUEif you get at least one six andFALSEotherwise. How would you give n the default value of 4?
sixes <- function(n = 4) {
win <- FALSE
for (i in 1:n) {
if (sample(1:6, size = 1) == 6) {
win <- TRUE
}
}
return(win)
}- Now write a program that uses your function
sixesfrom part (b), to simulate N plays of the game (each time you bet that you get at least 1 six in n rolls of a fair die). Your program should then determine the proportion of times you win the bet. This proportion is an estimate of the probability of getting at least one 6 in n rolls of a fair die. Run the program for n = 4 and N = 100, 1000, and 10000, conducting several runs for each N value. How does the variability of your results depend on N? The probability of getting no 6’s in n rolls of a fair die is \((5/6)^n\), so the probability of getting at least one is \(1 - (5/6)^n\). Modify your program so that it calculates the theoretical probability as well as the simulation estimate and prints the difference between them. How does the accuracy of your results depend on N? You may find thereplicatefunction useful here.
sixes_rep <- function(n = 4, r) {
# n is the number of dice
# r is the number of replicates (capital N in the text)
# observed probability of getting at least one 6
obs <- mean(replicate(r, sixes(n)))
# theoretical probability of getting at least one 6
theor <- 1 - (5/6)^n
difference <- obs - theor
# print the results
cat("Theoretical prob of at least one six in", n, "dice is:", theor, "\n")
cat("Empirical prob of at least one six after", r,"replications is:", obs, "\n")
cat("The difference is", difference ,"\n")
}
# test the function
set.seed(1)
sixes_rep(4, 100)## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 100 replications is: 0.48
## The difference is -0.03774691## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 100 replications is: 0.57
## The difference is 0.05225309## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 100 replications is: 0.53
## The difference is 0.01225309## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 1000 replications is: 0.526
## The difference is 0.008253086## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 1000 replications is: 0.511
## The difference is -0.006746914## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 1000 replications is: 0.506
## The difference is -0.01174691## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 10000 replications is: 0.5176
## The difference is -0.0001469136## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 10000 replications is: 0.5191
## The difference is 0.001353086## Theoretical prob of at least one six in 4 dice is: 0.5177469
## Empirical prob of at least one six after 10000 replications is: 0.5103
## The difference is -0.007446914As N increases, the variability of the results decreases and the accuracy increases.
Chapter 5, section 7, exercise 4
Consider the following program and its output:
# clear the workspace
rm(list=ls())
random.sum <- function(n) {
# sum of n random numbers
x[1:n] <- ceiling(10*runif(n))
cat("x:", x[1:n], "\n")
return(sum(x))
}
x <- rep(100, 10)
show(random.sum(10))## x: 3 5 10 1 1 8 5 2 5 7
## [1] 47## x: 4 5 9 4 4
## [1] 526Explain what is going wrong and how you would fix it.
Answer: The code is supposed to take the input argument n, generate n random numbers and then return the sum of those random numbers. However, since there already is an object named x in the environment, when n is less than the length of that object (in this case less than 10), the function replaces the first n numbers in x by the generated numbers but when it adds up all the numbers in x, it adds up not only those n numbers, but all numbers in x.
One way to fix the function:
random.sum <- function(n) {
x[1:n] <- ceiling(10*runif(n))
cat("x:", x[1:n], "\n")
return(sum(x[1:n]))
}
# show that it works
x <- rep(100, 10)
random.sum(10)## x: 7 1 10 10 9 6 5 1 3 4## [1] 56## x: 8 2 5 1 9## [1] 25Chapter 5, section 7, exercise 6
The Game of Life is a cellular automaton and was devised by the mathematician J.H. Conway in 1970. It is played on a grid of cells, each of which is either alive or dead. The grid of cells evolves in time and each cell interacts with its eight neighbours, which are the cells directly adjacent horizontally, vertically, and diagonally.
At each time step cells change as follows:
- A live cell with fewer than two neighbours dies of loneliness.
- A live cell with more than three neighbours dies of overcrowding.
- A live cell with two or three neighbours lives on to the next generation.
- A dead cell with exactly three neighbours comes to life.
The initial pattern constitutes the first generation of the system. The second generation is created by applying the above rules simultaneously to every cell in the first generation: births and deaths all happen simultaneously. The rules continue to be applied repeatedly to create further generations.
Theoretically the Game of Life is played on an infinite grid, but in practice we use a finite grid arranged as a torus. That is, if you are in the left-most column of the grid then your left-hand neighbours are in the right-most column, and if you are in the top row then your neighbours above are in the bottom row.
Here is an implementation of the Game of Life in R. The grid of cells is stored in a matrix A, where A[i,j] is 1 if cell (i, j) is alive and 0 otherwise.
neighbours <- function(A, i, j, n) {
# A is an n*n 0-1 matrix
# calculate number of neighbours of A[i,j]
.
.
.
}
# grid size
n <- 50
# initialise lattice
A <- matrix(round(runif(n^2)), n, n)
finished <- FALSE
while (!finished) {
# plot
plot(c(1,n), c(1,n), type = "n", xlab = "", ylab = "")
for (i in 1:n) {
for (j in 1:n) {
if (A[i,j] == 1) {
points(i, j)
}
}
}
# update
B <- A
for (i in 1:n) {
for (j in 1:n) {
nbrs <- neighbours(A, i, j, n)
if (A[i,j] == 1) {
if ((nbrs == 2) | (nbrs == 3)) {
B[i,j] <- 1
} else {
B[i,j] <- 0
}
} else {
if (nbrs == 3) {
B[i,j] <- 1
} else {
B[i,j] <- 0
}
}
}
}
A <- B
## continue?
#input <- readline("stop? ")
#if (input == "y") finished <- TRUE
}Note that this program contains an infinite loop! To stop it you will need to use the escape or stop button (Windows or Mac) or control-C (Unix). Alternatively, uncomment the last two lines. To get the program to run you will need to complete the function neighbours(A, i, j, n), which calculates the number of neighbours of cell (i, j).
Solution:
glidergun <- function(n) {
# sets the initial n*n lattice for a `glider gun`
# assumes n >= 40
A <- matrix(0, n, n)
A[31,13] <- 1
A[31,14] <- 1
A[32,12] <- 1
A[32,16] <- 1
A[33,11] <- 1
A[33,17] <- 1
A[33,25] <- 1
A[34,1] <- 1
A[34,2] <- 1
A[34,11] <- 1
A[34,15] <- 1
A[34,17] <- 1
A[34,18] <- 1
A[34,23] <- 1
A[34,25] <- 1
A[35,1] <- 1
A[35,2] <- 1
A[35,11] <- 1
A[35,17] <- 1
A[35,21] <- 1
A[35,22] <- 1
A[36,12] <- 1
A[36,16] <- 1
A[36,21] <- 1
A[36,22] <- 1
A[36,35] <- 1
A[36,36] <- 1
A[37,13] <- 1
A[37,14] <- 1
A[37,21] <- 1
A[37,22] <- 1
A[37,35] <- 1
A[37,36] <- 1
A[38,23] <- 1
A[38,25] <- 1
A[39,25] <- 1
return(A)
}
conway_plot <- function(A){
n <- nrow(A)
plot(c(1,n), c(1,n), type = "n", xlab = "", ylab = "")
for (i in 1:n) {
for (j in 1:n) {
if (A[i,j] == 1) {
points(i, j)
}
}
}
}
conway_update <- function(A){
n <- nrow(A)
B <- A
for (i in 1:n) {
for (j in 1:n) {
nbrs <- neighbors(A, i, j, n)
if (A[i,j] == 1) {
if ((nbrs == 2) | (nbrs == 3)) {
B[i,j] <- 1
} else {
B[i,j] <- 0
}
} else {
if (nbrs == 3) {
B[i,j] <- 1
} else {
B[i,j] <- 0
}
}
}
}
A <- B
A
}
# neighbors function
neighbors <- function(A, i, j, n) {
nbrs <- 0
x <- ((i-1):(i+1) - 1) %% n + 1
y <- ((j-1):(j+1) - 1) %% n + 1
for(k in x) {
for(l in y) {
if (A[k, l] == 1) {
nbrs <- nbrs + 1
}
}
}
return(nbrs - A[i, j])
}# If you wish to test the neighbors function, you can use this code
# which will run an infinite loop updating the plot until you enter 'y' to stop
A <- glidergun(50)
finished <- FALSE
while (!finished) {
# plot
conway_plot(A)
# update
A <- conway_update(A)
# continue
input <- readline("stop? ")
if (input == "y") finished <- TRUE
}The code below generates the plot after 100 updates of Conway’s game of life.
# grid size
n <- 50
# initialize lattice
A <- glidergun(n)
# code to run and stop after k iterations
k <- 100
for(i in 1:k) {
# update
A <- conway_update(A)
}
# plot
conway_plot(A)