R：低于可变分位数阈值的值的滚动平均值答案

【问题标题】：R: Rolling mean on values below variable quantile thresholdR：低于可变分位数阈值的值的滚动平均值
【发布时间】：2020-03-16 11:23:12
【问题描述】：

我有一个很大的data.table。我需要：

数值列的滚动分位数
应用于低于（移动）分位数阈值的数值列的滚动平均值

library(data.table)

set.seed(101)
data <- data.table(group=c(rep("A",10),rep("B",7)), value=rnorm(17))

> data
    group      value
 1:     A -0.3260365
 2:     A  0.5524619
 3:     A -0.6749438
 4:     A  0.2143595
 5:     A  0.3107692
 6:     A  1.1739663
 7:     A  0.6187899
 8:     A -0.1127343
 9:     A  0.9170283
10:     A -0.2232594
11:     B  0.5264481
12:     B -0.7948444
13:     B  1.4277555
14:     B -1.4668197
15:     B -0.2366834
16:     B -0.1933380
17:     B -0.8497547

为了确定“值”列的滚动分位数，我使用了 caTools 包中的函数 runquantile()，这需要大约 1 分钟才能在我的数据集上执行。

对于相同的滚动窗口（此处为k=4），当计算时间受到关注时，如何获得移动分位数以下的滚动平均值？在示例中，结果应该类似于列“mean_below_q”。

library(caTools)

data[,rolling_q := c(rep(NA,3),runquantile(value,k=4,0.4,endrule="trim")),group]

> data
    group      value   rolling_q mean_below_q
 1:     A -0.3260365          NA           NA
 2:     A  0.5524619          NA           NA
 3:     A -0.6749438          NA           NA
 4:     A  0.2143595 -0.21795730  -0.50049020
 5:     A  0.3107692  0.23364141  -0.23029220
 6:     A  1.1739663  0.23364141  -0.23029220
 7:     A  0.6187899  0.37237334   0.26256430
 8:     A -0.1127343  0.37237334   0.09901745
 9:     A  0.9170283  0.67843754   0.25302780
10:     A -0.2232594  0.03357052  -0.16799680
11:     B  0.5264481          NA           NA
12:     B -0.7948444          NA           NA
13:     B  1.4277555          NA           NA
14:     B -1.4668197 -0.53058593  -1.13083200
15:     B -0.2366834 -0.68321222  -1.13083200
16:     B -0.1933380 -0.22801430  -0.85175150
17:     B -0.8497547 -0.72714047  -1.15828700

非常感谢。

根据@chinsoon12 出色的脚本，它解决了这个问题，并且在大型数据集上大约需要 15 分钟。

# Option 1)

library(data.table)

set.seed(101)
data <- data.table(group=c(rep("A",10),rep("B",7),rep("C",13)), value=c(rnorm(17),rep(0,3),-0.9,rep(0,4),-0.5,-0.8,rep(0,3)))
sz <- 4L
prob <- 0.4

#see stats:::quantile.default
index <- 1 + (sz - 1) * prob
lo <- floor(index)
hi <- ceiling(index)
data[, c("sr", "er") := .(.I - sz + 1L, .I)]

data[, rolling_q := frollapply(value, sz, function(x) {
        x <- sort(x, partial = unique(c(lo, hi)))
        qs <- x[lo]
        i <- which(index > lo & x[hi] != qs)
        h <- (index - lo)[i]
        qs[i] <- (1 - h) * qs[i] + h * x[hi[i]]
        qs
        #.(q, sum(x[1L:lo]) / lo)
    }), group]

data[, mean_below_q := 
    data[data, on=.(group, er>=sr, er<=er), allow.cartesian=TRUE,
        by=.EACHI, mean(value[value<=i.rolling_q])]$V1
]

最后，我尝试将 caTools 包中的 runquantile() 与选项 1) 代码结合使用，用于滚动均值低于阈值。它似乎运行良好，在大型数据集上总共需要 2.5 分钟。

# rolling_q via runquantile & mean_below_q via Option 1)

library(data.table)
library(caTools)

set.seed(101)
data <- data.table(group=c(rep("A",10),rep("B",7),rep("C",13)), value=c(rnorm(17),rep(0,3),-0.9,rep(0,4),-0.5,-0.8,rep(0,3)))

data[,rolling_q := c(rep(NA,3),runquantile(value,k=4,0.4,endrule="trim")),group]

sz <- 4L
prob <- 0.4

index <- 1 + (sz - 1) * prob
lo <- floor(index)
hi <- ceiling(index)
data[, c("sr", "er") := .(.I - sz + 1L, .I)]

data[, mean_below_q := 
    data[data, on=.(group, er>=sr, er<=er), allow.cartesian=TRUE,
        by=.EACHI, mean(value[value<=i.rolling_q])]$V1
]

> data
    group      value   rolling_q sr er mean_below_q
 1:     A -0.3260365          NA -2  1           NA
 2:     A  0.5524619          NA -1  2           NA
 3:     A -0.6749438          NA  0  3           NA
 4:     A  0.2143595 -0.21795730  1  4  -0.50049017
 5:     A  0.3107692  0.23364141  2  5  -0.23029219
 6:     A  1.1739663  0.23364141  3  6  -0.23029219
 7:     A  0.6187899  0.37237334  4  7   0.26256434
 8:     A -0.1127343  0.37237334  5  8   0.09901745
 9:     A  0.9170283  0.67843754  6  9   0.25302777
