Оптимальный метод сравнения вектора чисел со значениями в другом векторе - PullRequest
4 голосов
/ 24 февраля 2011

Предположим, у меня есть два вектора значений:

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

И я хочу применить некоторую функцию FUN к каждому из входов a по отношению ко всему 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 и держусь подальше от циклов, и sapply кажется очень медленным. Есть ли другие альтернативы?

Что бы это ни стоило, я делаю это пару миллионов раз, когда length(b) всегда меньше length(a), а length(a) колеблется от 1 до 30.

Ответы [ 4 ]

4 голосов
/ 24 февраля 2011

Попробуйте это:

findInterval(a - 0.5, sort(b))

Улучшение скорости из a) избегая sort и b) избежания накладных расходов в findInterval и order с использованием более простых .Internal упаковщиков:

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 

Сложность определения findInterval2 и order2, вероятно, оправдана только если у вас есть куча итераций с довольно маленьким N.

Также время для большего 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 
3 голосов
/ 25 февраля 2011

Один из вариантов - использовать 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

@ Чарльз отмечает, что большую часть времени в примере 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
0 голосов
/ 30 марта 2011

Я бы очень опасался использовать внутреннюю часть 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 голосов
/ 25 февраля 2011

Просто дополнительное примечание: если вы знаете диапазон значений для каждого вектора, то, возможно, быстрее будет сначала рассчитать максимальное и минимальное значения, например,

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 

Однако, если вы не знаете диапазоны заранее или не можете сделать обоснованное предположение о них, то это, вероятно, будет медленнее.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...