logo标签列表

比较数字向量与另一个向量中的值的最佳方法

rvector
7
假设我有两个数值向量： 

a <- c(1,3,4,5,6,7,3)
b <- c(3,5,1,3,2)

我希望对每个输入的a应用一些函数FUN，并针对整个b进行操作，最有效的方法是什么。

更具体地说，在这种情况下，对于a中的每个元素，我想知道对于每个'a'值，有多少b中的元素大于或等于该值。天真的方法是执行以下操作：

sum(a < b)

当然，这样做是行不通的，因为它试图并行迭代每个向量，并给出警告：

长对象长度不是短对象长度的倍数

顺便说一下，该命令的输出是3。
然而，在我的情况下，我想要看到的输出是：
0 2 4 4 5 5 2

当然，我意识到可以使用for循环来实现这个功能：
out <- c()
for (i in a) {
    for (i in a) { out[length(out) + 1] = sum(b<i)}
}

同样的，我可以使用sapply这样做：
sapply(a, function(x)sum(b<x))

然而，我试图成为一名优秀的R程序员并且避免使用for循环，sapply看起来速度非常慢。是否有其他替代方案？

就我的情况而言，我需要执行几百万次操作，其中b的长度始终小于a的长度，而a的长度范围从1到30不等。

在您的一百万次迭代中，向量a和b是否都会变化，还是其中一个固定不变？ - Prasad Chalasani
两个向量在遍历所有数据时会多次生成。因此，它们确实会有所变化，每个向量大约有10000个不同的值。 - Pridkett
4个回答
6

试试这个：

findInterval(a - 0.5, sort(b))

通过避免使用 sort 和简化 .Internal 封装来避免在 findInterval 和 order 中的开销，从而提高速度：

order2 = function(x) .Internal(order(T, F, x))

findInterval2 = function(x, vec, rightmost.closed=F, all.inside=F) {
  nx <- length(x)
  index <- integer(nx)
  .C('find_interv_vec', xt=as.double(vec), n=length(vec),
    x=as.double(x), nx=nx, as.logical(rightmost.closed),
    as.logical(all.inside), index, DUP = FALSE, NAOK=T,
    PACKAGE='base')
  index
}

> system.time(for (i in 1:10000) findInterval(a - 0.5, sort(b)))
   user  system elapsed 
   1.22    0.00    1.22 
> system.time(for (i in 1:10000) sapply(a, function(x)sum(b<x)))
   user  system elapsed 
   0.79    0.00    0.78 
> system.time(for (i in 1:10000) rowSums(outer(a, b, ">")))
   user  system elapsed 
   0.72    0.00    0.72 
> system.time(for (i in 1:10000) findInterval(a - 0.5, b[order(b)]))
   user  system elapsed 
   0.42    0.00    0.42 
> system.time(for (i in 1:10000) findInterval2(a - 0.5, b[order2(b)]))
   user  system elapsed 
   0.16    0.00    0.15

如果你需要大量迭代，且N相对较小，那么定义findInterval2和order2可能是有必要的。

另外，针对更大的N，可以进行时间测试：

> a = rep(a, 100)
> b = rep(b, 100)
> system.time(for (i in 1:100) findInterval(a - 0.5, sort(b)))
   user  system elapsed 
   0.01    0.00    0.02 
> system.time(for (i in 1:100) sapply(a, function(x)sum(b<x)))
   user  system elapsed 
   0.67    0.00    0.68 
> system.time(for (i in 1:100) rowSums(outer(a, b, ">")))
   user  system elapsed 
   3.67    0.26    3.94 
> system.time(for (i in 1:100) findInterval(a - 0.5, b[order(b)]))
   user  system elapsed 
      0       0       0 
> system.time(for (i in 1:100) findInterval2(a - 0.5, b[order2(b)]))
   user  system elapsed 
      0       0       0
@Charles：这比OP的sapply解决方案慢--我通过执行system.time( replicate( 10000, ... ) )进行了测试。 - Prasad Chalasani
对于小的a/b，排序占据了运行时间的主导地位。在我的计算机上，b[order(b)]的基准测试速度约为sort(b)的4倍，并且比sapply原始版本快两倍。对于更大的a/b（rep（100）），使用findInterval可以获得数量级的改进。 - Charles
大部分剩余的开销都在findInterval中的is.sorted和is.na检查以及order中的各种检查等方面。如果您在这些函数周围定义更薄的包装器（例如，order2 = function(x) .Internal(order(T, F, x))），那么它会再快3倍左右。 - Charles
干得好，@Charles +1 - 一个很好的解决方案，使用了我很少遇到的函数。你能否发布你的findInterval2和order2函数来完善答案？ - Gavin Simpson
5

