向量化包含一個語句和函數的for循環

我正在嘗試vectorise的一個可重複的代碼示例。向量化包含一個語句和函數的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's sg（）。關鍵在於意識到它只是每一行中兩個值中最低的值。一旦出現了這種心理跳躍，就只剩下一個單獨的矢量來處理矢量化問題了。

# cutOff should be lower than the lowest of Sepal.Length & Sepal.Width 
    m <- pmin(iris$Sepal.Length, iris$Sepal.Width)

來源

2015-05-05 Henry E

我喜歡加上另一個答案：

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 
}

這種方法避免了一個或forouter環和比@ Nicola的方法快4倍時間：

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

來源

2015-05-05 16:14:13 sgibb

確實用'findInterval'（+1）。這也是我的出發點，但是我失敗了，結果出現了一個更加複雜的「切割」代碼。 – Henrik

這不會刪除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])) 
}

編輯注：由@nicola包含的做法，現在是最快的

來源

2015-05-05 14:14:17

雖然我喜歡@nicola智能解決方案，我更喜歡'dd'因爲'outer'是內存密集型很長'cutOffs'。 – ExperimenteR

Tx用於將我的解決方案包含在基準測試中。 – nicola

我想是這樣的：

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

來源

2015-05-05 14:14:23

您可以使用dplyr

library(dplyr) 
data_frame(y=cutOffs) %>% 
    rowwise() %>% 
    mutate(x = sum(iris$Sepal.Length > y & iris$Sepal.Width > y))

來源

2015-05-05 14:31:17 ExperimenteR

您可以使用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]))

來源

2015-05-05 14:42:16 nicola

基於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

即四個行與兩個值> 2，兩排，這兩個值> 4，et.c

來源

2015-05-05 16:02:43 Henrik

向量化包含一個語句和函數的for循環

回答

相關問題