Passing seed/Setting seed/ C within R code

Question

I am making C within R codes.

In my C code I am using rand() function to generate random number. The R-ext.pdf says that I have to set a seed using the commands;

  GetRNGstate();
  PutRNGstate();

Although I am using these commands above, I am still getting different values for the same seed. Could you give me any help?

The minimum example is:

In C:

# include <R.h>
# include <Rinternals.h>
# include <Rmath.h>
# include <R_ext/Linpack.h> 

 SEXP example(){

   SEXP output;
   GetRNGstate();
   PROTECT(output = allocVector(INTSXP, 1));
   INTEGER(output)[0] = rand() % 50;
   PutRNGstate();
   UNPROTECT(1);
   return(output);
 }

In R:

dyn.load("example.so")
## The following codes return different values at ever run 
set.seed(1)
.Call("example")

Thanks in advance.

Answer 1

That is a logic error in your thinking -- you correctly set the seed, initialize the R RNG from code ... but then call the system RNG instead of the R RNG.

Replace rand() with unif_rand() (or norm_rand()) and you should be set.

Rcpp makes all this easier and gives you vectorised access to draws from the various distribution function (but you can of course do all this by hand too in C if you prefer).

By using cppFunction() from Rcpp, we now also take care of RNGScope which in turn provides GetRNGstate() / PutRNGstate() (while older examples still show instantiation of RNGScope; adding it does no harm as it does an equivalent of reference counting).

So it really is a one-liner to define, auto-extend, compile, and load this:

R> cppFunction("double myrand() { return norm_rand(); }")
R> for (i in 1:5) { set.seed(42); cat(i, " -- ", myrand(), "\n") }
1  --  1.37096 
2  --  1.37096 
3  --  1.37096 
4  --  1.37096 
5  --  1.37096 
R> 

whereas without the re-seeding we get

R> for (i in 1:5) { cat(i, " -- ", myrand(), "\n") }
1  --  -0.564698 
2  --  0.363128 
3  --  0.632863 
4  --  0.404268 
5  --  -0.106125 
R> 

Lastly, if you really want you can of course continue to use rand() (but see the literature on its sucky performance) but then use its seeding function instead rather than R's.