从全局R进程中隔离本地环境的随机性

Question

我们可以使用set.seed()在R中设置随机种子,这具有全局效应.这是一个简单的例子来说明我的目标:

set.seed(0)
runif(1)
# [1] 0.8966972

set.seed(0)
f <- function() {
  # I do not want this random number to be affected by the global seed
  runif(1)
}
f()
# [1] 0.8966972

基本上我希望能够避免全局随机种子(即.Random.seed)在本地环境(例如R函数)中的影响,这样我就可以实现用户无法控制的某种随机性.例如,即使用户有set.seed(),他每次调用此函数时仍会得到不同的输出.

现在有两个实现.第一个依赖于set.seed(NULL)让R每次想要得到一些随机数时重新初始化随机种子:

createUniqueId <- function(bytes) {
  withPrivateSeed(
    paste(as.hexmode(sample(256, bytes, replace = TRUE) - 1), collapse = "")
  )
}
withPrivateSeed <- function(expr, seed = NULL) {
  oldSeed <- if (exists('.Random.seed', envir = .GlobalEnv, inherits = FALSE)) {
    get('.Random.seed', envir = .GlobalEnv, inherits = FALSE)
  }
  if (!is.null(oldSeed)) {
    on.exit(assign('.Random.seed', oldSeed, envir = .GlobalEnv), add = TRUE)
  }
  set.seed(seed)
  expr
}

你可以看到我得到不同的id字符串,即使我将种子设置为0,并且全局随机数流仍然可以重现:

> set.seed(0)
> runif(3)
[1] 0.8966972 0.2655087 0.3721239
> createUniqueId(4)
[1] "83a18600"
> runif(3)
[1] 0.5728534 0.9082078 0.2016819

> set.seed(0)
> runif(3)  # same
[1] 0.8966972 0.2655087 0.3721239
> createUniqueId(4)  # different
[1] "77cb3d91"
> runif(3)
[1] 0.5728534 0.9082078 0.2016819

> set.seed(0)
> runif(3)
[1] 0.8966972 0.2655087 0.3721239
> createUniqueId(4)
[1] "c41d61d8"
> runif(3)
[1] 0.5728534 0.9082078 0.2016819

第二种方案可以发现这里 Github上.它更复杂,基本思想是:

• 使用set.seed(NULL)(in .onLoad())在程序包启动期间初始化随机种子
• 将随机种子存储在单独的环境中(.globals$ownSeed)
• 每当我们想要生成随机数时:
  1. 将本地种子分配给全局随机种子
  2. 生成随机数
  3. 将新的全局种子(由于步骤2而更改)分配给本地种子
  4. 将全局种子恢复为其原始值

现在我的问题是这两种方法在理论上是否相同.第一种方法的随机性依赖于createUniqueId()调用时的当前时间和进程ID ,第二种方法似乎依赖于加载包时的时间和进程ID.对于第一种方法,是否有可能createUniqueId()在同一个R进程中同时发生两次调用,以便它们返回相同的id字符串？

更新

在下面的答案中,Robert Krzyzanowski提供了一些set.seed(NULL)可能导致严重ID冲突的经验证据.我做了一个简单的可视化:

createGlobalUniqueId <- function(bytes) {
  paste(as.hexmode(sample(256, bytes, replace = TRUE) - 1), collapse = "")
}
n <- 10000
length(unique(replicate(n, createGlobalUniqueId(5))))
length(unique(x <- replicate(n, createUniqueId(5))))
# denote duplicated values by 1, and unique ones by 0
png('rng-time.png', width = 4000, height = 400)
par(mar = c(4, 4, .1, .1), xaxs = 'i')
plot(1:n, duplicated(x), type = 'l')
dev.off()

当线到达图的顶部时,这意味着生成了重复值.但是,请注意这些重复不是连续的,即any(x[-1] == x[-n])通常FALSE.可能存在与系统时间相关的重复模式.由于我对基于时间的随机种子的工作方式缺乏了解,我无法进一步调查,但是你可以在这里和这里看到C源代码的相关部分.

Answer 1

我认为在你的函数中只有一个独立的RNG会很好,它不受全局种子的影响,但会拥有自己的种子.结果,randtoolbox提供此功能:

library(randtoolbox)
replicate(3, {
  set.seed(1)
  c(runif(1), WELL(3), runif(1))
})   
#            [,1]      [,2]      [,3]
#[1,] 0.265508663 0.2655087 0.2655087
#[2,] 0.481195594 0.3999952 0.9474923
#[3,] 0.003865934 0.6596869 0.4684255
#[4,] 0.484556709 0.9923884 0.1153879
#[5,] 0.372123900 0.3721239 0.3721239

顶部和底部行受种子影响,而中间行则"真正随机".

