A 3D Random walk simulation I was referring to the book Computational Physics by Nicholas J. Giordano and was studying the simple Random Walk model. (Pg. 169 - 173)

I consider a 3D Random Walk where the probability of moving in the positive direction is p and in the negative direction is q, in all the 3 coordinates and p + q = 1. In my model I assume:


$$ p_x=q_x =0.5,\space p_y=q_y=0.5 \space and \space p_z=q_z=0.5$$
$$l_x = l_y=l_z=\pm 1$$


where l is the step length. In my program I consider total number of walkers = N and total number of random steps taken by each walker is s. My goals are:


Plot a scatter plot : $$&lt;x_s^2&gt; \space vs. \space t$$
Find the slope from the regression line.


Here is my R code:

randomWalk &lt;- function(N, s) {
  S &lt;- 1:s
  av &lt;- numeric(s)
  for (r in 1:N) {
    x &lt;- 0
    y &lt;- 0
    z &lt;- 0
    for (i in 1:s) {
      if (runif(1)&lt;0.5) x &lt;- x+1
      else x &lt;- x-1
      if (runif(1)&lt;0.5) y &lt;- y+1
      else y &lt;- y-1
      if (runif(1)&lt;0.5) z &lt;- z+1
      else z &lt;- z-1
      av[i] &lt;- av[i] + (x^2+y^2+z^2)
    }
  }
  av &lt;- av/N
  plot(S, av,xlab="No. of steps(t)", ylab="&lt;r^2&gt;", main="Random Walk",cex=0.4, pch=15)
  abline(lm(av~(S)-1))
  print(lm(av~(S)-1))
}
randomWalk(500, 100) 


In the code I consider N=500 and s=100.

The plot I get is  

with value of slope approximately 3. But in page 173 there is a plot with the same condition but slope approximately 1.

My question is, Is there something wrong with my code or am I misunderstanding the plot given in page 173? A review will be helpful.
Quoting page 173:


  Another obvious generalization is to allow the walker to move in three dimensions [...] For this simulation we have restricted the steps to be of unit length along either ±x, ±y, or ±z.


I think the author meant that at each time step the walker will walk in one of three random directions (x or y or z), not all three. Without making drastic changes to your code, you could replace this section:

 if (runif(1)&lt;0.5) x &lt;- x+1 else x &lt;- x-1
 if (runif(1)&lt;0.5) y &lt;- y+1 else y &lt;- y-1
 if (runif(1)&lt;0.5) z &lt;- z+1 else z &lt;- z-1


with the following:

  direction &lt;- sample(c("x", "y", "z"), 1)
  if (direction == "x") { if (runif(1)&lt;0.5) x &lt;- x+1 else x &lt;- x-1 }
  if (direction == "y") { if (runif(1)&lt;0.5) y &lt;- y+1 else y &lt;- y-1 }
  if (direction == "z") { if (runif(1)&lt;0.5) z &lt;- z+1 else z &lt;- z-1 }


With this change, the slope is approximately 1, as shown in Figure 7.11.