[Rcpp-devel] pb rgamma?

[Sean O'Riordain]

The R definition of the gamma distribution is the "inverse-rate", refer
https://stat.ethz.ch/R-manual/R-patched/library/stats/html/GammaDist.html
Could the Rcpp gamma distribution be using the standard "rate"
parameterization as seen on the Wikipedia page
http://en.wikipedia.org/wiki/Gamma_distribution? Nb. the "inverse-rate"
parameterization is simply a substitution of newrate=1/rate...

> set.seed(7)
> rgamma(1,3,1/4)
[1] 29.69732

Sean


On 10 April 2015 at 10:58, Serguei Sokol <serguei.sokol at gmail.com> wrote:

> Hi,
>
> I have came across a case when rgamma of rcpp does not produce
> the same result as rgamma of R, especially when the parameter rate
> is different from 1. Is it a bug of one or another?
>
> Minimal example follows.
>
> file rgamma_case.cpp:
> 8<-------
> #include <Rcpp.h>
> using namespace Rcpp;
> // [[Rcpp::export]]
> NumericVector rgamma_case(double rate) {
>    RNGScope scope;
>    NumericVector x(1);
>    x=rgamma(1, 3, rate);
>    return(x);
> }
> /*** R
> # first, same result
> set.seed(7)
> rate=1
> (x <- rgamma_case(rate))
> set.seed(7)
> (r <- rgamma(1, 3, rate))
> stopifnot(x==r)
>
> # now, different result
> set.seed(7)
> rate=4
> (x <- rgamma_case(rate))
> set.seed(7)
> (r <- rgamma(1, 3, rate))
> stopifnot(x==r)
> */
>
> 8<-------
>
> On my system, I've got
> > sourceCpp("rgamma_case.cpp")
>
> > # first, same result
> > set.seed(7)
>
> > rate=1
>
> > (x <- rgamma_case(rate))
> [1] 7.42433
>
> > set.seed(7)
>
> > (r <- rgamma(1, 3, rate))
> [1] 7.42433
>
> > stopifnot(x==r)
>
> > # now, different result
> > set.seed(7)
>
> > rate=4
>
> > (x <- rgamma_case(rate))
> [1] 29.69732
>
> > set.seed(7)
>
> > (r <- rgamma(1, 3, rate))
> [1] 1.856083
>
> > stopifnot(x==r)
> Erreur : x == r n'est pas TRUE
>
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>



-- 
Kind regards,

Sean O'Riordain
seoriord at tcd.ie
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150410/62e81b4a/attachment.html>