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

Question

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

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

问题：

• 有没有办法以更有效的方式从同一数据中调用多个分位数？
• 我可以从 data.table 中提取“lapply”并像我做名字列表一样编写它吗？

我的带有微小合成数据的示例代码：

#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.5x 降低到 1x 以下。这表明 rollapply 函数也可能存在一些问题。

Answer 1

数据表优化

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

data.table我们可以在文档中阅读有关此优化的信息

?data.table.optimize

在那里我们有：

当 j 中的表达式仅包含这些函数 min、max、mean、median、var、sd、prod，例如 dt[、list(mean(x)、median(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 倍。

加速分位数

我们仍然可以尝试提高quantileR 中函数的速度，为此我的方法基本上是“使用源代码，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, ...) 
{
...
}

现在我们有了整个源代码，它很长，我们可以将它与 中的 R 中值源进行比较median.default。有了源代码，我们可以将其复制为用户定义的函数并对其进行分析（包含一小部分提供 的命名空间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否则，也可以使用 该包来执行类似的优化。