比使用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

@ db的方法:

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建议的那样,data.table包可能是一个很好的选择,如果findInterval它太慢,因为它有很多专注于排序数据的功能,但对我来说如何在这里应用它并不明显.