使用 r 中的公共随机数进行采样（高效！）

Question

有什么方法可以用R使用普通随机数进行抽样吗？

在很多情况下，您需要多次执行以下操作（例如，如果您想绘制 Monte Carlo 多个不同参数值的估计值）。首先，你从正态分布中抽取一万个变量，其次，你取这些样本的某个函数的平均值，返回单个浮点数。现在，如果我想更改一些参数，更改这两个函数中的任何一个，我将不得不一遍又一遍地重新执行这些步骤。

天真的方法是使用像 rnorm() 这样的函数一遍又一遍地对新抽取样本进行采样。一个不太天真的方法是使用一个不同的函数，它需要大量的普通随机数。但是，如果我使用这种方法，由于 R 主要使用按值传递语义，这里可能仍会进行大量复制。有哪些工具可以让我绕过这个问题并避免在第二种情况下进行所有这些复制？

Answer 1

我认为你在这里问两种类型的问题：

在编程上，我们能否以避开 R 的默认传递值的方式保留大量随机数据？
从数学上讲，如果我们拉取大量随机数据并从中零碎地挑选，我们可以任意更改拉取中使用的参数吗？

1 的答案是 "yes"：在 R 中可以使用按引用传递语义，但需要做更多的工作。我见过和玩过的所有实现都是使用 environments 或非 R 本地对象（C/C++ 指向结构的指针等）完成的。这是一个缓存大量随机 "normal" 数据并在每次调用时检查可用数据池的示例：

my_rnorm_builder <- function(deflen = 10000) {
  .cache <- numeric(0)
  .index <- 0L
  .deflen <- deflen
  check <- function(n) {
    if ((.index + n) > length(.cache)) {
      message("reloading") # this should not be here "in-production"
      l <- length(.cache)
      .cache <<- c(.cache[ .index + seq_len(l - .index) ],
                   rnorm(.deflen + n + l))
      .index <<- 0L
    }
  }
  function(n, mean = 0, sd = 1) {
    check(n)
    if (n > 0) {
      out <- mean + sd * .cache[ .index + seq_len(n) ]              
      .index <<- .index + as.integer(n)
      return(out)
    } else return(numeric(0))
  }
}

到目前为止，它对敌对用户或其他可能的错误没有弹性。它不保证可用剩余随机数的长度。（考虑到基准，进行这样的检查会使它的速度降低到合理的阈值以下。）

运行演示：

my_rnorm <- my_rnorm_builder(1e6)
# starts empty
get(".index", env=environment(my_rnorm))
# [1] 0
length(get(".cache", env=environment(my_rnorm)))
# [1] 0

set.seed(2)
my_rnorm(3) # should see "reloading"
# reloading
# [1] -0.8969145  0.1848492  1.5878453
my_rnorm(3) # should not see "reloading"
# [1] -1.13037567 -0.08025176  0.13242028
# prove that we've changed things internally
get(".index", env=environment(my_rnorm))
# [1] 6
length(get(".cache", env=environment(my_rnorm)))
# [1] 1000003

head(my_rnorm(1e6)) # should see "reloading"
# reloading
# [1]  0.7079547 -0.2396980  1.9844739 -0.1387870  0.4176508  0.9817528

让我们通过重新开始并重新设置我们的种子来确保 sigma*x+mu 的随机数缩放有意义：

# reload the definition of my_rnorm
my_rnorm <- my_rnorm_builder(1e6)
length(get(".cache", env=environment(my_rnorm)))
# [1] 0
set.seed(2)
my_rnorm(3) # should see "reloading"
# reloading
# [1] -0.8969145  0.1848492  1.5878453
my_rnorm(3, mean = 100) # should not see "reloading"
# [1]  98.86962  99.91975 100.13242

所以回答问题 2："yes"。快速检查发现最后三个数字确实是前一个区块中第二个 my_rnorm(3) 中的数字的“100 加”。因此，只需将 "normal" 个随机数移动 mu/sigma 即可。我们这样做的同时仍然使用随机数据的大型预拉缓存。

但值得吗？这本身就是幼稚test/comparison，欢迎提出建设性意见。

t(sapply(c(1,5,10,100,1000,10000), function(n) {
  s <- summary(microbenchmark::microbenchmark(
    base = rnorm(n),
    my   = my_rnorm(n),
    times = 10000, unit = "ns"
  ))
  c(n = n, setNames(s$median, s$expr))  
}))
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
# reloading
#          n   base    my
# [1,]     1   1100  1100
# [2,]     5   1400  1300
# [3,]    10   1600  1400
# [4,]   100   6400  2000
# [5,]  1000  53100  6600
# [6,] 10000 517000 49900

（所有中位数都以纳秒为单位。）因此，虽然 "smaller pulls done more frequently"（使用 rnorm）会从这种缓存中受益似乎很直观，但我无法解释为什么它不是很有用直到拉动 100 或更大。

这可以扩展到其他发行版吗？几乎可以确定。 "Uniform" 将是直截了当的（类似缩放和移位），但其他一些可能需要更多的微积分才能正确完成。（例如，如果不进行更多研究，"t" 分布如何改变预拉数据的自由度并不明显......如果这可能的话。虽然我在某些方面确实认为自己是统计学家, 我还没准备好在那个上申请 yes/no/maybe :-)

Answer 2

r2evans 关于是否值得的回答的补充？：我不这么认为，因为除了缓存随机抽取之外，还可以使用更快的 RNG。在这里，我将 dqrng 包中的 dqrnorm 添加到比较中：

dqrnorm 是 n <= 100
对于n > 100，缓存和dqrnorm具有可比性并且比rnorm

使用 r 中的公共随机数进行采样（高效！）

Sampling using common random numbers in r (efficiently!)

r

pass-by-reference

pass-by-value