<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:Garamond;
panose-1:2 2 4 4 3 3 1 1 8 3;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Balloon Text Char";
margin:0in;
margin-bottom:.0001pt;
font-size:8.0pt;
font-family:"Tahoma","sans-serif";}
span.EmailStyle17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.BalloonTextChar
{mso-style-name:"Balloon Text Char";
mso-style-priority:99;
mso-style-link:"Balloon Text";
font-family:"Tahoma","sans-serif";}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<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.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D"><o:p> </o:p></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:jpolmsted@gmail.com]
<br>
<b>Sent:</b> Tuesday, September 25, 2012 4:29 PM<br>
<b>To:</b> Jeffrey Pollock<br>
<b>Cc:</b> Goldfeld, Keith; rcpp-devel@lists.r-forge.r-project.org<br>
<b>Subject:</b> Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Jeff and List,<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></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/">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.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></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.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></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,<o:p></o:p></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><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:7.5pt;font-family:"Arial","sans-serif";color:#222222;background:white">becomes</span><o:p></o:p></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><o:p></o:p></p>
</div>
<p class="MsoNormal">I haven't ever timed this, though.<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">JPO<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<div>
<p class="MsoNormal"><o:p> </o:p></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>
<o:p></o:p></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:<o:p></o:p></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>
<o:p></o:p></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:<o:p></o:p></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.<o:p></o:p></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<o:p></o:p></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><o:p></o:p></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">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><o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</div>
</div>
</body>
</html>