<font>Jeff and List,</font><div><font><br></font></div><div><font>I was interested in the same fairly recently. The full shebang is here <</font><a href="http://www.rochester.edu/college/gradstudents/jolmsted/blog/2012/08/22/rcpp-rng-performance/">http://www.rochester.edu/college/gradstudents/jolmsted/blog/2012/08/22/rcpp-rng-performance/</a>>. I, unlike you, found some difference between the various approaches.</div>
<div><font><br></font></div><div><font>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.</font></div>
<div><font><br></font></div><div><font>For what it is worth, instead of an explicit "as<>()" conversion, I usually pull out the first element of the returned NumericVector. That is,</font></div><div><span style="color:rgb(34,34,34);font-size:10px;background-color:rgb(255,255,255)"><font face="courier new, monospace">rn[i] = as<double>(rnorm(1, 0.0, 1.0));</font></span><br style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:10px;background-color:rgb(255,255,255)">
</div><div><span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:10px;background-color:rgb(255,255,255)">becomes</span></div><div><span style="color:rgb(34,34,34);font-size:10px;background-color:rgb(255,255,255)"><font face="courier new, monospace">rn[i] = rnorm(1, 0.0, 1.0)[0];</font></span></div>
I haven't ever timed this, though.<div><br></div><div>JPO<br><div><span style="color:rgb(34,34,34);font-size:10px;background-color:rgb(255,255,255)"><font face="courier new, monospace"><br></font></span></div><div><span style="color:rgb(34,34,34);font-size:10px;background-color:rgb(255,255,255)"><font face="courier new, monospace"><br>
</font></span></div><div><div><br></div><font><font color="#666666"><font face="arial, helvetica, sans-serif">-------------------------------------------------------------------------</font><br><font face="garamond, serif">J. P. Olmsted<br>
<a href="mailto:j.p.olmsted@gmail.com" target="_blank">j.p.olmsted@gmail.com</a><br></font></font></font><font color="#666666"><a href="http://about.me/olmjo" target="_blank"><font face="garamond, serif">http://about.me/olmjo</font></a><font><br>
</font><font style="font-family:arial,helvetica,sans-serif" size="2"></font><font style="font-family:arial,helvetica,sans-serif" size="2">-------------------------------------------------------------------------</font></font><br>
<font style="font-family:arial,helvetica,sans-serif" size="2"><br></font><br><br>
<br><br><div class="gmail_quote">On Tue, Sep 25, 2012 at 4:08 PM, Jeffrey Pollock <span dir="ltr"><<a href="mailto:jeffpollock9@gmail.com" target="_blank">jeffpollock9@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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><br>
-----Original Message-----<br>
From: <a href="mailto:dmbates@gmail.com" target="_blank">dmbates@gmail.com</a> [mailto:<a href="mailto:dmbates@gmail.com" target="_blank">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" target="_blank">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><div>On Tue, Sep 25, 2012 at 1:53 PM, Rodney Sparapani <<a href="mailto:rsparapa@mcw.edu" target="_blank">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" target="_blank">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>
<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></blockquote></div><br></div></div>