<div dir="ltr"><div>The R definition of the gamma distribution is the "inverse-rate", refer <a href="https://stat.ethz.ch/R-manual/R-patched/library/stats/html/GammaDist.html">https://stat.ethz.ch/R-manual/R-patched/library/stats/html/GammaDist.html</a><br>Could the Rcpp gamma distribution be using the standard "rate" parameterization as seen on the Wikipedia page <a href="http://en.wikipedia.org/wiki/Gamma_distribution">http://en.wikipedia.org/wiki/Gamma_distribution</a>? Nb. the "inverse-rate" parameterization is simply a substitution of newrate=1/rate...<br><br>> set.seed(7)<br>> rgamma(1,3,1/4)<br>[1] 29.69732<br><br></div>Sean<br><br></div><div class="gmail_extra"><br><div class="gmail_quote">On 10 April 2015 at 10:58, Serguei Sokol <span dir="ltr"><<a href="mailto:serguei.sokol@gmail.com" target="_blank">serguei.sokol@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi,<br>
<br>
I have came across a case when rgamma of rcpp does not produce<br>
the same result as rgamma of R, especially when the parameter rate<br>
is different from 1. Is it a bug of one or another?<br>
<br>
Minimal example follows.<br>
<br>
file rgamma_case.cpp:<br>
8<-------<br>
#include <Rcpp.h><br>
using namespace Rcpp;<br>
// [[Rcpp::export]]<br>
NumericVector rgamma_case(double rate) {<br>
   RNGScope scope;<br>
   NumericVector x(1);<br>
   x=rgamma(1, 3, rate);<br>
   return(x);<br>
}<br>
/*** R<br>
# first, same result<br>
set.seed(7)<br>
rate=1<br>
(x <- rgamma_case(rate))<br>
set.seed(7)<br>
(r <- rgamma(1, 3, rate))<br>
stopifnot(x==r)<br>
<br>
# now, different result<br>
set.seed(7)<br>
rate=4<br>
(x <- rgamma_case(rate))<br>
set.seed(7)<br>
(r <- rgamma(1, 3, rate))<br>
stopifnot(x==r)<br>
*/<br>
<br>
8<-------<br>
<br>
On my system, I've got<br>
> sourceCpp("rgamma_case.cpp")<br>
<br>
> # first, same result<br>
> set.seed(7)<br>
<br>
> rate=1<br>
<br>
> (x <- rgamma_case(rate))<br>
[1] 7.42433<br>
<br>
> set.seed(7)<br>
<br>
> (r <- rgamma(1, 3, rate))<br>
[1] 7.42433<br>
<br>
> stopifnot(x==r)<br>
<br>
> # now, different result<br>
> set.seed(7)<br>
<br>
> rate=4<br>
<br>
> (x <- rgamma_case(rate))<br>
[1] 29.69732<br>
<br>
> set.seed(7)<br>
<br>
> (r <- rgamma(1, 3, rate))<br>
[1] 1.856083<br>
<br>
> stopifnot(x==r)<br>
Erreur : x == r n'est pas TRUE<br>
<br>
______________________________<u></u>_________________<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-<u></u>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-<u></u>project.org/cgi-bin/mailman/<u></u>listinfo/rcpp-devel</a><br>
</blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature">Kind regards,<div><br></div><div>Sean O'Riordain</div><div><a href="mailto:seoriord@tcd.ie" target="_blank">seoriord@tcd.ie</a></div><div><br></div><div><br></div></div>
</div>