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

【问题标题】：Generate identical random numbers in R and Julia在 R 和 Julia 中生成相同的随机数
【发布时间】：2015-06-11 12:11:44
【问题描述】：

我想在 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 问题 here 和 here。

我对感兴趣的人的用例：通过将输出与 R 中等效库的输出进行比较来测试需要随机数生成（例如统计引导）的新 Julia 代码。

【问题讨论】：

一种粗略的方法是预先生成所有引导程序副本（或者可能只是索引）并将它们存储在两个程序都可以使用的文件中。
这不是一个答案，但我猜种子变成 MT 库的初始状态的方式不一样。我认为答案可以而且必须在源代码中找到（对开源来说是的）。
@joran 同意，这就是我最终可能会做的事情。不过，这需要做一些工作（至少对我来说——我是 R 的相对新手），因为这意味着要同时更改 R 和 Julia 源以在文件中查找随机数。
@IainDunning 对我来说听起来很合理。我想我先在这里问一下，以防有人能在 5 分钟内回答我可能需要一整天的时间:-)
使用RCall 没有帮助？

标签： r random julia mersenne-twister

【解决方案1】：

这是一个老问题。

Paul Gilbert 在 1990 年代后期（！！）试图断言 R（当时的新人）中的模拟与 S-Plus（当时的现任）中的模拟结果相同时解决了同样的问题。

他的解决方案，仍然是 AFAICT 的黄金方法：用两种语言重新实现新代码，因为这是确保相同的种子、状态……以及任何其他影响它的唯一方法。

【讨论】：

换句话说，没有神奇的捷径 :-) 感谢您的回答，我会推迟一两天打勾，以防有人能想出对我特别有用的东西案例。

【解决方案2】：

见：

?set.seed

“梅森旋风”：来自松本和西村（1998 年）。一个扭曲的 GFSR，周期为 2^19937 - 1，在 623 个连续维度（在整个周期内）均等分布。 “种子”是一个 624 维的 32 位整数集合加上该集合中的当前位置。

您可能会看到是否可以链接到两种语言的相同 C 代码。如果您想查看列表/向量，请输入：

.Random.seed

【讨论】：

是的，“但是”如果这么简单，您也可以通过 GSL 的 Mersenne Twister 或其他设备获得相同的结果。通常，集合创建、操作等方面的微小局部差异会妨碍您。我只是写一个简单的例程......
如果有人想知道我们中的哪一个应该相信，……相信德克。它可能会为您节省很多时间。

【解决方案3】：

按照@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

【讨论】：

不错。所以有一个通过RCall 的快捷方式可能会被包装。但它只是强调了要点：如果您想要相同的 RNG 流，您真的想要运行相同的代码。我可能会从一个简单的生成器（比如 Marsaglia 的 KISS）开始，然后在两边编写代码
@DirkEddelbuettel ，我搜索了开源 Mersenne-Twisters 并在 Makoto Matsumoto's website 找到了示例（许多版本的源代码可供下载以及包含源代码的原始论文），R source code 和GSL。他们都有点不同。幸运的是 Julia 的 C 接口运行良好，R 提供了共享库等。
MT 太复杂了。在例如 Ziggurat（创建法线）中使用类似简单的全等线性生成器（用于创建制服；例如 KISS 之类的东西）。 Ziggurat 在 Julia 中使用，我有一个 RcppZiggurat 包。
另外，Julia 的默认全局 RNG 是最新的 dSFMT 库。 Base.Random.globalRNG() |> typeof 返回Base.Random.MersenneTwister。 R 的 Mersenne Twister 只是对原始论文稍作修改。我认为R's seeding 也有点过时了。
感谢您的回复。这很有趣，虽然我认为它比它需要的复杂一点是错误的吗？您不能只使用（来自 Julia）reval("set.seed(3)") 后跟 x = rcopy("runif(20)") 将数字导入 Julia 中的局部变量 x 吗？也许我错过了什么？