Indeed, changing the sample size from 5 to 500 results in;<br><br>> benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications = 1e5L)<br> test replications elapsed relative user.self sys.self user.child sys.child<br>
4 fun4() 100000 6.16 1.000 6.06 0 NA NA<br>3 fun3() 100000 6.19 1.005 5.91 0 NA NA<br>2 fun2() 100000 6.92 1.123 6.75 0 NA NA<br>
1 fun1() 100000 13.01 2.112 12.73 0 NA NA<br><br><div class="gmail_quote">On Tue, Sep 25, 2012 at 9:33 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">
<div link="blue" vlink="purple" lang="EN-US">
<div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">From:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> Jonathan Olmsted [mailto:<a href="mailto:jpolmsted@gmail.com" target="_blank">jpolmsted@gmail.com</a>]
<br>
<b>Sent:</b> Tuesday, September 25, 2012 4:29 PM<br>
<b>To:</b> Jeffrey Pollock</span></p><div><div class="h5"><br>
<b>Cc:</b> 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>
<b>Subject:</b> Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array<u></u><u></u></div></div><p></p><div><div class="h5">
<p class="MsoNormal"><u></u> <u></u></p>
<p class="MsoNormal">Jeff and List,<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">I was interested in the same fairly recently. The full shebang is here <<a href="http://www.rochester.edu/college/gradstudents/jolmsted/blog/2012/08/22/rcpp-rng-performance/" target="_blank">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">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.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">For what it is worth, instead of an explicit "as<>()" conversion, I usually pull out the first element of the returned NumericVector. That is,<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:7.5pt;font-family:"Courier New";color:#222222;background:white">rn[i] = as<double>(rnorm(1, 0.0, 1.0));</span><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:7.5pt;font-family:"Arial","sans-serif";color:#222222;background:white">becomes</span><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:7.5pt;font-family:"Courier New";color:#222222;background:white">rn[i] = rnorm(1, 0.0, 1.0)[0];</span><u></u><u></u></p>
</div>
<p class="MsoNormal">I haven't ever timed this, though.<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">JPO<u></u><u></u></p>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<p class="MsoNormal" style="margin-bottom:12.0pt"><span style="font-family:"Arial","sans-serif";color:#666666">-------------------------------------------------------------------------</span><span style="color:#666666"><br>
</span><span style="font-family:"Garamond","serif";color:#666666">J. P. Olmsted<br>
<a href="mailto:j.p.olmsted@gmail.com" target="_blank">j.p.olmsted@gmail.com</a><br>
</span><span style="color:#666666"><a href="http://about.me/olmjo" target="_blank"><span style="font-family:"Garamond","serif"">http://about.me/olmjo</span></a><br>
</span><span style="font-size:10.0pt;font-family:"Arial","sans-serif";color:#666666">-------------------------------------------------------------------------</span><br>
<span style="font-size:10.0pt;font-family:"Arial","sans-serif""><br>
</span><br>
<br>
<br>
<u></u><u></u></p>
<div>
<p class="MsoNormal">On Tue, Sep 25, 2012 at 4:08 PM, Jeffrey Pollock <<a href="mailto:jeffpollock9@gmail.com" target="_blank">jeffpollock9@gmail.com</a>> wrote:<u></u><u></u></p>
<p class="MsoNormal" style="margin-bottom:12.0pt">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>
<u></u><u></u></p>
<div>
<p class="MsoNormal">On Tue, Sep 25, 2012 at 8:19 PM, Goldfeld, Keith <<a href="mailto:Keith.Goldfeld@nyumc.org" target="_blank">Keith.Goldfeld@nyumc.org</a>> wrote:<u></u><u></u></p>
<p class="MsoNormal">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.<u></u><u></u></p>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt"><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<u></u><u></u></p>
</div>
<div>
<div>
<p class="MsoNormal">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><u></u><u></u></p>
</div>
</div>
</div>
<p class="MsoNormal"><br>
<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><u></u><u></u></p>
</div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
</div>
</div></div></div>
</div>
</blockquote></div><br>