在R和Julia中生成相同的随机数

Question

我希望在R和Julia中生成相同的随机数。两种语言似乎都默认使用Mersenne-Twister库，但是在Julia 1.0.0中：

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

在 R 中会产生 0.811...。

set.seed(3)
runif(1)

产生0.168。

有什么想法吗？

相关的 SO 问题在这里和这里。

我要使用它的人可能会感兴趣：通过将输出与R中等效库的输出进行比较，测试需要随机数生成（例如统计Bootstrap）的新Julia代码。

- Colin T Bowers

---

一种简单的方法是预先生成所有引导复制品（或仅生成索引），并将它们存储在一个文件中，两个程序都可以使用该文件。 - joran

这不是一个答案，但我猜测种子被转换为MT库的初始状态的方式不同。我认为答案可以在源代码中找到（开源万岁）。 - IainDunning

@joran 同意，这可能是我最终要做的事情。不过这需要一些工作（至少对我来说 - 我在 R 中是相对新手），因为它意味着修改 R 和 Julia 源代码以查找文件中的随机数。 - Colin T Bowers

我的32位Linux Julia平台上，srand(3); rand()生成了0.8116984049958615。 - rickhg12hs

@rickhg12hs 我使用的是Ubuntu 14.04 64位版本。我应该在问题中提到这一点。但我猜这种差异只是强调Dirk的答案中没有魔法捷径这一点…… - Colin T Bowers

Answer 1

这是一个老问题。

在上世纪90年代末，Paul Gilbert曾尝试证明在R（当时的新手）中进行的模拟能够得到与S-Plus（当时的老牌软件）相同的结果，解决了同样的问题。

他提出的方案至今仍然是最佳实践：在两种语言中都重新编写代码以确保种子、状态等完全相同。

- Dirk Eddelbuettel

---

换句话说，没有魔法快捷方式 :-) 谢谢回答，我会等一两天再决定采纳，以防有人能够提出适用于我的特殊情况的解决方案。 - Colin T Bowers

Answer 2

跟随@Khashaa的RCall建议，我们可以清晰地设置种子并从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

并且来自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

** 编辑 **

根据@ColinTBowers的建议，这是一种更简单/更清晰的方法，可在Julia中访问R个随机数。

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

---

不错。那么通过 RCall 可能可以快捷地实现包装。但这也强调了一个主要的观点：如果你想要相同的 RNG 流，你 确实 需要运行相同的代码。我可能会从一个简单的生成器（比如 Marsaglia 的 KISS）开始，在两侧编写相应的代码。 - Dirk Eddelbuettel

@DirkEddelbuettel，我搜索了开源的Mersenne-Twisters，并在Makoto Matsumoto的网站上找到了示例（包括多个版本的源代码下载和原始论文中包含的源代码），R源代码和GSL。它们都有一些不同。幸运的是，Julia的C接口很好用，而R提供了共享库等。 - rickhg12hs

MT太过复杂。可以使用像简单同余线性生成器（用于创建均匀分布，例如像KISS这样的东西）在Ziggurat中创建正态分布。Julia中使用了Ziggurat，并且我在R中有一个名为RcppZiggurat的包。 - Dirk Eddelbuettel

此外，Julia的默认全局RNG是最新的dSFMT库。Base.Random.globalRNG() |> typeof 返回 Base.Random.MersenneTwister。R的Mersenne Twister仅略微修改了原始论文。我认为R的种子也有点问题。 - rickhg12hs

谢谢您的回复。这很有趣，但我想它比必要的要复杂一些，我错了吗？你不能只使用（从Julia）reval("set.seed(3)")然后跟着x = rcopy("runif(20)")将数字导入到Julia的本地变量x中吗？也许我漏掉了什么？ - Colin T Bowers

@ColinTBowers 可能性很大。我是一个RCall新手，只是报告了我尝试过的第一件事情，并且它起作用了。如果有更简单和/或更可靠的方法将R随机数传递到Julia中，那我完全赞成！ - rickhg12hs

Answer 3

请参见：

?set.seed

"Mersenne-Twister"是Matsumoto和Nishimura（1998）发明的一种扭曲的GFSR，其周期为2^19937-1，在整个周期内具有623个连续维度的等分布性。 "种子"是一个由32位整数组成的624维集合以及该集合中的当前位置。

如果您想查看列表/向量，请尝试从两种语言中链接到相同的C代码。

.Random.seed

- IRTFM

---

是的，“但是”如果这么容易，您也可以通过GSL的Mersenne Twister或其他方式获得相同的结果。通常，集合创建、操作等方面的小局部差异会成为障碍。我只会编写一个简单的例程... - Dirk Eddelbuettel

如果有人想知道该相信我们中的谁，那就相信Dirk吧。这样做可能会为你节省很多时间。 - IRTFM