Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Generate identical random numbers in R and Julia

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.

like image 239
Colin T Bowers Avatar asked Apr 07 '15 01:04

Colin T Bowers


People also ask

How do I use random seed in Julia?

In Julia, you can set a seed to the random number generator using the srand() function. The code example below sets the seed to 1234. Generating a random variable with rand(1) after setting the seed to 1234 will always generate the same number, i.e. it will always return 0.5908446386657102.

How do you generate a random matrix in Julia?

Use rand(Uniform(-2,1),3,3) or rand(Normal(x,y), 3,3) to get normally distributed values.


2 Answers

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.

like image 189
Dirk Eddelbuettel Avatar answered Sep 19 '22 16:09

Dirk Eddelbuettel


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
like image 27
rickhg12hs Avatar answered Sep 19 '22 16:09

rickhg12hs