[Rcpp-devel] pb rgamma?

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 <mailto: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