一种选择是使用outer()将二元运算符函数>应用于a和b：

> outer(a, b, ">")
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] FALSE FALSE FALSE FALSE FALSE
[2,] FALSE FALSE  TRUE FALSE  TRUE
[3,]  TRUE FALSE  TRUE  TRUE  TRUE
[4,]  TRUE FALSE  TRUE  TRUE  TRUE
[5,]  TRUE  TRUE  TRUE  TRUE  TRUE
[6,]  TRUE  TRUE  TRUE  TRUE  TRUE
[7,] FALSE FALSE  TRUE FALSE  TRUE

答案可以通过上面结果的行求和得出：
> rowSums(outer(a, b, ">"))
[1] 0 2 4 4 5 5 2

对于这个数据集，这个解决方案比findIntervals()稍微快一些，但差别不大。
> system.time(replicate(1000, findInterval(a - 0.5, sort(b))))
   user  system elapsed 
  0.131   0.000   0.132 
> system.time(replicate(1000, rowSums(outer(a, b, ">"))))
   user  system elapsed 
  0.078   0.000   0.079

这个版本比sapply()的版本稍微快一点，但是差别不大：

> system.time(replicate(1000, sapply(a, function(x)sum(b<x))))
   user  system elapsed 
  0.082   0.000   0.082

@Charles指出，findInterval()示例中大部分时间都用于sort()，可以通过使用order()来避免这种情况。当这样做时，findInterval()解决方案比outer()解决方案更快：

> system.time(replicate(1000, findInterval(a - 0.5, b[order(b)])))
   user  system elapsed 
  0.049   0.000   0.049
1

在生产代码中使用 R 的内部机制时，我会非常谨慎。这些内部机制在不同版本之间很容易发生变化。

sort.int 比 sort 更快 - 而且 b[order(b)] 比 sort.int(b) 更快，这很奇怪。R 肯定可以改进其排序算法...

除非您使用 R 的内部机制，否则似乎使用 vapply 实际上更快：

> system.time(for (i in 1:10000) findInterval(a - 0.5, sort(b)))
   user  system elapsed 
   0.99    0.00    0.98 
> system.time(for (i in 1:10000) findInterval(a - 0.5, sort.int(b)))
   user  system elapsed 
    0.8     0.0     0.8 
> system.time(for (i in 1:10000) findInterval(a - 0.5, b[order(b)]))
   user  system elapsed 
   0.32    0.00    0.32 
> system.time(for (i in 1:10000) sapply(a, function(x)sum(b<x)))
   user  system elapsed 
   0.61    0.00    0.59 
> system.time(for (i in 1:10000) vapply(a, function(x)sum(b<x), 0L))
   user  system elapsed 
   0.18    0.00    0.19
0

仅作为补充说明：如果您知道每个向量的值范围，那么先计算最大值和最小值可能会更快，例如：

order2 = function(x) .Internal(order(T, F, x))
findInterval2 = function(x, vec, rightmost.closed=F, all.inside=F) {
  nx <- length(x)
  index <- integer(nx)
  .C('find_interv_vec', xt=as.double(vec), n=length(vec),
    x=as.double(x), nx=nx, as.logical(rightmost.closed),
    as.logical(all.inside), index, DUP = FALSE, NAOK=T,
    PACKAGE='base')
  index
}

f <- function(a, b) {
  # set up vars
  a.length <- length(a)
  b.length <- length(b)
  b.sorted <- b[order2(b)]
  b.min <- b.sorted[1]
  b.max <- b.sorted[b.length]
  results <- integer(a.length)

  # pre-process minimums
  v.min <- which(a <= b.min)

  # pre-process maximums
  v.max <- which(a > b.max)
  results[v.max] <- b.max

  # compare the rest
  ind <- c(v.min, v.max)
  results[-ind] <- findInterval2(a[-ind] - 0.5, b.sorted)
  results
}

这将给出以下时间

> N <- 10
> n <- 1e5
> b <- runif(n, 0, 100)
> a <- runif(n, 40, 60) # NB smaller range of values than b
> summary( replicate(N, system.time(findInterval2(a - 0.5, b[order2(b)]))[3]) )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0300  0.0300  0.0400  0.0390  0.0475  0.0500 
> summary( replicate(N, system.time(f(a, b))[3]) )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.010   0.030   0.030   0.027   0.030   0.040

然而，如果您事先不知道范围，或者无法对其进行合理猜测，则这可能会更慢。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接