arrays - 将函数应用于多维数组 : R vs MATLAB

Question

这个问题可以认为与this one相关，这帮助我提高了 R 在计算大数组均值时的性能。不幸的是，在这种情况下，我试图应用更复杂的东西(比如分位数计算)。

我有一个包含超过 4000 万个元素的 4 维数组，我想计算特定维度上的第 66 个百分位数。这里有 MATLAB 代码:

> n = randn(100, 50, 100, 20);
> tic; q = quantile(n, 0.66, 4); toc
Elapsed time is 0.440824 seconds.

让我们在 R 中做一些类似的事情。

> n = array(rnorm(100*50*100*20), dim = c(100,50,100,20))
> start = Sys.time(); q = apply(n, 1:3, quantile, .66); print(Sys.time() - start)
Time difference of 1.600693 mins

我知道 MATLAB wrt R 的性能更好，但在这种情况下我不知道该怎么做。可能我只需要等待 2 分钟而不是一秒钟...... 我希望有人可以建议我任何改善运行时间的方法， 无论如何，先谢谢你...

更新 我已将一些建议应用到评论中并减少了运行时间:

> start = Sys.time(); q = apply(n, 1:3, quantile, .66, names = FALSE); print(Sys.time() - start)
Time difference of 33.42773 secs

我们距离 MATLAB 性能还有很长的路要走，但至少我学到了一些东西。

更新 我在这里提出了一些与讨论的“分位数”功能相关的进步 here .我上面显示的相同代码的运行时间已从 33 秒变为 5 秒...

Answer 1

的 RcppOctave 包调用 GNU Octave API函数； 如果您还不知道 GNU Octave，它与 Matlab非常相似，旨在实现完全兼容。

这几乎和 Matlab direct 一样快...

library(RcppOctave)

set.seed(1)
n = array(rnorm(100*50*100*20), dim = c(100,50,100,20))

system.time( s <- octave_style_quantile(n, .66, 4) )
##    user  system elapsed 
##   0.526   0.048   0.574

# *R* `quantile` argument `type=5` is the method that matches matlab.
system.time( r <- apply(n, 1:3, quantile, .66, names = FALSE, type=5) )
##    user  system elapsed 
##  23.308   0.029  23.327

dim(r)
## [1] 100  50 100

identical(r,s)
## [1] TRUE

---

Octave 的相当直接的翻译 quantile.m 到 R.

octave_style_quantile <- function (x, p=NULL, dim=NULL) {
  if ( is.null(p) ) p <- c(0.00, 0.25, 0.50, 0.75, 1.00)

  if ( is.null(dim) ) {
    ## Find the first non-singleton dimension.
    dim <- which(dim(x) > 1)[1];
  }

  stopifnot( is.numeric(x)||is.logical(x),
             is.numeric(p),
             dim <= length(dim(x)) )

  ## Set the permutation vector.
  perm <- seq_along(dim(x))
  perm[1] <- dim
  perm[dim] <- 1

  ## Permute dim to the 1st index.
  x <- aperm(x, perm);

  ## Save the size of the permuted x N-d array.
  sx = dim(x);

  ## Reshape to a 2-d array.
  dim(x) <- c( sx[1], prod(sx[-1]) );

  ## Calculate the quantiles.
  q = .CallOctave("quantile",x,p)

  ## Return the shape to the original N-d array.
  dim(q) <- c( length(p), sx[-1] )

  ## Permute the 1st index back to dim.
  q = aperm(q, perm);
  if( any(dim(q)==1) ) dim(q) <- dim(q)[-which(dim(q)==1)]
  q
}