I was interested to see if there was any real speed difference between the different methods suggested, and it looks like... there isn't... <br><br>library(inline)<br>library(Rcpp)<br>library(rbenchmark)<br><br>fun1 <- cxxfunction(body = '<br>
RNGScope scope; <br> NumericVector rn(5);<br> for (int i = 0; i < 5; i++) {<br> rn[i] = as<double>(rnorm(1, 0.0, 1.0)); <br> }<br> return rn;<br>
', plugin = "Rcpp")<br><br>fun2 <- cxxfunction(body = '<br> RNGScope scope; <br> NumericVector rn(5);<br> for (int i = 0; i < 5; i++) {<br>
rn[i] = Rf_rnorm(0.0, 1.0);<br> }<br> return rn;<br> ', plugin = "Rcpp")<br><br>fun3 <- cxxfunction(body = '<br> RNGScope scope; <br>
NumericVector rn(5);<br> for (int i = 0; i < 5; i++) {<br> rn[i] = norm_rand();<br> }<br> return rn;<br> ', plugin = "Rcpp")<br>
<br>fun4 <- cxxfunction(body = '<br> RNGScope scope; <br> return rnorm(5, 0.0, 1.0);<br> ', plugin = "Rcpp")<br><br>set.seed(1)<br>print(fun1())<br><br>set.seed(1)<br>
print(fun2())<br><br>set.seed(1)<br>print(fun3())<br><br>set.seed(1)<br>print(fun4())<br><br>benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications = 1e6L)<br><br>output;<br><br>> set.seed(1)<br>
> print(fun1())<br>[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078<br><br>> set.seed(1)<br>> print(fun2())<br>[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078<br><br>> set.seed(1)<br>> print(fun3())<br>
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078<br><br>> set.seed(1)<br>> print(fun4())<br>[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078<br><br>> benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications = 1e6L)<br>
test replications elapsed relative user.self sys.self user.child sys.child<br>3 fun3() 1000000 7.15 1.000 7.01 0.00 NA NA<br>4 fun4() 1000000 7.33 1.025 7.27 0.01 NA NA<br>
2 fun2() 1000000 7.53 1.053 7.41 0.00 NA NA<br>1 fun1() 1000000 7.99 1.117 7.91 0.00 NA NA<br><br><br><div class="gmail_quote">On Tue, Sep 25, 2012 at 8:19 PM, Goldfeld, Keith <span dir="ltr"><<a href="mailto:Keith.Goldfeld@nyumc.org" target="_blank">Keith.Goldfeld@nyumc.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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.<br>
<div class="im HOEnZb"><br>
-----Original Message-----<br>
From: <a href="mailto:dmbates@gmail.com">dmbates@gmail.com</a> [mailto:<a href="mailto:dmbates@gmail.com">dmbates@gmail.com</a>] On Behalf Of Douglas Bates<br>
Sent: Tuesday, September 25, 2012 3:13 PM<br>
To: Rodney Sparapani<br>
Cc: Goldfeld, Keith; <a href="mailto:rcpp-devel@lists.r-forge.r-project.org">rcpp-devel@lists.r-forge.r-project.org</a><br>
Subject: Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array<br>
<br>
</div><div class="HOEnZb"><div class="h5">On Tue, Sep 25, 2012 at 1:53 PM, Rodney Sparapani <<a href="mailto:rsparapa@mcw.edu">rsparapa@mcw.edu</a>> wrote:<br>
> On 09/25/2012 01:37 PM, Goldfeld, Keith wrote:<br>
><br>
>>> code<- 'Rcpp::RNGScope scope;<br>
>><br>
>><br>
>> + Rcpp::NumericVector rn(5);<br>
>><br>
>> + for (int i=0; i< 5; i++) {<br>
>><br>
>> + rn(i) = rnorm(1,0,1);<br>
>><br>
>> + }<br>
>><br>
>> + return Rcpp::wrap(rn);<br>
>><br>
>> + '<br>
>><br>
>>><br>
>><br>
>>> fun<- cxxfunction(body=code, plugin="Rcpp")<br>
><br>
><br>
> You just need an explicit type conversion like so:<br>
><br>
> funk2 <- cxxfunction(body='RNGScope scope;<br>
> NumericVector rn(5);<br>
> for (int i=0; i < 5; i++) rn[i] = as<double>(rnorm(1,0,1));<br>
> return wrap(rn);',<br>
> includes="using namespace Rcpp;\n", plugin="Rcpp")<br>
><br>
> funk2()<br>
<br>
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.<br>
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.<br>
<br>
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.<br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br>