Q1. Consider creating a game map of the classic game called `Minesweeper’ when complete information of mine locations is provided. Suppose we have an \(n \times m\) binary matrix which shows where mines are located (1 = mine, 0 = nothing). A game map displays the number of mines in the surrounding 8 pixels at each pixel. Using the following two approaches, write your R functions that produce a game map when an \(n\times m\) binary matrix for the mine location is given.

i) Scan all pixels from (1, 1) to (n, m) using for loops, and check how many mines are in the surrounding 8 pixels at every location.

Before creating the for-loop code, a sample matrix is made (‘sample_mat’) so as to check whether the code works without no error.

set.seed(2011142127)

a <- sample(c(0,1), size = 35, replace = TRUE, prob = c(0.7,0.3))

sample_mat <- matrix(a, ncol = 7, byrow = TRUE)
sample_mat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0    1    0    0    0    0
## [2,]    1    0    0    0    0    0    0
## [3,]    1    1    0    0    1    1    0
## [4,]    0    0    1    1    0    0    0
## [5,]    0    0    1    0    0    0    0

The purpose of this code is to count the number of the mines surrounding a pixel. The pixel of the input matrix is 1 when there exists a mine, and 0 otherwise. So, we can obtain the number of mines around the pixel through extracting a 3*3 matrix whose center is the very pixel and summing all of the numbers of the 3*3 matrix apart from the center pixel.

However, this method returns an error without manipulating the input matrix, because the pixels in the edges can’t count the number of mines outside the matrix. So, 0 padding should be performed in all of the edges before counting the mines. Since 0 implies there is no mine in this area, the zero padding makes sense.

Let us see how the function ‘minesweeper_forloop’ work. The number of rows and columns are defined as ‘r’ and ‘c’ respectively, and a matrix ‘mat_pad’, which has 2 more columns and rows than the input matrix ‘mat’, is defined by having all of the elements be 0. The reason for creating ‘mat_pad’ is to add 0’s on every edge of the ‘mat’ matrix. Then, define another matrix ‘mat_result’ whose size is identical to the original matrix ‘mat’.

Using a nested for loop, insert each element of ‘mat’ into the ‘mat_pad’ starting from mat_pad[2,2] to mat_pad[r+1, c+1]. Then, the edges of mat_pad is still 0.

The next step is divided into 2 steps. For one thing, number 99 is put into the pixel of ‘mat_result’ where a mine exists in the corresponding pixel of ‘mat’. Second, if the pixel of ‘mat’ is empty (the number of it is 0), add up all of the surrounding numbers indicating how many numbers of mine there are and put that number to the corresponding pixel of ‘mat_result’. These processes are written through if-else statement in the second nested for-loop.

There is a nested for-loop (3rd) in the else-tatement, working as counting the surrounding 8 pixels. Inside the loop, another if-else statement exists so as to avoid the counting of the current pixel. In this problem, the value of the empty space is 0 so simply adding up all of the 9 elements of the 3*3 matrix formed by the 3rd nested for-loop has no problem. However, there should be a case where the number of the vacant pixel is assigned a non-zero value, such as NAN or NULL. Thus, the if-tatement containing ‘next’ command is added in order to handle such cases. After counting all of the surrounding digits, the total value saved in the ‘temp’ variable is transferred to the corresponding pixel of ‘mat_result’. At last, the final ‘mat_result’ matrix is returned.

minesweeper_forloop <- function(mat){
  
  r <- nrow(mat)
  c <- ncol(mat)
  mat_pad <- matrix(0, nrow = r+2, ncol = c+2)
  mat_result <- matrix(0, nrow = r, ncol = c)
  
  # zero padding
  for(i in 1:r){
    for(j in 1:c){
      mat_pad[i+1, j+1] <- mat[i,j]
    }
  }
  
  for(i in 1:r){
    for(j in 1:c){
      
      # put 99 to the mine pixel
      if(mat_pad[i+1, j+1] == 1){
        mat_result[i,j] <- 99
      }
      
      # add all of the surrounding numbers 
      else{
        temp <- 0
        
        for(k in 1:3){ # to sum up the enclosing pixels: 3*3
          for(l in 1:3){
            if(k == 2 & l == 2){
              next # not counting the current space
            }
            else{
              temp <- temp + mat_pad[i+k-1,j+l-1]
            }
          }
        }
        # indicates how many mines are in the surroundings
        mat_result[i,j] <- temp
        
      }
    }
  }
  
  return(mat_result)
}

The following code shows how the function works with the sample matrix.

sample_mat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0    1    0    0    0    0
## [2,]    1    0    0    0    0    0    0
## [3,]    1    1    0    0    1    1    0
## [4,]    0    0    1    1    0    0    0
## [5,]    0    0    1    0    0    0    0
minesweeper_forloop(sample_mat)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    2   99    1    0    0    0
## [2,]   99    4    2    2    2    2    1
## [3,]   99   99    3    3   99   99    1
## [4,]    2    4   99   99    3    2    1
## [5,]    0    2   99    3    1    0    0

ii) Consider summing up all the matrices that have shifted the original binary matrix by one pixel in all 8 directions.

(1,1) (1,2) (1,3)
(2,1) ★ (2,3)
(3,1) (3,2) (3,3)

Take the left picture as an example. When there is a mine (expressed as a black star) at (2,2) in this example, the number of all the pixels surrounding the star should rise to 1. In other words, the number of (3,3) pixel should increase by 1 because of the mine in (2,2), and the number of (2,1) pixel should increase by 1 as well for the same reason, and so on. Likewise, one mine affects all of the surrounding 8 spaces. So, the result of the originally empty pixel is same as that of a matrix made by adding up all of the 8 matrices which have moved to each 8 directions. Then, insert the number 99 to the mine pixel. At last, return the matrix whose edges are eliminated (to match the size of the input matrix).

The number of rows and columns of input matrix is designated as ‘r’ and ‘c’ like the minesweeper_forloop function. As mentioned in the question, this code does not make use of any for-oop; rather, it will use ‘rbind’ and ‘cbind’ to move the input matrix and pad it with zeros.

The following conditions will be assumed in the implementation of the code; when a matrix is received in the function, it is padded with 0’s along the edges. Moving the matrix to right is same as adding a row of 0 on the edges of top and bottom side and then adding 2 successive columns of 0’s on the left side. Moving the matrix to upper right means adding 2 consecutive columns of 0’s on the left side and then add 2 consecutive rows of 0’s on the bottom side. This methodology is applied to the other 6 directions. Then, define a matrix which is the summation of the 8 matrices. Then, put the number 99 to the pixels corresponding to the mine pixel. At last, return the matrix, omitting the 4 edges.

minesweeper_sumup <- function(mat){
  
  r <- nrow(mat)
  c <- ncol(mat)
  
  # indicate 8 directions respectively
  mat_right <- rbind(rbind(0, cbind(matrix(0, ncol = 2, nrow = r), mat)), 0)
  mat_upright <- rbind(cbind(matrix(0, ncol = 2, nrow = r), mat),
                       matrix(0, ncol = c+2, nrow = 2))
  mat_up <- cbind(cbind(0, rbind(mat, matrix(0, ncol = c, nrow = 2))), 0)
  mat_upleft <- rbind(cbind(mat, matrix(0, ncol = 2, nrow = r)),
                      matrix(0, ncol = c+2, nrow = 2))
  mat_left <- rbind(rbind(0, cbind(mat, matrix(0, ncol = 2, nrow = r))),0)
  mat_downleft <- rbind(matrix(0, ncol = c + 2, nrow = 2),
                        cbind(mat, matrix(0, ncol = 2, nrow = r)))
  mat_down <- cbind(cbind(0, rbind(matrix(0, ncol = c, nrow = 2), mat)), 0)
  mat_downright <- rbind(matrix(0, ncol = c+2, nrow = 2),
                         cbind(matrix(0, ncol = 2, nrow = r), mat))
  
  # add all of the 8 direction matrices
  mat_sum <- mat_right + mat_upright + mat_up + mat_upleft + mat_left + mat_downleft + mat_down + mat_downright
  
  # simply pad 0's to the original matrix's edge
  mat_pad <- rbind(rbind(0, cbind(cbind(0, mat), 0)), 0)
  # to find where the mines exist easily
  mat_sum[mat_pad == 1] <- 99
  
  mat_sum <- mat_sum[-c(1,r+2), -c(1,c+2)]
  # clear out the redundant edges of mat_sum matrix
  return(mat_sum)
  
}

The following code shows how the function works with the sample matrix.

sample_mat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0    1    0    0    0    0
## [2,]    1    0    0    0    0    0    0
## [3,]    1    1    0    0    1    1    0
## [4,]    0    0    1    1    0    0    0
## [5,]    0    0    1    0    0    0    0
minesweeper_sumup(sample_mat)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    2   99    1    0    0    0
## [2,]   99    4    2    2    2    2    1
## [3,]   99   99    3    3   99   99    1
## [4,]    2    4   99   99    3    2    1
## [5,]    0    2   99    3    1    0    0

Q2. Estimate \(\pi\) as follows.
Let N be the total number of uniform random numbers used for simulation in both cases. As N changes as N = 300, 600, 900, 1200, 1500, compare the two methods by calculating how close your estimates are to the true \(\pi\). Overlay time series plots (ts.plot), and keep the vertical axis within a reasonable range for clear visualization. Which one converges to \(\pi\) faster? Discuss your results.

i) Use the fact that the area of a circle of radius 1 is \(\pi\), as shown in the lecture notes. You need two uniform random numbers for one iteration of simulation.

In this case, the value of \(\pi\) is approximated based on the following assumptions. Let the (x, y) be a coordinate in a 2*2 square, and each x and y are independent and uniformly distributed. Then, the probability that the point is in a circle of radius 1 is \(\pi\)/4. Then, the number of points inside the circle divided by the number of all the points in the square is the estimation of \(\pi\)/4. So, the value π is estimated by the estimated probability multiplied by 4.

At first, generate N numbers (the number N is given in the problem) from the uniform distribution (unif(0,1)) and save it to ‘u1’. After doing some manipulation, N/2 numbers will be picked up among them. The reason for assigning N numbers first is to make the ith number of u1 in this function identical to that of the function of Q2-ii). Perform the same process and save it to ‘u2’. The ‘u1’ and ‘u2’ are modified as x = 2u1-1 and y = 2u2-1, which still follow uniform distribution unif(-1,1).

Next, define a zero vector ‘one’ whose length is N/2. When the square of ith number of x and y returns a value lower than or equal to 1, the ith value of the ‘one’ vector is shifted to 1. Repeat this process for N/2 times and add up all of the elements of the ‘one’ vector. After the summation is divided by N/2, it returns the estimated probability of being inside the circle. Thus, this estimated probability multiplied by 4 is the estimated \(\pi\).

pi_circle <- function(N){
  set.seed(2011142127) # when defining function, set seed number inside the function
  
  # u1, u2 ~ unif(0,1)
  u1 <- runif(N) # assign N random numbers and then
  u2 <- runif(N) # select N/2 of them from the first
  
  # x,y ~ unif(-1,1)
  x <- 2*u1 - 1
  y <- 2*u2 - 1
  
  one <- rep(0, N/2) #zero vector at first
  
  for(i in 1:(N/2)){
    if(x[i]^2 + y[i]^2 <= 1){# if each x and y is in the radius-one circle,
      one[i] <- 1 # add 1
    }
    else{
      next
    }
  }# 'one' vector works as an indicator whether the point is in the circle
  
  # summation of 'one' vector divided by total case 'N/2' is pi/4.
  return(4*sum(one)/(N/2))
}

This logic is included in the function pi_circle. When N is 300, 600, 900, 1200, and 1500, the estimated \(\pi\) is as follows.

pi_circle(300)
## [1] 2.88
pi_circle(600)
## [1] 3.106667
pi_circle(900)
## [1] 2.986667
pi_circle(1200)
## [1] 3.12
pi_circle(1500)
## [1] 3.178667

ii) Use the fact that the volume of a sphere of radius 1 is 4\(\pi\)/3. You need three uniform random numbers for one iteration of simulation.

The key point of the estimation here is almost same as the previous problem, apart from the fact that it is proceeded in a 3 dimensional space. Let (x,y,z) be a coordinate of 222 cube, and each x,y and z are independent and uniformly distributed (unif(-1,1)). Then the probability that the point lies in a sphere of radius 1 is (4\(\pi\)/3)/8 = \(\pi\)/6. Then, the number of points inside the sphere divided by the number of all the points in the cube is the estimated value of \(\pi\)/6. So, the value π is estimated by the estimated value multiplied by 6.

The logic of ‘pi_sphere’ function that generates the estimated \(\pi\) in 3 dimension is almost identical to that of ‘pi_circle’ function, except for some notations and additional variables. At first, 3 uniformly distributed variables u1, u2 and u3 is made (u1, u2, u3 ~ unif(0,1)), whose length is N. Then, change them so that they still follow uniform distribution (x,y, z ~ unif(-1,1)). Define an indicator vector ‘one’ with length ‘N/3’ and add 1 to the ith space if the sum of the square of ith x, y, and z is lower than or equal to 1. The added value of ‘one’ vector divided by N/3 is the estimated probability \(\pi/6\), so the estimated \(\pi\) is obtained simply by multiplying 6 to the value.

pi_sphere <- function(N){
  set.seed(2011142127)
  
  # u1, u2, u3 ~ unif(0,1)
  u1 <- runif(N) # assign N random numbers and then
  u2 <- runif(N) # choose N/3 of them from the first
  u3 <- runif(N)
  
  # x,y, z ~ unif(-1,1)
  x <- 2*u1 - 1
  y <- 2*u2 - 1
  z <- 2*u3 - 1
  
  one <- rep(0, N/3) #zero vector
  
  for(i in 1:(N/3)){# if each x, y and z are in the radius-one sphere,
    if(x[i]^2 + y[i]^2 + z[i]^2 <= 1){
      one[i] <- 1 # add 1
    }
    else{
      next
    }
  }# 'one' vector works as an indicator whether the point is in the sphere
  
  # summation of 'one' vector divided by total case 'N/3' is pi/6.
  return(6*sum(one)/(N/3))
}

This logic is included in the function pi_sphere. When N is 300, 600, 900, 1200, and 1500, the estimated \(\pi\) is as follows.

pi_sphere(300)
## [1] 2.88
pi_sphere(600)
## [1] 3.24
pi_sphere(900)
## [1] 3.14
pi_sphere(1200)
## [1] 2.985
pi_sphere(1500)
## [1] 3.144