Can I have multiple independent random number generators in R like I can in C++?

Question

I have many independent random processes (arrival processes say) that require me to generate random numbers. I want to use common random numbers for each of these processes that I can compare how different policies perform when controlling these policies.

I want Process A to be governed by Generator A (using seed A) I want Process B to be governed by Generator B (using seed B) ..

and so on.

Is this possible to implement in R. I can't find anyone who has done it. I have tried. Forgive me if this is a repeated question.

Thanks

Jak

Answer 1

This is something that I've occassionally wanted to do - and haven't yet come up with much better than the following kludge (which is only really useful if you're using just 1 or 2 different random distributions, as you have to write a function for each:

#Make a list of seeds - generalises to mkore than 2
seed <- list(NA,NA)
set.seed(1)
seed[[1]] <- .Random.seed
set.seed(2)
seed[[2]] <- .Random.seed

my_runif <- function(...,which.seed=1)
{
  .Random.seed <<- seed[[which.seed]]
  x <-runif(...)
  seed[[which.seed]] <<- .Random.seed
  x
}

##Print some data for comparison
> set.seed(1); runif(10)
 [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968 0.94467527 0.66079779 0.629114040.06178627
> set.seed(2); runif(10)
 [1] 0.1848823 0.7023740 0.5733263 0.1680519 0.9438393 0.9434750 0.1291590 0.8334488 0.4680185 0.5499837

#Test
> my_runif(1,which.seed=1)
[1] 0.2655087
> my_runif(1,which.seed=1)
[1] 0.3721239
> my_runif(1,which.seed=1)
[1] 0.5728534
> my_runif(1,which.seed=2)
[1] 0.1848823
> my_runif(1,which.seed=1)
[1] 0.9082078

I'd imagine that the <<- will break if you call my_runif from inside another function.

fortunes::fortune("<<-")

ETA: The following might be more robust

my_runif <- function(...,which.seed=1)
{
  assign(".Random.seed", seed[[which.seed]], envir = .GlobalEnv)
  x <-runif(...)
  seed <- seed #Bring into local envir
  seed[[which.seed]] <- .Random.seed
  assign("seed", seed, envir = .GlobalEnv)
  x
}

Answer 2

Well the good news is that you already do -- see help(RNGkind):

 The currently available RNG kinds are given below.  ‘kind’ is
 partially matched to this list.  The default is
 ‘"Mersenne-Twister"’.

 ‘"Wichmann-Hill"’ [...]

 ‘"Marsaglia-Multicarry"’: [...]

 ‘"Super-Duper"’: [...]

 ‘"Mersenne-Twister"’: [...]

 ‘"Knuth-TAOCP-2002"’: [...]

 ‘"Knuth-TAOCP"’: [...]

 ‘"L'Ecuyer-CMRG"’: 

 ‘"user-supplied"’: Use a user-supplied generator.  See
      ‘Random.user’ for details.

and user-supplied lets you use your own.

And for N(0,1), you also have

 ‘normal.kind’ can be ‘"Kinderman-Ramage"’, ‘"Buggy
 Kinderman-Ramage"’ (not for ‘set.seed’), ‘"Ahrens-Dieter"’,
 ‘"Box-Muller"’, ‘"Inversion"’ (the default), or ‘"user-supplied"’.
 (For inversion, see the reference in ‘qnorm’.)  [...]

For parallel work, see the (excellent) vignette of the parallel package that came with R. There are existing generators for multiple threads/cores/... etc.

Last but not least, R is of course extensible and you could for example use Rcpp where we have a few posts on random numbers over at the Rcpp Gallery site.