[R] Is R the right choice for simulating first passage times of random walks?

[Dennis Murphy]

Hi Steve:

Very, very nice. Thanks for the useful Rcpp script. I'm not surprised
that a C++ version blows my humble little R function out of the water
:) I noticed that the R function ran a lot more slowly when the
sojourns were very long. It suggests that algorithms that entail
conditional iteration are quite likely to be better off written in a
compiled programming language that can communicate with R. It also
shows off the capabilities of the Rcpp package.

Best regards,
Dennis

On Sun, Jul 31, 2011 at 8:32 PM, Steve Lianoglou
<mailinglist.honeypot at gmail.com> wrote:
> Hi,
>
> I haven't been following this thread very closely, but I'm getting the
> impression that the "inner loop" that's killing you folks here looks
> quite simple (assuming it is the one provided below).
>
> How about trying to write the of this `f4` function below using the
> rcpp/inline combo. The C/C++ you will need to write looks to be quite
> trivial, let's change f4 to accept an x argument as a vector:
>
> I've defined f4 in the same way as Dennis did:
>
>> f4 <- function()
>>  {
>>     x <- sample(c(-1L,1L), 1)
>>
>>      if (x >= 0 ) {return(1)} else {
>>           csum <- x
>>           len <- 1
>>               while(csum < 0) {
>>                   csum <- csum + sample(c(-1, 1), 1)
>>                   len <- len + 1
>>                                          }     }
>>      len
>>  }
>
> Now, let's do some inline/c++ mojo:
>
> library(inline)
> inc <- "
> #include <stdio.h>
> #include <stdlib.h>
> #include <time.h>
> "
>
> fxx <-cxxfunction(includes=inc, plugin="Rcpp", body="
>  int len = 1;
>  int x = ((rand() % 2 ) == 0) ? 1 : -1;
>  int csum = x;
>
>  while (csum < 0) {
>    x = ((rand() % 2 ) == 0) ? 1 : -1;
>    len++;
>    csum = csum + x;
>  }
>
>  return wrap(len);
> ")
>
> Assuming I've faithfully translated this into c++, the timings aren't
> all that comparable.
>
> Doing 500 replicates with the pure R version:
>
> set.seed(123)
> system.time(out <- replicate(500, f4()))
>   user  system elapsed
>  31.525   0.120  32.510
>
> Doing 10,000 replicates using the fxx function doesn't even break a sweat:
>
> system.time(outxx <- replicate(10000, fxx()))
>   user  system elapsed
>  0.371   0.001   0.373
>
> range(out)
> [1]       1 1994308
>
> range(outxx)
> [1]        1 11909394
>
> Hope I'm not too off of the mark, here.
> -steve
> --
> Steve Lianoglou
> Graduate Student: Computational Systems Biology
>  | Memorial Sloan-Kettering Cancer Center
>  | Weill Medical College of Cornell University
> Contact Info: http://cbio.mskcc.org/~lianos/contact
>