[Rcpp-devel] pb rgamma?

[Serguei Sokol]

Le 10/04/2015 12:41, Sean O'Riordain a écrit :
> 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
Thanks Sean. However, I am a little bit puzzled by
a decision to implement a different algorithm in Rcpp
than in R. Isn't it contradictory with the spirit of sugar?

Serguei.

>
> 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
>
>     _________________________________________________
>     Rcpp-devel mailing list
>     Rcpp-devel at lists.r-forge.r-__project.org <mailto:Rcpp-devel at lists.r-forge.r-project.org>
>     https://lists.r-forge.r-__project.org/cgi-bin/mailman/__listinfo/rcpp-devel <https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel>
>
>
>
>
> --
> Kind regards,
>
> Sean O'Riordain
> seoriord at tcd.ie <mailto:seoriord at tcd.ie>
>
>