[R] cycling k times a realization of a random walk.....problems..

[Gavin Simpson]

On Fri, 2009-12-04 at 05:54 -0800, Federico Bonofiglio wrote:
> hello R-masters.
> 
> i have an R-issue here that i don't know if you'd wish to help me
> about it:
> 
> briefly i'd like to generate many (say hundred) realizations of a
> random walk, execute a few operations on each of them (mean time of
> return), and graph each realization on the same plot.
> IN OTHER WORDS I'D LIKE TO IMPOSE A LOOPING CYCLE TO THE COMMAND NOT
> THE ARGUMENT OF THE COMMAND.
> 
> for some of these questions i have already a partial answer: my main
> problem here is automatizing the process for 100 times.
> 
> my random walk generating function is:
> 
> rwalk <- function(n,p, x0=0)
>          {
>           x <- rbinom(n,1,p)
>           x[which(x==0)] <- -1
>           y<-c(x0,x)
>           y <- cumsum(y)
>            }
<snip />
 #NOTE: I'VE ALSO TRIED TO REITERATE THE WALK WITH THE rep()
>     # t<-rep(rwalk(1e+04,0.5,0),100)
>      # then of course i have the problem of splicing this
> ridicolously   #one milion long vector in 100 bits of tenthousand, and
> still goig through all the calculations and plotting i meant
> before....

One way of doing away with the loop, is to just get all the rbinoms you
want and format it appropriately, before cumsum: Try this:

rwalk <- function(n, p, x0 = 0, sim = 1)
         {
x <- rbinom(n * sim, 1, p)
x[!x] <- -1
x <- matrix(x, nrow = sim, byrow = TRUE)
x <- cbind(rep(x0, sim), x)
y <- t(apply(x, 1, cumsum))
return(y)
}

sim is the extra argument, indicating the number of walks you want. Set
it to 100 to get what you need. Each row in the returned matrix is a
random walk of the form you describe.

You can then use the apply() function to compute your statistics on each
row (walk) and aggregate the results, e.g.

> set.seed(123)
> system.time(walks <- rwalk(10000, p = 0.5, x0 = 0, sim = 100))
   user  system elapsed 
  0.747   0.010   0.759 
> tmp[1:10, 1:10] # first 10 walks, first 10 data in each
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    1    0    1    0   -1    0   -1   -2    -3
 [2,]    0   -1   -2   -3   -4   -3   -4   -5   -6    -7
 [3,]    0    1    2    1    0   -1    0   -1   -2    -3
 [4,]    0    1    0   -1    0   -1    0   -1    0    -1
 [5,]    0    1    0   -1    0   -1   -2   -3   -4    -3
 [6,]    0   -1   -2   -3   -4   -5   -4   -5   -4    -3
 [7,]    0   -1    0    1    0   -1   -2   -1    0    -1
 [8,]    0    1    0   -1   -2   -3   -2   -3   -2    -1
 [9,]    0   -1   -2   -3   -4   -3   -4   -5   -4    -3
[10,]    0   -1   -2   -3   -2   -3   -4   -5   -4    -3
> apply(tmp[1:10, ], 1, function(x) which(x == 0))
[[1]]
[1] 1 3 5 7

[[2]]
[1] 1

[[3]]
[1] 1 5 7

[[4]]
[1] 1 3 5 7 9

[[5]]
[1] 1 3 5

[[6]]
[1] 1

[[7]]
[1]  1  3  5  9 11 17 23 27 29

[[8]]
[1] 1 3

[[9]]
[1] 1

[[10]]
[1] 1

That might make it easier to work with as you only need to generate the
walks once and then develop each summary function independently. If you
start with a small example, say n = 10 and sim = 10, you should be able
to follow the computations and check they are correct.

HTH

G
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%