There does appear to be a difference between the native Rf_rnorm and the R rnorm function when the sample sizes get large. In fact, the Rf_rnorm appears to be about twice as fast (using benchmark). Not surprisingly, doing a loop in R is about 50 times slower.
From: Jonathan Olmsted [mailto:jpolms...@gmail.com] Sent: Tuesday, September 25, 2012 4:29 PM To: Jeffrey Pollock Cc: Goldfeld, Keith; rcpp-devel@lists.r-forge.r-project.org Subject: Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array Jeff and List, I was interested in the same fairly recently. The full shebang is here <http://www.rochester.edu/college/gradstudents/jolmsted/blog/2012/08/22/rcpp-rng-performance/>. I, unlike you, found some difference between the various approaches. If I had to guess why we come up with different results, I'd say that the time for any one call for any of your inline functions is 90% overhead and the differences between the various approaches, while present, is on a vastly different (and smaller) scale than the overhead itself due to only 5 draws being taken. Of course, how any of these timings inform one's approach depends on the specifics of how they'll set up their analysis. For what it is worth, instead of an explicit "as<>()" conversion, I usually pull out the first element of the returned NumericVector. That is, rn[i] = as<double>(rnorm(1, 0.0, 1.0)); becomes rn[i] = rnorm(1, 0.0, 1.0)[0]; I haven't ever timed this, though. JPO ------------------------------------------------------------------------- J. P. Olmsted j.p.olms...@gmail.com<mailto:j.p.olms...@gmail.com> http://about.me/olmjo ------------------------------------------------------------------------- On Tue, Sep 25, 2012 at 4:08 PM, Jeffrey Pollock <jeffpollo...@gmail.com<mailto:jeffpollo...@gmail.com>> wrote: I was interested to see if there was any real speed difference between the different methods suggested, and it looks like... there isn't... library(inline) library(Rcpp) library(rbenchmark) fun1 <- cxxfunction(body = ' RNGScope scope; NumericVector rn(5); for (int i = 0; i < 5; i++) { rn[i] = as<double>(rnorm(1, 0.0, 1.0)); } return rn; ', plugin = "Rcpp") fun2 <- cxxfunction(body = ' RNGScope scope; NumericVector rn(5); for (int i = 0; i < 5; i++) { rn[i] = Rf_rnorm(0.0, 1.0); } return rn; ', plugin = "Rcpp") fun3 <- cxxfunction(body = ' RNGScope scope; NumericVector rn(5); for (int i = 0; i < 5; i++) { rn[i] = norm_rand(); } return rn; ', plugin = "Rcpp") fun4 <- cxxfunction(body = ' RNGScope scope; return rnorm(5, 0.0, 1.0); ', plugin = "Rcpp") set.seed(1) print(fun1()) set.seed(1) print(fun2()) set.seed(1) print(fun3()) set.seed(1) print(fun4()) benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications = 1e6L) output; > set.seed(1) > print(fun1()) [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 > set.seed(1) > print(fun2()) [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 > set.seed(1) > print(fun3()) [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 > set.seed(1) > print(fun4()) [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 > benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications = > 1e6L) test replications elapsed relative user.self sys.self user.child sys.child 3 fun3() 1000000 7.15 1.000 7.01 0.00 NA NA 4 fun4() 1000000 7.33 1.025 7.27 0.01 NA NA 2 fun2() 1000000 7.53 1.053 7.41 0.00 NA NA 1 fun1() 1000000 7.99 1.117 7.91 0.00 NA NA On Tue, Sep 25, 2012 at 8:19 PM, Goldfeld, Keith <keith.goldf...@nyumc.org<mailto:keith.goldf...@nyumc.org>> wrote: Yes - I agree that iterating is not the best way to go, but I am using this for an agent based simulation where iteration is really my only option. I won't be generating more than a few thousand records, so the time shouldn't be much of a factor. Of course, in R it is painfully slow - so that is why I am going the Rcpp route. -----Original Message----- From: dmba...@gmail.com<mailto:dmba...@gmail.com> [mailto:dmba...@gmail.com<mailto:dmba...@gmail.com>] On Behalf Of Douglas Bates Sent: Tuesday, September 25, 2012 3:13 PM To: Rodney Sparapani Cc: Goldfeld, Keith; rcpp-devel@lists.r-forge.r-project.org<mailto:rcpp-devel@lists.r-forge.r-project.org> Subject: Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array On Tue, Sep 25, 2012 at 1:53 PM, Rodney Sparapani <rspar...@mcw.edu<mailto:rspar...@mcw.edu>> wrote: > On 09/25/2012 01:37 PM, Goldfeld, Keith wrote: > >>> code<- 'Rcpp::RNGScope scope; >> >> >> + Rcpp::NumericVector rn(5); >> >> + for (int i=0; i< 5; i++) { >> >> + rn(i) = rnorm(1,0,1); >> >> + } >> >> + return Rcpp::wrap(rn); >> >> + ' >> >>> >> >>> fun<- cxxfunction(body=code, plugin="Rcpp") > > > You just need an explicit type conversion like so: > > funk2 <- cxxfunction(body='RNGScope scope; > NumericVector rn(5); > for (int i=0; i < 5; i++) rn[i] = as<double>(rnorm(1,0,1)); > return wrap(rn);', > includes="using namespace Rcpp;\n", plugin="Rcpp") > > funk2() That works but you wouldn't want to use it as a general model because you are creating vectors of length 1 then dereferencing the first element in each of these vectors and copying it to another vector. Obviously with n=5 the difference in execution time will be minimal but with n=5 million it won't be. Using the Rcpp sugar function rnorm will be the easiest approach unless, like me, you get a little queasy when using Rcpp sugar functions just because it is so hard to pick them apart and see what is actually happening. I would probably end up going back to the R API and calling norm_rand or Rf_rnorm. That latter returns a double from two double arguments, mu and sigma. _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org<mailto:Rcpp-devel@lists.r-forge.r-project.org> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org<mailto:Rcpp-devel@lists.r-forge.r-project.org> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
