Why isn't my random walk simulation working correctly?

Question

I have written the following code to simulate an unbiased random walk on Z^2. With probability 1/4, the "destination" is supposed to move one unit up, left, right, or down. So I made "destination" a matrix with two columns, one for the x-coordinate and one for the y-coordinate, and increment/decrement the appropriate coordinate as according to the value of runif(1).

N_trials <- 10
N_steps <- 10
destination <- matrix(0,N_trials,2)
for(n in 1:N_steps) {

    p <- runif(1)
    if(p < 1/4) {
      destination[n,1] <- destination[n,1] - 1
    }
    else if(p < 1/2) {
      destination[n,1] <- destination[n,1] + 1
    }
    else if(p < 3/4) {
      destination[n,2] <- destination[n,2] + 1
    }
    else if(p < 1) {
      destination[n,2] <- destination[n,2] - 1
    }
  }

However, the process never seems to move out of the set {(0,0),(1,0),(-1,0),(0,1),(0,-1)}. Why is this? Is there an error in the logic of my code?

Answer 1

Rather than using loops, you can vectorize the random walk.

The idea is to first create a matrix of possible steps:

steps <- matrix(c(0,0,-1,1,-1,1,0,0),nrow = 4)

which is:

     [,1] [,2]
[1,]    0   -1
[2,]    0    1
[3,]   -1    0
[4,]    1    0

Then you can feed random subscripts into it:

steps[sample(1:4,10,replace = TRUE),]

for example will create a matrix of 9 rows where each row is randomly chosen from the steps matrix.

If you rbind this with c(0,0) as a starting position, and then take the cumulative sum (cumsum) of each column, you have your walk. You can wrap this all in a function:

rand.walk <- function(n){
  steps <- matrix(c(0,0,-1,1,-1,1,0,0),nrow = 4)
  walk <- steps[sample(1:4,n,replace = TRUE),]
  walk <-rbind(c(0,0),walk)
  apply(walk,2,cumsum)
}

For example, plot(rand.walk(1000),type = 'l') produces a graph which looks something like: