滚动将数据表中的多个分位数应用于多个列

Question

背景：
我可以使用 data.table（见附件）从我的数据中获取多个时刻，但这需要很长时间。我在想，对表格进行排序以获得特定百分位数的过程对于找到几个百分位数会更有效。

像中位数这样的一次性统计数据需要 1.79 毫秒，而非中位数分位数在 122.8 毫秒时需要 68 倍的时间。必须有一种方法来减少计算时间。

问题：

• 有没有一种方法可以更有效地从同一数据中调用多个分位数？
• 我可以将“lapply”从 data.table 中提取出来，然后像我做名单一样组成它吗？

我的示例代码包含少量合成数据：

#libraries
library(data.table)      #data.table
library(zoo)             #roll apply

#reproducibility
set.seed(45L)

#make data
DT<-data.table(V1=c(1L,2L),
               V2=LETTERS[1:3],
               V3=round(rnorm(300),4),
               V4=round(runif(150),4),
               V5=1:1200)
DT

#get names
my_col_list <- names(DT)[c(3,4)]

#make new variable names
my_name_list1 <- paste0(my_col_list, "_", "33rd_pctile")
my_name_list2 <- paste0(my_col_list, "_", "77rd_pctile")

#compute values
for(i in 1:length(my_col_list)){
     #first 
     DT[, (my_name_list1[i]) := unlist(lapply(.SD,
                                          function(x) rollapply(x,
                                                                7,
                                                                quantile,
                                                                fill = NA,
                                                                probs = 1/3)), 
                                   recursive = F),
        .SDcols = my_col_list[i]]
     #second
     DT[, (my_name_list2[i]) := unlist(lapply(.SD,
                                          function(x) rollapply(x,
                                                                7,
                                                                quantile,
                                                                fill = NA,
                                                                probs = 7/9)), 
                                   recursive = F),
        .SDcols = my_col_list[i]]
}

#display it
head(DT,10)

对一次性统计数据与多次统计数据进行微基准测试表明，分位数很昂贵。

res2 <- microbenchmark(          DT[, (my_name_list1[i]) := unlist(lapply(.SD,
                                                                          function(x) rollapply(x,
                                                                                                7,
                                                                                                mean,
                                                                                                fill = NA)), 
                                                                   recursive = F),
                                    .SDcols = my_col_list[i]],
                                 times = 5)

表示平均需要大约 1.75 毫秒（中位数为 1.79 秒）

> res2
Unit: milliseconds
                                                                                                                                            expr
 DT[, `:=`((my_name_list1[i]), unlist(lapply(.SD, function(x) rollapply(x,      7, mean, fill = NA)), recursive = F)), .SDcols = my_col_list[i]]
      min       lq     mean   median       uq     max neval
 1.465779 1.509114 1.754145 1.618591 1.712103 2.46514     5

但计算分位数需要 100 倍

res3 <- microbenchmark(          DT[, (my_name_list1[i]) := unlist(lapply(.SD,
                                                                          function(x) rollapply(x,
                                                                                                7,
                                                                                                quantile,
                                                                                                fill = NA,
                                                                                                probs = 1/3)), 
                                                                   recursive = F),
                                    .SDcols = my_col_list[i]],
                                 times = 5)

res3

和

> res3
Unit: milliseconds
                                                                                                                                                             expr
 DT[, `:=`((my_name_list1[i]), unlist(lapply(.SD, function(x) rollapply(x,      7, quantile, fill = NA, probs = 1/3)), recursive = F)), .SDcols = my_col_list[i]]
      min       lq     mean   median       uq      max neval
 118.5833 119.2896 122.8432 124.0168 124.4183 127.9082     5

更新：

• “分位数”的中位数大约需要 128 秒，而“中位数” 花费少得多。它们不是一回事。
• 遍历“分位数”的“类型”的 9 个选项给出平均值 时间在 129 毫秒和 157 毫秒之间。这里没有“易赢”。
• “WGCNA”包需要来自 bioconductor 的“GO.db”，这不是 使用“install.packages”命令安装。还需要包 未与“WGCNA”或“GO.db”一起安装的“impute”。还 “预处理核心”。
• 使用（最终可以工作的）WGCNA 包减少了
  的平均时间 滚动分位数到 68.1 毫秒。这是大约一半的时间，但它是 不到 1/70 的时间。
