我有一个大小为10000 x 100的矩阵和一个长度为100的向量。我想将一个自定义函数percentry应用于矩阵的每一列,该函数接受向量参数和标量参数,这样在迭代j时,与percentry一起使用的参数是矩阵的列j和向量的条目j-。有没有办法使用apply函数之一来执行此操作?
这是我的密码。它运行,但没有返回正确的结果。
percentile <- function(x, v){
length(x[x <= v]) / length(x)
}
X <- matrix(runif(10000 * 100), nrow = 10000, ncol = 100)
y <- runif(100)
result <- apply(X, 2, percentile, v = y)
我一直在使用的解决方法是将y附加到X,然后重新编写百分比函数,如下所示。
X <- rbind(X, y)
percentile2 <- function(x){
v <- x[length(x)]
x <- x[-length(x)]
length(x[x <= v]) / length(x)
}
result <- apply(X, 2, percentile2)
这段代码确实返回了正确的结果,但我更喜欢更优雅的代码。
如果您了解R
是矢量化的,并且知道正确的函数,那么您可以完全避免循环,并在一行相对简单的代码中完成整个操作。。。
colSums( t( t( X ) <= y ) ) / nrow( X )
通过向量化,R将在X
的每一列中循环使用y
中的每个元素(默认情况下,它将在各行中循环使用,因此我们使用转置函数t
将列转换为行,应用逻辑比较<=
,然后再次转置回
由于TRUE
和FALSE
的值分别为1和0,我们可以使用colSums
有效地获得每列中满足条件的行数,然后将每列除以总行数(记住回收规则!)。结果完全一样。。。。
res1 <- apply(X2, 2, percentile2)
res2 <- colSums( t( t( X ) <= y ) ) / nrow( X )
identical( res1 , res2 )
[1] TRUE
很明显,因为这不使用任何R循环,所以速度更快(在这个小矩阵上大约是10倍)。
更好的方法是这样使用rowMeans
(感谢@flodel):
rowMeans( t(X) <= y )
我认为最简单、最清晰的方法是使用for
循环:
result2 <- numeric(ncol(X))
for (i in seq_len(ncol(X))) {
result2[i] <- sum(X[,i] <= y[i])
}
result2 <- result2 / nrow(X)
我能想到的最快、最短的解决方案是:
result1 <- rowSums(t(X) <= y) / nrow(X)
SimonO101在他的回答中解释了这是如何工作的。正如我所说,它很快。然而,缺点是不太清楚这里到底计算了什么,尽管您可以通过将这段代码放在一个命名良好的函数中来解决这个问题。
flodel还提出了一种使用CCD_ 11的解决方案,CCD_。然而,要做到这一点,您首先需要将每个列或矩阵放入list
或data.frame
:中
result3 <- mapply(percentile, as.data.frame(X), y)
就速度而言(请参阅下面的一些基准测试),for循环并没有那么糟糕,而且它比使用apply
更快(至少在这种情况下)。使用rowSums
和矢量回收的技巧更快,是使用apply
的解决方案的10倍以上。
> X <- matrix(rnorm(10000 * 100), nrow = 10000, ncol = 100)
> y <- runif(100)
>
> system.time({result1 <- rowSums(t(X) <= y) / nrow(X)})
user system elapsed
0.020 0.000 0.018
>
> system.time({
+ X2 <- rbind(X, y)
+ percentile2 <- function(x){
+ v <- x[length(x)]
+ x <- x[-length(x)]
+ length(x[x <= v]) / length(x)
+ }
+ result <- apply(X2, 2, percentile2)
+ })
user system elapsed
0.252 0.000 0.249
>
>
> system.time({
+ result2 <- numeric(ncol(X))
+ for (i in seq_len(ncol(X))) {
+ result2[i] <- sum(X[,i] <= y[i])
+ }
+ result2 <- result2 / nrow(X)
+ })
user system elapsed
0.024 0.000 0.024
>
> system.time({
+ result3 <- mapply(percentile, as.data.frame(X), y)
+ })
user system elapsed
0.076 0.000 0.073
>
> all(result2 == result1)
[1] TRUE
> all(result2 == result)
[1] TRUE
> all(result3 == result)
[1] TRUE