10:     A -0.2232594  0.03357052  7 10  -0.16799684
11:     B  0.5264481          NA  8 11           NA
12:     B -0.7948444          NA  9 12           NA
13:     B  1.4277555          NA 10 13           NA
14:     B -1.4668197 -0.53058593 11 14  -1.13083206
15:     B -0.2366834 -0.68321222 12 15  -1.13083206
16:     B -0.1933380 -0.22801430 13 16  -0.85175154
17:     B -0.8497547 -0.72714047 14 17  -1.15828722
18:     C  0.0000000          NA 15 18           NA
19:     C  0.0000000          NA 16 19           NA
20:     C  0.0000000          NA 17 20           NA
21:     C -0.9000000  0.00000000 18 21  -0.22500000
22:     C  0.0000000  0.00000000 19 22  -0.22500000
23:     C  0.0000000  0.00000000 20 23  -0.22500000
24:     C  0.0000000  0.00000000 21 24  -0.22500000
25:     C  0.0000000  0.00000000 22 25   0.00000000
26:     C -0.5000000  0.00000000 23 26  -0.12500000
27:     C -0.8000000 -0.40000000 24 27  -0.65000000
28:     C  0.0000000 -0.40000000 25 28  -0.65000000
29:     C  0.0000000 -0.40000000 26 29  -0.65000000
30:     C  0.0000000  0.00000000 27 30  -0.20000000
    group      value   rolling_q sr er mean_below_q

【问题讨论】：

标签： r data.table mean quantile rolling-computation

【解决方案1】：

我认为有一个错字，因为数据显示的平均值低于百分位而不是高于百分位。

这里有 2 个选项：

1) 使用 frollapply（因为它返回长度为 1 的值），因此您需要另一个非等值连接来计算低于百分位的值的平均值。

data[, rolling_q := frollapply(value, sz, function(x) {
        x <- sort(x, partial = unique(c(lo, hi)))
        qs <- x[lo]
        i <- which(index > lo & x[hi] != qs)
        h <- (index - lo)[i]
        qs[i] <- (1 - h) * qs[i] + h * x[hi[i]]
        qs
        #.(q, sum(x[1L:lo]) / lo)
    }), group]

data[, mean_below_q := 
    data[data, on=.(group, er>=sr, er<=er), allow.cartesian=TRUE,
        by=.EACHI, mean(value[value<=i.rolling_q])]$V1
]

2) 使用 1 个非等值连接计算两个值

data[, c("rolling_q", "mean_below_q") := 
    .SD[.SD, on=.(group, er>=sr, er<=er), allow.cartesian=TRUE,
        by=.EACHI, {
            if (.N >= sz) {
                x <- x.value
                x <- sort(x, partial = unique(c(lo, hi)))
                qs <- x[lo]
                i <- which(index > lo & x[hi] != qs)
                h <- (index - lo)[i]
                qs[i] <- (1 - h) * qs[i] + h * x[hi[i]]
                qs

                .(qs, sum(x[1L:lo]) / lo)
            } else 
                .(NA_real_, NA_real_)
        }][, (1L:3L) := NULL]
]

输出：

    group      value sr er   rolling_q mean_below_q
 1:     A -0.3260365 -2  1          NA           NA
 2:     A  0.5524619 -1  2          NA           NA
 3:     A -0.6749438  0  3          NA           NA
 4:     A  0.2143595  1  4 -0.21795730  -0.50049017
 5:     A  0.3107692  2  5  0.23364141  -0.23029219
 6:     A  1.1739663  3  6  0.23364141  -0.23029219
 7:     A  0.6187899  4  7  0.37237334   0.26256434
 8:     A -0.1127343  5  8  0.37237334   0.09901745
 9:     A  0.9170283  6  9  0.67843754   0.25302777
10:     A -0.2232594  7 10  0.03357052  -0.16799684
11:     B  0.5264481  8 11          NA           NA
12:     B -0.7948444  9 12          NA           NA
13:     B  1.4277555 10 13          NA           NA
14:     B -1.4668197 11 14 -0.53058593  -1.13083206
15:     B -0.2366834 12 15 -0.68321222  -1.13083206
16:     B -0.1933380 13 16 -0.22801430  -0.85175154
17:     B -0.8497547 14 17 -0.72714047  -1.15828722

数据：

library(data.table)

set.seed(101)
data <- data.table(group=c(rep("A",10),rep("B",7)), value=rnorm(17))
sz <- 4L
prob <- 0.4

#see stats:::quantile.default
index <- 1 + (sz - 1) * prob
lo <- floor(index)
hi <- ceiling(index)
data[, c("sr", "er") := .(.I - sz + 1L, .I)]

【讨论】：

感谢您的回答。不幸的是，在更广泛的数据集上，我发现了上述代码的一些意外行为。我已经相应地编辑了我的问题。我实在想不通哪里出了问题。
另外，我确实犯了一个错误，因为我的意思是“低于”分位数阈值。感谢您指出这一点。
@Chris437 我已经修复了有多个零的情况。让我们关注选项 1，因为选项 2 太慢了。你能检查一下选项 1 是否还有其他错误吗？
感谢您的更新。选项 1) 现在正在运行整个数据集。但是，当应用于我在上面问题中给出的最后一个示例时，它需要 14.5 分钟并且仍然会为 mean_below_q 产生错误的值。我认为对rolling_q（60 秒）使用非常快的runquantile 函数并通过您的脚本确定 mean_below_q 可能是有意义的。但是，仍然存在一个错误（参见最后一个示例，将prob 更改为0.2）？
我没有测试包