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>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] 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 > 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> > 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 > 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