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.

8
  • 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. Commented Apr 7, 2015 at 2: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). Commented Apr 7, 2015 at 2: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. Commented Apr 7, 2015 at 3: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 :-) Commented Apr 7, 2015 at 3:37
  • Using RCall doesnt help? Commented Apr 7, 2015 at 4:20

3 Answers 3

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.

Sign up to request clarification or add additional context in comments.

Comments

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

6 Comments

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
@DirkEddelbuettel , I searched for open-source Mersenne-Twisters and found examples at Makoto Matsumoto's website (many versions of source code for download and the original paper that included source code), R source code, and GSL. They are all a little different. Fortunately Julia's C interface works well and R provides a shared library, etc.
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.
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 is a bit whacked as well.
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?
|
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

2 Comments

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...
In case anyone wonders which one of us to believe,... believe Dirk. It'll probably save you a lot of time.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.