将函数应用于多维数组：R 与 MATLAB

这个问题可以被认为与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 分钟而不是一秒钟...... 我希望有人能建议我任何提高跑步时间的方法， 不管怎样，先谢谢你了……

UPDATE我已将一些建议应用到评论中，并减少了运行时间：

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

我们距离 MATLAB 的表现还很远，但至少我学到了一些东西。

UPDATE我在这里提出了与讨论的“分位数”函数相关的一些进展here。我上面显示的相同代码的运行时间已从 33 秒缩短到 5 秒......

The 八度音程包调用GNU 倍频程API函数； 如果你还不知道GNU 倍频程, it is very如同Matlab并致力于完全兼容。

这几乎和直接 Matlab 一样快......

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 的相当直接的翻译分位数.m to 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
}