18

I'd like to generate identical random numbers in R and Julia. Both languages appear to use the Mersenne-Twister library by default, however in Julia 1.0.0:

julia> using Random
julia> Random.seed!(3)
julia> rand()
0.8116984049958615

Produces 0.811..., while in R:

set.seed(3)
runif(1)

produces 0.168.

Any ideas?

Related SO questions here and here.

My use case for those who are interested: Testing new Julia code that requires random number generation (e.g. statistical bootstrapping) by comparing output to that from equivalent libraries in R.

Julia Learner
  • 2,754
  • 15
  • 35
Colin T Bowers
  • 18,106
  • 8
  • 61
  • 89
  • 2
    A crude way would be to generate all the bootstrap replicates (or perhaps just the indices) up front and store them in a file that both programs could use. – joran Apr 07 '15 at 02:08
  • 1
    This isn't an answer, but I'm guessing the way the seed is turned into the initial state for the MT library isn't the same. I assume the answers can, and must, be found in the source (yay for open source). – IainDunning Apr 07 '15 at 02:45
  • @joran Agreed, and this is what I may end up doing. There is a bit of work to this though (for me at least - I'm a relative novice in R) as it implies altering both the R and Julia source to look for random numbers in the file. – Colin T Bowers Apr 07 '15 at 03:35
  • @IainDunning Sounds reasonable to me. I thought I'd ask here first just in case someone can answer in 5 minutes what could take me a full day :-) – Colin T Bowers Apr 07 '15 at 03:37
  • Using `RCall` doesnt help? – Khashaa Apr 07 '15 at 04:20
  • @Khashaa `RCall` certainly helps me transmit data between the two (say, for example, if I generate a random vector of numbers that I want to use as the source of randomness in both languages), but, as with the suggestion of joran, it still implies I'll need to edit the R source code of interest to point to that random vector. Admittedly, this is exactly what I'll probably end up doing :-) – Colin T Bowers Apr 07 '15 at 07:17
  • FYI, `srand(3); rand()` on my 32-bit Linux Julia platform produces `0.8116984049958615`. – rickhg12hs Apr 08 '15 at 19:52
  • @rickhg12hs I'm on Ubuntu 14.04 64-bit. Probably should have mentioned that in question. But I guess that difference just emphasizes the point in Dirk's answer - that there is no magic short-cut here... – Colin T Bowers Apr 10 '15 at 00:43

3 Answers3

8

That is an old problem.

Paul Gilbert addressed the same issue in the late 1990s (!!) when trying to assert that simulations in R (then then newcomer) gave the same result as those in S-Plus (then the incumbent).

His solution, and still the golden approach AFAICT: re-implement in fresh code in both languages as the this the only way to ensure identical seeding, state, ... and whatever else affects it.

Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725
  • In other words, no magic shortcut :-) Thanks for answering, I'll hold off on giving the tick for a day or two in case someone can come up with something that works in my particular special case. – Colin T Bowers Apr 07 '15 at 03:48
5

Pursuing the RCall suggestion made by @Khashaa, it's clear that you can set the seed and get the random numbers from R.

julia> using RCall

julia> RCall.reval("set.seed(3)")
RCall.NilSxp(16777344,Ptr{Void} @0x0a4b6330)

julia> a = zeros(Float64,20);

julia> unsafe_copy!(pointer(a), RCall.reval("runif(20)").pv, 20)
Ptr{Float64} @0x972f4860

julia> map(x -> @printf("%20.15f\n", x), a);
   0.168041526339948
   0.807516399072483
   0.384942351374775
   0.327734317164868
   0.602100674761459
   0.604394054040313
   0.124633444240317
   0.294600924244151
   0.577609919011593
   0.630979274399579
   0.512015897547826
   0.505023914156482
   0.534035353455693
   0.557249435689300
   0.867919487645850
   0.829708693316206
   0.111449153395370
   0.703688358888030
   0.897488264366984
   0.279732553754002

and from R:

> options(digits=15)
> set.seed(3)
> runif(20)
 [1] 0.168041526339948 0.807516399072483 0.384942351374775 0.327734317164868
 [5] 0.602100674761459 0.604394054040313 0.124633444240317 0.294600924244151
 [9] 0.577609919011593 0.630979274399579 0.512015897547826 0.505023914156482
