Векторизация цикла for, содержащего оператор и функцию
Воспроизводимый пример кода, который я пытаюсь прорисовать.
cutOffs <- seq(1,10,0.2)
plotOutput <- matrix(nrow=length(cutOffs), ncol=2)
colnames(plotOutput) <- c("x","y")
plotOutput[,"y"] <- cutOffs
for(plotPoint in 1:length(cutOffs))
{
plotOutput[plotPoint, "x"] <-
nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] &
iris$Sepal.Width > cutOffs[plotPoint]), ])
}
plotOutput
В частности, я хочу узнать, есть ли способ для векторизации этой части.
nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] &
iris$Sepal.Width > cutOffs[plotPoint]), ])
Скажем, я должен был использовать библиотеку plyr или какую-то форму приложения, вероятно, не так много ускоряется, и это действительно то, что я ищу. В сущности, я пытаюсь понять, есть ли какая-то техника для векторизации, которую я пропустил или сумел пропустить во время поиска.
UPDATE:
Unit: milliseconds
expr min lq mean median uq max neval
op() 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700 1
jr() 3976.53088 3976.53088 3976.53088 3976.53088 3976.53088 3976.53088 1
dd() 4253.21050 4253.21050 4253.21050 4253.21050 4253.21050 4253.21050 1
exp() 5085.45331 5085.45331 5085.45331 5085.45331 5085.45331 5085.45331 1
nic() 8719.82043 8719.82043 8719.82043 8719.82043 8719.82043 8719.82043 1
sg() 16.66177 16.66177 16.66177 16.66177 16.66177 16.66177 1
Более реалистичное приближение того, что я на самом деле делаю, это
# generate data
numObs <- 1e5
iris <- data.frame( Sepal.Length = sample(1:numObs), Sepal.Width = sample(1:numObs) )
cutOffs <- 1:(numObs*0.01)
plotOutput <- matrix(nrow=length(cutOffs), ncol=2)
colnames(plotOutput) <- c("x","y")
plotOutput[,"y"] <- cutOffs
за которым следует какой-либо конкретный метод, который предпочитают.
В общем случае он будет использоваться на наборах данных с 50 000 - 200 000 точек.
Был большой прыжок от использования
sum(Sepal.Length > cutOffs[plotPoint] & Sepal.Width > cutOffs[plotPoint])
это то, чего я раньше не встречал, как более оптимальный подход.
В то же время лучшим ответом является sgibb sg(). Ключ понимает, что это имеет значение только самое низкое из двух значений в каждой строке. После того, как этот умственный скачок был сделан, только один вектор, оставленный для обработки и векторизации, достаточно прост.
# cutOff should be lower than the lowest of Sepal.Length & Sepal.Width
m <- pmin(iris$Sepal.Length, iris$Sepal.Width)
Ответы
Ответ 1
Мне нравится добавить еще один ответ:
sg <- function() {
# cutOff should be lower than the lowest of Sepal.Length & Sepal.Width
m <- pmin(iris$Sepal.Length, iris$Sepal.Width)
ms <- sort.int(m)
# use `findInterval` to find all the indices
# (equal to "how many numbers below") lower than the threshold
plotOutput[,"x"] <- length(ms)-findInterval(cutOffs, ms)
plotOutput
}
Этот подход позволяет избежать цикла for
или outer
и в 4 раза быстрее, чем подход @nicola:
microbenchmark(sg(), nic(), dd())
#Unit: microseconds
# expr min lq mean median uq max neval
# sg() 88.726 104.5805 127.3172 123.2895 144.2690 232.441 100
# nic() 474.315 526.7780 625.0021 602.3685 706.7530 997.412 100
# dd() 669.841 736.7800 887.4873 847.7730 976.6445 2800.930 100
identical(sg(), dd())
# [1] TRUE
Ответ 2
Вы можете использовать outer
:
plotOutput[,"x"]<-colSums(outer(1:nrow(iris),1:length(cutOffs),function(x,y) iris$Sepal.Length[x] > cutOffs[y] & iris$Sepal.Width[x] > cutOffs[y]))
Ответ 3
Это не удаляет цикл for
, но я предполагаю, что он даст вам некоторое ускорение - не стесняйтесь сравнивать и дайте нам знать, как он сравнивается с вашими реальными данными:
for(i in seq_along(cutOffs)) {
x <- cutOffs[i]
plotOutput[i, "x"] <- with(iris, sum(Sepal.Length > x & Sepal.Width > x))
}
Здесь немного теста с использованием данных образца (который, возможно, крошечный, но может дать некоторые указания):
library(microbenchmark)
microbenchmark(op(), jr(), dd(), exp(), nic())
Unit: microseconds
expr min lq median uq max neval
op() 6745.428 7079.8185 7378.9330 9188.0175 11936.173 100
jr() 1335.931 1405.2030 1466.9180 1728.6595 4692.748 100
dd() 684.786 711.6005 758.7395 923.6670 4473.725 100
exp() 1928.083 2066.0395 2165.6985 2392.7030 5392.475 100
nic() 383.007 402.5495 439.3835 541.6395 851.488 100
Функции, используемые в эталонном тесте, определяются следующим образом:
op <- function(){
for(plotPoint in 1:length(cutOffs))
{
plotOutput[plotPoint, "x"] <-
nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] &
iris$Sepal.Width > cutOffs[plotPoint]), ])
}
plotOutput
}
jr <- function() {
cbind(x = sapply(cutOffs, counts), y = plotOutput[,"y"])
}
dd <- function() {
for(i in seq_along(cutOffs)) {
x <- cutOffs[i]
plotOutput[i, "x"] <- with(iris, sum(Sepal.Length > x & Sepal.Width > x))
}
plotOutput
}
exp <- function() {
data_frame(y=cutOffs) %>%
rowwise() %>%
mutate(x = sum(iris$Sepal.Length > y & iris$Sepal.Width > y))
}
nic <- function() {
plotOutput[,"x"]<-colSums(outer(1:nrow(iris),1:length(cutOffs),function(x,y) iris$Sepal.Length[x] > cutOffs[y] & iris$Sepal.Width[x] > cutOffs[y]))
}
Отредактируйте примечание: включенный подход by @nicola, который теперь самый быстрый
Ответ 4
Вы можете использовать dplyr
library(dplyr)
data_frame(y=cutOffs) %>%
rowwise() %>%
mutate(x = sum(iris$Sepal.Length > y & iris$Sepal.Width > y))
Ответ 5
Я думаю, что-то вроде:
counts <- function(x) sum(iris$Sepal.Length > x & iris$Sepal.Width > x )
cbind(x = sapply(cutOffs, counts), y = plotOutput[,"y"])
и просто проверить:
res <- cbind(x=sapply(cutOffs,counts), y=plotOutput[,"y"])
identical(plotOutput,res)
[1] TRUE
Ответ 6
Другая возможность, основанная на pmin
, cut
и table
brk <- c(cutOffs, Inf)
rev(cumsum(rev(table(cut(pmin(iris$Sepal.Length, iris$Sepal.Width), brk)))))
Небольшой пример, который может быть проще использовать, если вы хотите работать через код "изнутри":
set.seed(1)
df <- data.frame(x = sample(1:10, 6), y = sample(1:10, 6))
cutOffs <- seq(from = 2, to = 8, by = 2)
brk <- c(cutOffs, Inf)
rev(cumsum(rev(table(cut(pmin(df$x, df$y), brk)))))
# (2,4] (4,6] (6,8] (8,Inf]
# 4 2 1 0
I.e., четыре строки с обоими значениями > 2, две строки с обоими значениями > 4, et.c