San*_*ado 10 arrays r data.table tidyverse
我经常使用大型 3D 数组（纬度、经度和时间），例如 720x1440x480 的大小。通常，我需要对每个纬度和经度随时间进行操作，例如，获取平均值（导致 2D 数组）或及时获取滚动平均值（导致 3D 数组），或更复杂的函数。
我的问题是：哪种包装（或方式）最有效和最快？
我知道一个选项是基础 R，使用 apply 函数和滚动函数与提供 rollapply 函数的包 zoo 混合。另一种方式是使用 tidyverse，另一种方式是使用 data.table。以及这些包之间的组合。但有没有最快的？
例如，如果我有这个数据立方体：
data <- array(rnorm(721*1440*480),dim = c(721,1440,480))
哪些维度是纬度、经度和时间，如下所示：
lat <- seq(from = -90, to = 90, by = 0.25)
lon <- seq(from = 0, to = 359.75, by = 0.25)
time <- seq(from = as.Date('1980-01-01'), by = 'month', length.out = 480)
我通常需要做这样的事情（这是在基础 R + 动物园）：
# Average in time
average_data <- apply(data, 1:2, mean)
# Rolling mean, width of window = 3
library(zoo)
rolling_mean <- function(x){
return(rollapply(data = x, width = 3, by = 1, FUN = mean))
}
rolling_mean_data <- apply(X = data, MARGIN = 1:2,
FUN = rolling_mean)
rolling_mean_data <- aperm(a = rolling_mean_data,perm = c(2,3,1))
函数可能改变并不总是意味着，也可能是其他统计数据，如标准偏差或与时间序列的相关性。
那么，进行此类微积分的最快方法是什么？
通常data.table是相当快和内存效率。您可能需要检查?data.table::frollapply并注意，它as.data.table为类“数组”提供了 S3 方法。
下面将您提供的基本 R 代码片段与其data.table对应代码片段进行比较（我对 tidyverse 不太熟悉）。我只关注计算本身而不关注输出格式 - 因此对于data.table解决方案，输出不会转换回数组。这是关于此的一些信息 - 也aperm像上面一样使用。
基准测试是使用简化的数据集创建的（计算完成后，我将使用数据集的结果进行更新）。
编辑： 刚刚使用整个数据集更新了基准测试结果。
library(data.table)
library(microbenchmark)
library(zoo)
lat <- seq(from = -90, to = 90, by = 0.25)
lon <- seq(from = 0, to = 359.75, by = 0.25)
time <- seq(from = as.Date('1980-01-01'), by = 'month', length.out = 480)
# data <- array(rnorm(721 * 1440 * 480), dim = c(721, 1440, 480), dimnames = list(lat = lat, lon = lon, time = time))
data <- array(rnorm(72 * 144 * 48), dim = c(72, 144, 48), dimnames = list(lat = NULL, lon = NULL, time = NULL)) # reduced dataset
DT <- as.data.table(data)
setorder(DT, time) # as per @jangorecki we need to ensure that DT is sorted by time
rolling_mean <- function(x) {
return(rollapply(
data = x,
width = 3,
by = 1,
FUN = mean
))
}
microbenchmark(
"mean - base R" = {
apply(data, 1:2, mean)
},
"mean - data.table" = {
DT[, .("mean_value" = mean(value)), by = .(lat, lon)]
},
times = 1L
)
Unit: seconds
expr min lq mean median uq max neval
mean - base R 33.152048 33.152048 33.152048 33.152048 33.152048 33.152048 1
mean - data.table 8.287379 8.287379 8.287379 8.287379 8.287379 8.287379 1
microbenchmark(
"rollmean - base R" = {
apply(X = data,
MARGIN = 1:2,
FUN = rolling_mean)
},
"rollmean - data.table" = {
DT[, rollmean := frollmean(value, n = 3), by = .(lat, lon)]
},
times = 1L
)
Unit: seconds
expr min lq mean median uq max neval
rollmean - base R 477.26644 477.26644 477.26644 477.26644 477.26644 477.26644 1
rollmean - data.table 48.80027 48.80027 48.80027 48.80027 48.80027 48.80027 1
如您所见，使用data.table可显着更快地提供结果。
为了获得数组的有效方法，我们可以转置数组并使用colSums(...). 这比（我最喜欢的）data.table更快，并且不需要从数组强制到类似 data.frame 的对象。感谢@ismirsehregal 提供的基准测试结构。
请注意，这是较小的数据集。虽然从数组 -> data.table 转换会产生转换时间损失，但该data.table方法在某些时候可能会变得更快，因为它可以利用与该方法不同的并行核心colMeans。我没有足够的 RAM 来使用更大的数据集。
colMeans(aperm(data, c(3L, 1L, 2L)))
all.equal(apply(data, 1:2, mean), colMeans(aperm(data, c(3L, 1L, 2L))))
## [1] TRUE
microbenchmark::microbenchmark(
"mean - base R" = apply(data, 1:2, mean),
"mean - data.table" = DT[, .("mean_value" = mean(value)), by = .(lat, lon)],
"mean - base colMeans" = colMeans(aperm(data, c(3L, 1L, 2L)))
)
## Unit: milliseconds
## expr min lq mean median uq max neval
## mean - base R 48.6293 52.30380 57.611552 54.8006 61.2880 146.4658 100
## mean - data.table 7.5610 8.95255 13.004569 9.7459 17.4888 28.7324 100
## mean - base colMeans 6.1713 6.45505 8.421273 6.6429 7.7849 99.3952 100
对于高效的滚动手段来说，zoo::roll_applyr()是直观但不快速。我们可以使用许多包解决方案，例如data.table::frollmean()或RcppRoll::roll_mean()。另请参阅此处了解滚动平均值的各种方法：
或者@JanGorecki 关于性能的问题：
data.table::frollmean()除了使用with之外，我还从第一个链接介绍了 @Matti Pastell 的解决方案和 @pipefish 的解决方案apply：
##@Matti Pastell's
base_ma = function(x, n = 5){stats::filter(x, rep(1 / n, n), sides = 2)}
##@pipefish's
base_ma2 = function(x, n = 5L) {
cx <- c(0,cumsum(x))
(cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n
}
microbenchmark(
"rollmean - apply zoo" = {
apply(X = data,
MARGIN = 1:2,
FUN = rolling_mean)
},
"rollmean - apply dt.frollmean" = {
apply(X = data,
MARGIN = 1:2,
FUN = frollmean, n = 3L)
},
"rollmean - data.table" = {
DT[, rollmean := frollmean(value, n = 3), by = .(lat, lon)]
},
"rollmean - apply stats::filter" = {
apply(X = data,
MARGIN = 1:2,
FUN = base_ma, n = 3L)
},
"rollmean - apply w_cumsum" = {
apply(X = data,
MARGIN = 1:2,
FUN = base_ma2, n = 3L)
},
times = 1L
)
Unit: milliseconds
expr min lq mean median uq max neval
rollmean - apply zoo 2823.6218 2823.6218 2823.6218 2823.6218 2823.6218 2823.6218 1
rollmean - apply dt.frollmean 331.7787 331.7787 331.7787 331.7787 331.7787 331.7787 1
rollmean - data.table 358.0787 358.0787 358.0787 358.0787 358.0787 358.0787 1
rollmean - apply stats::filter 435.5238 435.5238 435.5238 435.5238 435.5238 435.5238 1
rollmean - apply w_cumsum 148.7428 148.7428 148.7428 148.7428 148.7428 148.7428 1