How can I replicate random draws in RcppArmadillo?

Question

Here's a C++ function to draw N independent normal deviates with mean zero and standard deviation s

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::export]] 
List rnorm_cpp(double s, int N){

    arma::colvec epsilon = s * arma::randn(N);
    return List::create(Named("e") = epsilon);

}

And here's a (virtually identical) R version

rnormR <- function(s, N){

  epsilon <- rnorm(N, mean = 0, sd = s)
  return(list(e = epsilon))

}    

After sourcing rnorm_cpp and rnormR I ran the following:

set.seed(1234)
fooR <- rnormR(s = 5, N = 10)

set.seed(1234)
barR <- rnormR(s = 5, N = 10)

set.seed(1234)
fooCpp <- rnorm_cpp(s = 5, N = 10)

set.seed(1234)
barCpp <- rnorm_cpp(s = 5, N = 10)

Finally, I ran identical and got the following results:

> identical(fooR, barR)
[1] TRUE
> identical(barR, fooCpp)
[1] FALSE
> identical(fooCpp, barCpp)
[1] FALSE

I was expecting to get TRUE for all three of these. How can I: (1) replicate random draws across calls to rnorm_cpp and (2) get identical draws for calls to rnormR and rnorm_cpp?

Answer 1

The function arma::randn() is not connected to the R RNGs, so calling set.seed() has no influence on it.

What we do in Rcpp is to take advantage of the fine R API which permits us to access the same RNGs from R and C++. And by being careful with RNGScope instances (which get inserted automagically) the RNG state is always correct between R and C++.

But you simply cannot assume that any other third-party RNG (here: Arma's) was automatically aligned as well. Moreover, in this particular case, Conrad's documentation for Armadillo is clear:

To change the seed, use the std::srand() function

To clarify (Hi, @DWin) here is full R and C++ example:

R> set.seed(42); rnorm(5)           ## Five N(0,1) draws in R
[1]  1.3710 -0.5647  0.3631  0.6329  0.4043
R> cppFunction('NumericVector foo(int n) { return rnorm(n); }')
R> set.seed(42); foo(5)             ## Five N(0,1) draws from C++ fun.
[1]  1.3710 -0.5647  0.3631  0.6329  0.4043
R> 

We get the same numbers via R and C++ as we a) seed the RNGs identically and b) do in fact call the same RNG provided by R.