• 使用“RollingWindow”包中的“RollingMedian”得到 169.8 微秒（又名 0.1698 毫秒），这快了很多，但是 不是任意分位数。
• 使用“perccal”似乎将分位数的计算降低到 ~145
  微秒。在 rollapply 中，这会将计算时间降低到 15.3 毫秒，这是一个 8.5 倍的提升。我不确定还有多少 这块石头里有血要挤出来。

想法：

• “perccal”方法似乎只使用单个内核。那里 可能是一些“并行”过程，可以让我拆分摘要 针对不同核心的不同变量。这可能会给一些 加速。
• 随着我们向数据中添加更多项，加速开始降低。 增加到 9600 行可将加速从约 8.5 倍降低到 1 倍以下。 这表明 rollapply 函数也可能存在一些问题。

Answer 1

数据表优化

你说得对，中位数在这种情况下特别快，这是因为它运行的是专门的 C 代码，而不是纯 R 代码的分位数函数。

我们可以在data.table 的文档中阅读到这种优化

?data.table.optimize

我们有：

当 j 中的表达式仅包含这些函数 min、max、 平均值、中位数、var、sd、prod，例如 dt[、list(mean(x)、 中位数(x), min(y), max(y)), by=z]，它们得到了非常有效的优化 使用我们所说的 GForce。这些函数被 gmean 替换， gmedian, gmin, gmax 代替

他们举了一个例子来说明中位数情况下的速度提升：

# Generate a big data.table with a relatively many columns
set.seed(1L)
dt = lapply(1:20, function(x) sample(c(-100:100), 5e6L, TRUE))
setDT(dt)[, id := sample(1e5, 5e6, TRUE)]
print(object.size(dt), units="Mb") # 400MB, not huge, but will do

# GForce
options(datatable.optimize = 2L) # optimisation 'on'
system.time(ans1 <- dt[, lapply(.SD, median), by=id])
system.time(ans2 <- dt[, lapply(.SD, function(x) as.numeric(stats::median(x))), by=id])
identical(ans1, ans2)

在我的系统上，R 内部版本比 data.table 版本慢大约 44 倍。

加速分位数

我们仍然可以尝试提高 R 中 quantile 函数的速度，为此我的方法基本上是“使用源代码，Luke”并查看分位数函数。查看源代码，我们得到标准的泛型函数：

>> quantile
function (x, ...) 
UseMethod("quantile")
<bytecode: 0x0000000009154c78>
<environment: namespace:stats>

我们可以进一步追踪：

>> methods(quantile)
[1] quantile.default* quantile.ecdf*    quantile.POSIXt*  quantile.zoo     
see '?methods' for accessing help and source code

并查看默认功能。

>> stats:::quantile.default
function (x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE, 
    type = 7, ...) 
{
...
}

现在我们有了整个源，它很长，我们可以将它与median.default 中的 R 中值源进行比较。使用源代码，我们可以将其复制为用户定义的函数并对其进行分析（包括为format_perc 提供命名空间的一小部分内容），从中可以清楚地看出只有两行相关，即排序和输出格式，排序与中值函数非常相似，可能很难改进。但是，可以通过将其注释掉来完全跳过格式。

fast.quant <- function (x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE, 
                        type = 7, ...) 
{
  if (is.factor(x)) {
    ...

    ...
    if (names && np > 0L) {
      #names(qs) <- stats:::format_perc(probs)
    }
    ...
}

总而言之，这个修复将运行时间减少了一半，它仍然不是优化的中值，但很可能在不离开 R 的情况下很难获得更好的性能。

有可能，甚至很可能，data.table 中的优化也可以用来帮助进行分位数计算，因为 data.table 也在 C 中实现了排序。然而，人们仍然希望利用仅需要部分排序的优势。否则，Rcpp 包也可以用于执行类似的优化。