使用带有 NA 和 mcmapply 的 rbinom 进行优化

Question

我正在尝试使用两个数值向量从二项式分布中采样数据，每个向量都大约有 800 万个元素。它们看起来像这样：


b = runif(n = 8000000, min = 0, max = 1)
d=rpois( 8000000, lambda =0.1)

我想对 b 和 d 的每个元素应用这个函数

rbinom_NA = function(x , y) {
  result = rbinom(n = x,
                  size = 1,
                  prob = y) 
  
  if (length(result) == 0)
    return(NA)
  return(result)
}

我用mcmapply做的。你可以在我的电脑上看到我需要多长时间。

rbinom_result = vector(mode = "list", length = 8000000)
rbinom_result = mcmapply(
  d,
  b,
  FUN = rbinom_NA, 
  mc.cores = detectCores()-1) #run on Ubuntu, 32 cores and 64Gb memory
# system.time result:
#  user  system elapsed 
# 90.631 116.584 154.353

然后我计算之前结果的每个元素的平均值：

final_result= vapply(rbinom_result, mean, FUN.VALUE = 1)
# system.time result:
# user  system elapsed 
# 23.372   0.050  23.412

如何让这两位代码运行更快？预先感谢您的帮助。

编辑： 我需要分别计算 rbinom_result 的每个元素的平均值，这可能是 NA 或可变数量的整数，具体取决于传递给的数字rbinom_NA() 由向量 d。所以实际上我会用从不同的 lambda 创建的 d 向量进行计算，所以 rbinom_result 现在看起来像这样：

#if lamb=0.01
head(rbinom_result)
$pos1
[1] NA

$pos2
[1] NA

$pos3
[1] NA
#etc

#if lamb=5
head(rbinom_result)

 $pos1
 [1] 0 0 0 0 0 0 0
 
 $pos2
 [1] 1 1 1 0 0 1 1
 
 $pos3
 [1] 0 0 0 0 0 0
# etc

然后最后我想将 final_result 整理成一个数据框，其中包含每个 lambda 值的结果（lambda 是列，有 8e6 行）。所以整个事情看起来像这样：

library(parallel)
lambdas = c(0.01,5)#just two here, will be more
final_result = vector(mode = "list", length = length(lambdas))

for (lamb in lambdas){
b = runif(n = 8e6, min = 0, max = 1)
d=rpois( 8e6, lambda =lamb)

rbinom_NA = function(x , y) {
  result = rbinom(n = x,
                  size = 1,
                  prob = y) 
  
  if (length(result) == 0)
    return(NA)
  return(result)
}


rbinom_result = vector(mode = "list", length = 8e6)
rbinom_result = mcmapply(
  d,
  b,
  FUN = rbinom_NA, 
  mc.cores = detectCores()-1) #run on Ubuntu, 32 cores and 64Gb memory

names(rbinom_result) = paste0("pos",1:8e6)

final_result[[as.character(lamb)]]= vapply(rbinom_result, mean, FUN.VALUE = 1)

}
final_result = do.call("cbind", final_result)

最后：

 head(final_result)
     0.01         5
pos1   NA 0.0000000
pos2   NA 0.7142857
pos3   NA 0.0000000
pos4   NA 0.6000000
pos5   NA 1.0000000
pos6   NA 0.1666667
#etc

Answer 1

您可以通过涉及 rbinom:

的单个命令为每个 lambda 获取 final_result

b <- runif(8e6)
d <- rpois(8e6, 0.1)
system.time(final_result <- rbinom(length(b), d, b)/d)
#   user  system elapsed 
#  0.327   0.058   0.384

之所以有效，是因为以下两行遵循原子 x 和 y 的相同分布：

mean(rbinom(x, 1, y))
rbinom(1, x, y)/x

不同之处在于，只需将 rbinom 中的 n（第一个）参数从1 到 length(x)。此外，如果 x = 0、rbinom(1, x, y)/x 为 NA，则 NA 行为将保持不变。

对于多个 lambda:

library(parallel)

fSim <- function(lambda, n){
  b <- runif(n)
  d <- rpois(n, lambda)
  return(rbinom(length(b), d, b)/d)
}

lambdas <- c(0.1, 5)
system.time(final_result <- mcmapply(fSim, lambdas, 8e6, mc.cores = min(length(lambdas), detectCores() - 1)))
#   user  system elapsed 
#  2.340   0.931   2.031

使用带有 NA 和 mcmapply 的 rbinom 进行优化

Optimization using rbinom with NAs and mcmapply

optimization

r