5

I'm building a biased correlated random walk, and I've managed to build the RW, and bias it for westerly movement.

The issue: I need the walk to be bound on one (or all) side(s).

The current code is:

walk<-function(n.times){ 
   plot(524058:542800,2799758:2818500,type="n",
      xlab="Easting",ylab="Northing")#aren‌​a 
      y<-2815550 ##startY 
      x<-542800 #startX
      N<-4000
      E<-4000
      points(x,y,pch=16,col="red",cex=1)
      for (i in 1:n.times) {
          yi <- sample(c(N,N/2,N/4,N/8,N/12,N/16,
                      0,-N,-N/2,-N/4,-N/8,-N/12,-N/16),1)      
          xi<-sample(c(E,E/12,E/16, 0,-E,-E/2,-E/4,-E/8,-E/12,-E/16),1)       
          lines(c(x,x+xi),c(y,y+yi),col="blue")
          x<-x+xi
          y<-y+yi 
      }
   }
   iterations<-125 
   walk(iterations) 

So far the closest I've come is using

 if(y>2818500 | y<2799758 | x>542800 | x<524058)  break 

which simply stops the walk if it leaves the arena.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
Jesse001
  • 924
  • 1
  • 13
  • 37
  • 1
    I am not familiar with this but i suspect some more info might be needed to get a useful answer. Can you share more code - on how you build the RW. What do you want to happen when the RW would cross the boundary? – user20650 Oct 15 '14 at 19:36
  • When it intersects a boundary I would like it to continue the walk for the remaining number of iterations (basically just bounce off the wall and keep going) – Jesse001 Oct 15 '14 at 19:39
  • when it bounces off the wall does it reflect at an opposite angle or does it return the way it came? – Ben Bolker Oct 15 '14 at 19:44
  • The walk code itself is simple: > walk<-function(n.times){ > plot(524058:542800,2799758:2818500,type="n",xlab="Easting",ylab="Northing")#arena >y<-2815550 ##startY >x<-542800 #startX >N<-4000 >E<-4000 >points(x,y,pch=16,col="red",cex=1) >for (i in 1:n.times) >{yi <- sample(c(N,N/2,N/4,N/8,N/12,N/16, 0,-N,-N/2,-N/4,-N/8,-N/12,-N/16),1) >xi<-sample(c(E,E/12,E/16, 0,-E,-E/2,-E/4,-E/8,-E/12,-E/16),1) >lines(c(x,x+xi),c(y,y+yi),col="blue") >x<-x+xi >y<-y+yi >iterations<-125 >walk(iterations) – Jesse001 Oct 15 '14 at 19:47
  • Sorry, I don't know how to make code read pretty in this forum yet – Jesse001 Oct 15 '14 at 19:48
  • @Ben Bolker: I'd like it to reflect randomly. Basically hitting a wall would cause it to invoke the sample function again – Jesse001 Oct 15 '14 at 19:49
  • try `repeat { ... ; if () break }` – Ben Bolker Oct 15 '14 at 20:02
  • Thank you very much for your suggestion. For some reason I can't seem to get it to recognize that line. It keeps saying: "Error:'...'used in an incorrect context". when I removed the and replaced it with a logical operator I got the results below. I tried modifying the code in two ways: 1) repeat{if)xi>5424058)break} which puts it in an infinite loop 2) repeat{if(xi<525000)break} which had no effect on the walk did I misunderstand your comment? If so, how would I better implement it? Again I appreciate your assistance – Jesse001 Oct 15 '14 at 20:55
  • the `...` is just shorthand for "put some appropriate stuff here" – Ben Bolker Oct 15 '14 at 21:30

1 Answers1

9

A slightly cleaned-up version of the function: the primary change is the addition of a repeat {} loop that will re-pick the step until the new location is within the limits (could also use a while() {} loop).

update: didn't read the problem statement carefully enough, forgot the bias. This code incorporates the bias in the same way that the OP's code does. For N-S movement the mean step length is 0; for the E-W movement, by leaving out some of the positive step possibilities, we get mean(steps.x) equal to -0.0875; since the step possibilities are uniformly sampled, the walk drifts to the left by an average of 0.0875*stepsize[1] units per step.

 walk <- function(n.times=125,
               xlim=c(524058,542800),
               ylim=c(2799758,2818500),
               start=c(542800,2815550),
               stepsize=c(4000,4000)) {
    ## blank plot of arena
    plot(c(0,0),type="n",xlim=xlim,ylim=ylim,
           xlab="Easting",ylab="Northing") 
    ## extract starting point
    x <- start[1]
    y <- start[2]
    ## define potential step sizes
    steps <- 1/c(1,2,4,8,12,16)
    ## all possible positive or negative steps for N-S movement
    steps.y <- c(steps,-steps,0)
    ## bias E-W movement by leaving out some positive steps
    steps.x <- c(steps[c(1,5,6)],-steps,0)
    ## plot starting location
    points(x,y,pch=16,col="red",cex=1)
    for (i in 1:n.times) {
        repeat {
           ## pick jump sizes
           xi <- stepsize[1]*sample(steps.x,1)
           yi <- stepsize[2]*sample(steps.y,1)
           ## new candidate locations
           newx <- x+xi
           newy <- y+yi
           ## IF new locations are within bounds, then
           ##    break out of the repeat{} loop (otherwise
           ##    try again)
           if (newx>xlim[1] && newx<xlim[2] &&
               newy>ylim[1] && newy<ylim[2]) break
        }
        lines(c(x,newx),c(y,newy),col="blue") ## draw move
        ## set new location to candidate location
        x <- newx
        y <- newy
    }
}
set.seed(101)
walk(1000)

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Fantastic Ben thank you very much! If I could inquire a bit (I'm sorry but I don't read code that well). In the step function I'm failing to see where the bias comes in/how it works. I'll need to adjust this quite a few times so it's critical I figure this out. your help is greatly appreciated! – Jesse001 Oct 16 '14 at 17:29
  • I just incorporated the bias in the same way that your code did. For N-S movement the mean step length is 0; for the E-W movement, by leaving out some of the positive step possibilities, we get `mean(steps.x)` equal to -0.0875; since the step possibilities are uniformly sampled, the walk drifts to the left by an average of 0.0875*stepsize[1] units per step. – Ben Bolker Oct 16 '14 at 17:39
  • This is amazing. Ben, may as ask you to comment your code. I'm not understanding the process inside the loop. – Paulo E. Cardoso Oct 16 '14 at 20:18
  • Ben, Thank you for your kindness. What happens when the `if` statement `== T`? – Paulo E. Cardoso Oct 16 '14 at 22:56