1

我有一个大小为 10000 x 100 的矩阵和一个长度为 100 的向量。我想对矩阵的每一列应用一个自定义函数percentile，它接受一个向量参数和一个标量参数，这样在迭代j，与percentile一起使用的参数是矩阵的第j列和向量的第j列。有没有办法使用其中一个应用功能来做到这一点？

这是我的代码。它运行，但不返回正确的结果。

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，然后重新编写 percentile 函数，如下所示。

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)

这段代码确实返回了正确的结果，但我更喜欢更优雅的东西。

4

2 回答 2

2

如果您了解它R是矢量化的并且知道正确的功能，您可以完全避免循环，并在一个相对简单的行中完成整个事情......

 colSums(  t( t( X ) <= y ) ) / nrow( X )

y通过矢量化，R 将在每一列中回收每个元素X（默认情况下，它将跨行执行此操作，因此我们使用转置函数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  )
于 2013-08-23T09:55:53.257 回答
2

我认为最简单和最清晰的方法是使用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 在他的回答中解释了这是如何工作的。正如我所说，它很快。但是，缺点是不太清楚这里到底计算了什么，尽管您可以通过将这段代码放在一个命名良好的函数中来解决这个问题。

florel 还建议了mapply一个apply可以在多个向量上工作的解决方案。但是，要使其工作，您首先需要将每个列或矩阵放在 a listor中data.frame：

result3 <- mapply(percentile, as.data.frame(X), y)

在速度方面（请参阅下面的一些基准测试）for-loop 并没有那么糟糕，而且它比使用apply（至少在这种情况下）更快。向量回收的技巧rowSums更快，比使用apply.

> 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
于 2013-08-23T09:19:41.110 回答