有没有办法通过R使用常见的随机数进行采样？
在很多情况下，您需要多次执行以下操作(例如，如果您想在许多不同的参数值处绘制蒙特卡罗估计值)。首先，您从正态分布中采样一万个变量，然后，取这些样本的某个函数的平均值，返回单个 float 。现在，如果我想更改一些参数，更改这两个函数中的任何一个，我将不得不一遍又一遍地重新执行这些步骤。
最简单的方法是使用诸如 rnorm() 之类的函数一遍又一遍地对新的抽奖进行采样。一种不太幼稚的方法是使用不同的函数来获取大量常见随机数。但是，如果我使用这种方法，由于 R 主要使用按值传递语义，因此这里可能仍然会进行大量复制。有哪些工具可以让我解决这个问题并避免在第二种情况下进行所有这些复制？
最佳答案
我认为您在这里问了两种类型的问题:
以编程方式，我们能否以绕过 R 默认值传递的方式保留大量随机数据？
从数学上讲，如果我们提取大量随机数据并从中进行零碎挑选，我们可以任意更改提取中使用的参数吗？
问题 1 的答案是"is":R 中可以实现引用传递语义，但需要更多工作。我见过和使用过的所有实现都是通过环境或非 R 原生对象(指向结构体的 C/C++ 指针等)完成的。下面是一个示例，它缓存大量随机“正常”数据并在每次调用时检查可用数据池:
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:"is"。快速检查表明，最后三个数字确实是前一个 block 中第二个 my_rnorm(3) 中的数字“100 加”。因此，只需将“正常”随机数移动 mu/sigma 即可。我们在执行此操作时仍然使用大量预先提取的随机数据缓存。
但这值得吗？这本身就是一个幼稚的测试/比较，欢迎提出建设性的建议。
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
(所有中位数均以纳秒为单位。)因此，虽然“更频繁地进行较小的拉取”(使用 rnorm)会从这种缓存中受益，这似乎很直观，但我无法解释为什么会这样在拉动 100 或更大之前不是很有帮助。
这可以扩展到其他发行版吗？几乎可以确定。 “统一”是直接的(类似于缩放和移位)，但其他一些可能需要更多的微积分才能正确完成。 (例如，如果没有更多的研究，“t”分布如何改变预先提取的数据的自由度，这一点并不明显……如果可能的话。虽然我确实在某些方面认为自己是一名统计学家，但我我还没准备好对此做出肯定/否定/也许的回答:-)
关于r - 使用 r 中的常见随机数进行采样(高效!)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/57260745/