基于此,这是您的功能的实现:

sample_WELL <- function(n, size=n) {
  findInterval(WELL(size), 0:n/n)
}

createUniqueId_WELL <- function(bytes) {
  paste(as.hexmode(sample_WELL(256, bytes) - 1), collapse = "")
}

length(unique(replicate(10000, createUniqueId_WELL(5))))
#[1] 10000

# independency on the seed: 
set.seed(1)
x <- replicate(100, createGlobalUniqueId(5))
x_WELL <- replicate(100, createUniqueId_WELL(5))
set.seed(1)
y <- replicate(100, createGlobalUniqueId(5))
y_WELL <- replicate(100, createUniqueId_WELL(5))
sum(x==y)
#[1] 100
sum(x_WELL==y_WELL)
#[1] 0

编辑

为了理解为什么我们得到重复的密钥,我们应该看一下调用时会发生什么set.seed(NULL).所有与RNG相关的代码都是用C语言编写的,所以我们应该直接转到svn.r-project.org/R/trunk/src/main/RNG.c并参考该函数do_setseed.如果seed = NULL那么明确TimeToSeed被称为.有一条评论声明它应该位于datetime.c中,但是,它可以在svn.r-project.org/R/trunk/src/main/times.c中找到.

导航R源可能很困难,所以我在这里粘贴函数:

/* For RNG.c, main.c, mkdtemp.c */
attribute_hidden
unsigned int TimeToSeed(void)
{
    unsigned int seed, pid = getpid();
#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_REALTIME)
    {
    struct timespec tp;
    clock_gettime(CLOCK_REALTIME, &tp);
    seed = (unsigned int)(((uint_least64_t) tp.tv_nsec << 16) ^ tp.tv_sec);
    }
#elif defined(HAVE_GETTIMEOFDAY)
    {
    struct timeval tv;
    gettimeofday (&tv, NULL);
    seed = (unsigned int)(((uint_least64_t) tv.tv_usec << 16) ^ tv.tv_sec);
    }
#else
    /* C89, so must work */
    seed = (Int32) time(NULL);
#endif
    seed ^= (pid <<16);
    return seed;
}

所以每次我们打电话时set.seed(NULL),R都会执行以下步骤:

1. 以秒和纳秒为单位获取当前时间(如果可能,以#if defined块为单位的平台依赖性)
2. 将位移应用于纳秒,并将位'xor'es应用于秒
3. 将位移应用于pid并将其与前一个结果进行"xor"运算
4. 将结果设置为新种子

好吧,现在几乎可以看出,当生成的种子发生碰撞时,我们会得到重复的值.我的猜测是当两个调用落在1秒内时会发生这种情况,因此tv_sec是常量.为了证实这一点,我引入了一个滞后:

createUniqueIdWithLag <- function(bytes, lag) {
  Sys.sleep(lag)
  createUniqueId(bytes)
}
lags <- 1 / 10 ^ (1:5)
sapply(lags, function(x) length(unique(replicate(n, createUniqueIdWithLag(5, x)))))
[1] 1000 1000  996  992  990

令人困惑的是,即使滞后与纳秒相比仍然很大,我们仍然会发生碰撞!让我们进一步深入研究,我为种子写了一个"调试模拟器":

emulate_seed <- function() {
  tv <- as.numeric(system('echo $(($(date +%s%N)))', intern = TRUE))
  pid <- Sys.getpid()
  tv_nsec <- tv %% 1e9
  tv_sec <- tv %/% 1e9
  seed <- bitwXor(bitwShiftL(tv_nsec, 16), tv_sec)
  seed <- bitwXor(bitwShiftL(pid, 16), seed)
  c(seed, tv_nsec, tv_sec, pid)
}

z <- replicate(1000, emulate_seed())
sapply(1:4, function(i) length(unique(z[i, ])))
# unique seeds, nanosecs, secs, pids:
#[1]  941 1000   36    1

这实在令人困惑:纳秒都是独一无二的,并不能保证最终种子的唯一性.为了证明这种效果,这里有一个重复:

#            [,1]        [,2] 
#[1,] -1654969360 -1654969360
#[2,]   135644672   962643456
#[3,]  1397894128  1397894128 
#[4,]        2057        2057
bitwShiftL(135644672, 16)
#[1] -973078528
bitwShiftL(962643456, 16)
#[1] -973078528

最后一点:这两个数字的二进制表示和移位是

00001000000101011100011000000000 << 16 => 1100011000000000 + 16 zeroes
00111001011000001100011000000000 << 16 => 1100011000000000 + 16 zeroes

所以,是的,这实际上是一次不必要的碰撞.

好吧,毕竟,最后的结论是:set.seed(NULL)容易受到高负荷的影响,并且在处理多个连续呼叫时不能保证没有冲突!