[Rcpp-devel] pb rgamma?

[Serguei Sokol]

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