[13] 0.534035353455693 0.557249435689300 0.867919487645850 0.829708693316206
[17] 0.111449153395370 0.703688358888030 0.897488264366984 0.279732553754002

** EDIT **

Per the suggestion by @ColinTBowers, here's a simpler/cleaner way to access R random numbers from Julia.

julia> using RCall

julia> reval("set.seed(3)");

julia> a = rcopy("runif(20)");

julia> map(x -> @printf("%20.15f\n", x), a);
   0.168041526339948
   0.807516399072483
   0.384942351374775
   0.327734317164868
   0.602100674761459
   0.604394054040313
   0.124633444240317
   0.294600924244151
   0.577609919011593
   0.630979274399579
   0.512015897547826
   0.505023914156482
   0.534035353455693
   0.557249435689300
   0.867919487645850
   0.829708693316206
   0.111449153395370
   0.703688358888030
   0.897488264366984
   0.279732553754002
rickhg12hs
  • 10,638
  • 6
  • 24
  • 42
  • 1
    Nice. So there is a shortcut via `RCall` which could probably be wrapped. But it just underlines the main point: if you want the same RNG stream, you *really* want the same code to run. I might start from a simple generator (say Marsaglia's KISS) and just code those up on both sides – Dirk Eddelbuettel Apr 12 '15 at 14:31
  • @DirkEddelbuettel , I searched for open-source Mersenne-Twisters and found examples at [Makoto Matsumoto's website](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) (many versions of source code for download and the original paper that included source code), [R source code](https://github.com/wch/r-source/blob/07ac1e21db175ac877530a5d0105906911e56c18/src/main/RNG.c#L561), and [GSL](https://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html). They are all a little different. Fortunately Julia's C interface works well and R provides a shared library, etc. – rickhg12hs Apr 12 '15 at 14:53
  • MT is too complicated. Use something like the simple congruent linear generators (for creating uniforms; eg something like KISS) in eg Ziggurat (which creates Normals). Ziggurat is used in Julia, and I have a package RcppZiggurat for R. – Dirk Eddelbuettel Apr 12 '15 at 15:00
  • Also, Julia's default global RNG is the latest dSFMT library. `Base.Random.globalRNG() |> typeof` returns `Base.Random.MersenneTwister`. R's Mersenne Twister is just slightly modified from the original paper. I think [R's seeding](https://github.com/wch/r-source/blob/07ac1e21db175ac877530a5d0105906911e56c18/src/main/RNG.c#L278) is a bit whacked as well. – rickhg12hs Apr 12 '15 at 15:12
  • Thanks for responding. This is interesting, although am I wrong in thinking it is a bit more complicated than it needs to be? Couldn't you just use (from Julia) `reval("set.seed(3)")` followed by `x = rcopy("runif(20)")` to import the numbers into the local variable `x` in Julia? Perhaps I'm missing something? – Colin T Bowers Apr 13 '15 at 02:02
  • @ColinTBowers Quite possible. I'm an `RCall` noob and just reported the first thing I tried and it worked. If there is a simpler and/or more robust way to get `R` random numbers into `Julia`, I'm all for it! – rickhg12hs Apr 13 '15 at 07:38
2

See:

?set.seed

"Mersenne-Twister": From Matsumoto and Nishimura (1998). A twisted GFSR with period 2^19937 - 1 and equidistribution in 623 consecutive dimensions (over the whole period). The ‘seed’ is a 624-dimensional set of 32-bit integers plus a current position in that set.

And you might see if you can link to the same C code from both languages. If you want to see the list/vector, type:

.Random.seed
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Yes, "but" if it were that easy you also get the same results via the GSL's Mersenne Twister or other. Usuaully, small local difference in set creation, manipulation etc pp get in the way. I'd just write a simple routine... – Dirk Eddelbuettel Apr 07 '15 at 11:17
  • 4
    In case anyone wonders which one of us to believe,... believe Dirk. It'll probably save you a lot of time. – IRTFM Apr 07 '15 at 15:26