mercredi 29 juillet 2020

Is there a best practice for using non-R RNGs in Rcpp code?

I am writing an Rcpp function wrapped by an R function. The Rcpp function does some parallel random sampling, necessitating that I not use R's RNG. (Note: this is not a question about random seeds for parallel processes; I've got that handled.) I would like users to be able to call set.seed() to maintain reproducibility, even though I'm using a different RNG. A pseudocode example of what I am doing is below. (I can edit to make working code if needed.)

My questions are: (a) Is there something obviously dangerous to the below approach? (It feels hacky.) (b) Is there a more well-established method?

FWIW: there is a similar question whose answer describes but does not show what I am trying to do, seeding a call to an RNG with R's RNG. But I can't find an explicit example of someone doing this.

R wrapper function

sample_wrapper <- function() {
  # get current state of RNG stream
  # first entry of .Random.seed is an integer representing the algorithm used 
  # second entry is current position in RNG stream
  # subsequent entries are pseudorandom numbers
  seed_pos <- .Random.seed[2]

  seed <- .Random.seed[seed_pos + 2]

  out <- sample_cpp(seed = seed)

  # move R's position in the RNG stream forward by 1 with a throw away sample
  runif(1)

  # return the output
  out
}

Rcpp sampling function

sample_cpp(const int seed){
  // seed my own rng
  MyRngClass my_rng = MyRng(seed);

  // do my own sample
  double result = MySampleFun(rng = my_rng);

  // return result
  return result;
}

Example usage

set.seed(123)

sample_wrapper()



Aucun commentaire:

Enregistrer un commentaire