r - 比使用 for 循环更快的编码

Question

假设我有以下数据框

set.seed(36)    

n <- 300

dat <- data.frame(x = round(runif(n,0,200)), y =  round(runif(n, 0, 500)))
d <- dat[order(dat$y),]

对于 d$y<=300 的每个值，我必须创建一个变量 res其中分子是指标 (d$x <= d$y[i]) 的总和分母是指标 (d$y >= d$y[i]) 的总和。我已经将代码写在for loop中:

res <- NULL

for( i in seq_len(sum(d$y<=300)) ){

    numerator   <- sum(d$x <= d$y[i])
    denominator <- sum(d$y >= d$y[i])

    res[i] <- numerator / denominator
 } 

但我担心的是 x 的观察数量何时和y大，即数据框的行数增加，for loop会慢慢地工作。此外，如果我模拟数据 1000 次并且每次都运行 for loop ，程序效率会很低。

什么是更有效的代码解决方案？

Answer 1

这取决于 d 已经按原样排序:

# example data
set.seed(36)    
n <- 1e5
dat <- data.frame(x = round(runif(n,0,200)), y =  round(runif(n, 0, 500)))
d <- dat[order(dat$y),]

我的建议(感谢@alexis_laz 提供分母):

system.time(res3 <- {
  xs <- sort(d$x)                                 # sorted x
  yt <- d$y[d$y <= 300]                           # truncated y

  num = findInterval(yt, xs)
  den = length(d$y) - match(yt, d$y) + 1L

  num/den
})
#    user  system elapsed 
#       0       0       0 

OP的方法:

system.time(res <- {
  res <- NULL
  for( i in seq_len(sum(d$y<=300)) ){
    numerator   <- sum(d$x <= d$y[i])
    denominator <- sum(d$y >= d$y[i])
    res[i] <- numerator / denominator
  }
  res
})
#    user  system elapsed 
#   50.77    1.13   52.10 

# verify it matched
all.equal(res,res3) # TRUE

@d.b 的方法:

system.time(res2 <- {
  numerator = rowSums(outer(d$y, d$x, ">="))
  denominator = rowSums(outer(d$y, d$y, "<="))
  res2 = numerator/denominator

  res2 = res2[d$y <= 300]
  res2
})
# Error: cannot allocate vector of size 74.5 Gb

# ^ This error is common when using outer() on large-ish problems

矢量化。通常，如果可以矢量化，R 中的任务会更快。与有序向量相关的关键函数的名称很困惑(findInterval、sort、order 和 cut)，但幸运的是它们都适用于向量。

连续与离散。无论数据是连续的还是具有质量点/重复值，上面的 match 应该是计算分母的快速方法。如果数据是连续的(因此没有重复)，则分母可以是 seq(length(xs), length = length(yt), by=-1)。如果它是完全离散的并且有大量重复(如此处的示例)，则可能还有某种方法可以使其更快，也许像其中之一:

 den2 <- inverse.rle(with(rle(yt), list(
    values = length(xs) - length(yt) +  rev(cumsum(rev(lengths))), 
    lengths = lengths)))

 tab  <- unname(table(yt))
 den3 <- rep(rev(cumsum(rev(tab))) + length(xs) - length(yt), tab)

 # verify
 all.equal(den,den2) # TRUE
 all.equal(den,den3) # TRUE

findInterval 仍然适用于连续数据的分子。我猜这对于这里考虑的重复值情况并不理想(因为我们冗余地寻找许多重复 yt 值的间隔)。类似的加速想法可能也适用。

其他选项。 正如 @chinsoon 所建议的，如果 findInterval 太慢，则 data.table 包可能是一个不错的选择，因为它有很多功能集中于已排序的数据，但对我来说如何应用它并不明显。