R - 计算前 k 个非 NA 值的滚动平均值答案

【问题标题】：R - Calculate rolling mean of previous k non-NA valuesR - 计算前 k 个非 NA 值的滚动平均值
【发布时间】：2021-07-03 15:58:27
【问题描述】：

我正在尝试计算 dplyr/tidyverse 框架内之前 k 个非 NA 值的滚动平均值。我编写了一个似乎可以工作的函数，但想知道是否已经有某个包中的函数（这可能比我的尝试更有效率）正在执行此操作。示例数据集：

tmp.df <- data.frame(
  x = c(NA, 1, 2, NA, 3, 4, 5, NA, NA, NA, 6, 7, NA)
)

假设我想要前 3 个非 NA 值的滚动平均值。那么输出y应该是：

y 的前 5 个元素是 NAs，因为 x 第一次有 3 个之前的非 NA 值位于第 6 行，这 3 个元素的平均值为 2。下一个 y 元素是不言自明的。第 9 行得到 4，因为 x 的前 3 个非 NA 值位于第 5、6 和 7 行，依此类推。

我的尝试是这样的：

roll_mean_previous_k <- function(x, k){
  
  require(dplyr)
  
  res                      <- NA
  lagged_vector            <- dplyr::lag(x)
  lagged_vector_without_na <- lagged_vector[!is.na(lagged_vector)]
  previous_k_values        <- tail(lagged_vector_without_na, k)
  
  if (length(previous_k_values) >= k) res <- mean(previous_k_values)
  
  res
  
}

按如下方式使用（使用slider 包中的slide_dbl 函数）：

library(dplyr)

tmp.df %>% 
  mutate(
    y = slider::slide_dbl(x, roll_mean_previous_k, k = 3, .before = Inf)
  )

提供所需的输出。但是，我想知道是否有现成的（如前所述）更有效的方法来做到这一点。我应该提到，我分别从zoo 和RcppRoll 包中知道rollmean 和roll_mean，但除非我弄错了，否则它们似乎在固定滚动窗口上工作，可以选择处理@987654338 @ 值（例如忽略它们）。就我而言，我想“扩展”我的窗口以包含 k 非 NA 值。

欢迎提出任何想法/建议。

编辑 - 模拟结果

感谢所有贡献者。首先，我没有提到我的数据集确实更大并且经常运行，因此任何性能改进都是最受欢迎的。因此，在决定接受哪个答案之前，我运行了以下模拟来检查执行时间。请注意，某些答案需要稍作调整才能返回所需的输出，但如果您认为您的解决方案被歪曲（因此效率低于预期），请随时告诉我，我会相应地进行编辑。我在下面的回答中使用了G. Grothendieck 的技巧，以消除对if-else 检查滞后、非NA 向量长度的需要。

所以这是模拟代码：

library(tidyverse)
library(runner)
library(zoo)
library(slider)
library(purrr)
library(microbenchmark)

set.seed(20211004)
test_vector <- sample(x = 100, size = 1000, replace = TRUE)
test_vector[sample(1000, size = 250)] <- NA

# Based on GoGonzo's answer and the runner package
f_runner <- function(z, k){
  
  runner(
    x = z, 
    f = function(x) {
      mean(`length<-`(tail(na.omit(head(x, -1)), k), k)) 
    }
  )
  
}

# Based on my inital answer (but simplified), also mentioned by GoGonzo 
f_slider <- function(z, k){
  
  slide_dbl(
    z,
    function(x) {
      mean(`length<-`(tail(na.omit(head(x, -1)), k), k)) 
    },
    .before = Inf
  )
}

# Based on helios' answer. Return the correct results but with a warning.
f_helios <- function(z, k){
  
    reduced_vec <-  na.omit(z)
    unique_means <-  rollapply(reduced_vec, width = k, mean)
    
    start <-  which(!is.na(z))[k] + 1
    repeater <-  which(is.na(z)) + 1
    repeater_cut <-  repeater[(repeater > start-1) & (repeater <= length(z))]
    
    final <- as.numeric(rep(NA, length(z)))
    index <-  start:length(z)
    final[setdiff(index, repeater_cut)] <- unique_means
    final[(start):length(final)] <- na.locf(final)
    final
}

# Based on G. Grothendieck's answer (but I couldn't get it to run with the performance improvements)
f_zoo <- function(z, k){
  
  rollapplyr(
    z, 
    seq_along(z), 
    function(x, k){
      mean(`length<-`(tail(na.omit(head(x, -1)), k), k)) 
    },
    k)

}

# Based on AnilGoyal's answer
f_purrr <- function(z, k){
  
    map_dbl(
      seq_along(z), 
      ~ ifelse(
        length(tail(na.omit(z[1:(.x -1)]), k)) == k,
        mean(tail(na.omit(z[1:(.x -1)]), k)), 
        NA
        )
      )

}

# Check if all are identical #
all(
  sapply(
    list(
      # f_helios(test_vector, 10),
      f_purrr(test_vector, 10),
      f_runner(test_vector, 10),
      f_zoo(test_vector, 10)
    ), 
    FUN = identical, 
    f_slider(test_vector, 10),
  )
)

# Run benchmarking #
microbenchmark(
  # f_helios(test_vector, 10),
  f_purrr(test_vector, 10),
  f_runner(test_vector, 10),
  f_slider(test_vector, 10),
  f_zoo(test_vector, 10)
)

结果：

Unit: milliseconds
                      expr     min       lq     mean   median       uq      max neval  cld
  f_purrr(test_vector, 10) 31.9377 37.79045 39.64343 38.53030 39.65085 104.9613   100   c 
 f_runner(test_vector, 10) 23.7419 24.25170 29.12785 29.23515 30.32485  98.7239   100  b  
 f_slider(test_vector, 10) 20.6797 21.71945 24.93189 26.52460 27.67250  32.1847   100 a   
    f_zoo(test_vector, 10) 43.4041 48.95725 52.64707 49.59475 50.75450 122.0793   100    d

基于上述，除非代码可以进一步改进，否则slider 和runner 解决方案似乎更快。任何最终建议都非常受欢迎。

非常感谢您的宝贵时间！！

【问题讨论】：

标签： r dplyr na rolling-computation

【解决方案1】：

使用runner，它将类似于mean 的3 元素tail 非na 值的窗口。您可以使用滑块获得相同的结果

library(runner)
tmp.df <- data.frame(
  x = c(NA, 1, 2, NA, 3, 4, 5, NA, NA, NA, 6, 7, NA)
)

# using runner
tmp.df$y_runner <- runner(
  x = tmp.df$x, 
  f = function(x) {
    mean(
      tail(
        x[!is.na(x)],
        3
      )
    )
  }
)

# using slider
tmp.df$y_slider <- slider::slide_dbl(
  tmp.df$x, 
  function(x) {
    mean(
      tail(
        x[!is.na(x)],
        3
      )
    )
  }, 
  .before = Inf
)

tmp.df


#    x    y_runner y_slider
# 1  NA      NaN      NaN
# 2   1      1.0      1.0
# 3   2      1.5      1.5
# 4  NA      1.5      1.5
# 5   3      2.0      2.0
# 6   4      3.0      3.0
# 7   5      4.0      4.0
# 8  NA      4.0      4.0
# 9  NA      4.0      4.0
# 10 NA      4.0      4.0
# 11  6      5.0      5.0
# 12  7      6.0      6.0
# 13 NA      6.0      6.0

【讨论】：

您的包裹runner() 是非常好的GoGonzo。 +1
非常感谢。请注意，runner::runner() 目前可能不起作用，请改用library(runner);runner(...)。这将很快在下一个版本中修复
感谢您的回复。 runner 包看起来确实很有趣，尽管你的输出并不是我想要的（x 没有滞后，前 5 个元素应该是 NAs）。我已经修改了runner 中的函数来完成它，它似乎可以工作。我将运行一个模拟来测试执行时间。谢谢。

【解决方案2】：

rollapplyr. 关于问题中关于 rollmean 的评论，zoo 也有 rollappy 和 rollapplyr（右对齐），它们通过指定向量允许输入的每个组件具有不同的宽度（和偏移量）（就像我们在这里所做的那样）或宽度列表——请参阅 ?rollapply 了解更多信息。我们在下面使用了一个相对简单的宽度向量，并展示了一些改进的宽度向量，它们运行得更快。

操作创建一个 Mean 函数，该函数接受一个向量，删除最后一个元素和所有 NA，并根据需要将剩下的最后 k 个元素扩展到具有 NA 的 k 个元素。最后取其平均值。我们使用 rollapplyr 将其应用于宽度为 seq_along(x) 的 x。

性能改进

用折叠包中的 na_rm 替换 na.omit
将 rollapplyr 的第二个参数替换为此处显示的代码。这里的想法是，NA 的 k+1 个最长游程的长度之和加上 k+1 形成了我们需要考虑的元素数量的界限。当我尝试使用 1300 行（由问题中的 100 个数据副本形成）并且没有添加太多额外代码时，这（加上使用 na_rm）的运行速度比问题中的代码快 25%。
```
pmin(with(rle(is.na(x)), sum(tail(sort(lengths[values]), k+1)))+k+1, seq_along(x))
```
用 w 替换 rollapplyr 的第二个参数，此处显示 w。这里的想法是使用 findInterval 找到元素 k 非 NA 的背面，这提供了更紧密的界限。当尝试使用相同的 1300 行以增加 2 行代码为代价时，这个（加上使用 na_rm）的运行速度几乎是问题中代码的两倍。
```
tt <- length(x) - rev(cumsum(rev(!is.na(x))))
w <- seq_along(tt) - findInterval(tt - k - 1, tt)
```

代码。使用问题中的数据，根据我的基准测试，下面的代码（不使用上述改进）比问题中的代码运行得稍快（不是很多），它只是两行代码。

library(dplyr)
library(zoo)

Mean <- function(x, k) mean(`length<-`(tail(na.omit(head(x, -1)), k), k))
tmp.df %>% mutate(y = rollapplyr(x, seq_along(x), Mean, k = 3))

给予：

【讨论】：

感谢您的回答，尤其是length<- 的技巧。但是，我无法应用您建议的性能改进。我已经发布了一些模拟结果，如果您认为您的建议会进一步减少执行时间，如果您能相应地修改代码，我将不胜感激。非常感谢您的时间和回答！

【解决方案3】：

由于我不知道在任何标准库中计算输出的现成方法，我想出了下面的实现roll_mean_k_efficient，这似乎大大加快了你的计算速度。请注意，此实现使用了 zoo 包中的 rollapply 和 na.locf 方法。

rm(list = ls())

library("zoo")
library("rbenchmark")
library("dplyr")

x = rep(c(NA, 1, 2, NA, 3, 4, 5, NA, NA, NA, 6, 7, NA), 100)

# your sample (extended)
tmp.df <- data.frame(
  x = rep(c(NA, 1, 2, NA, 3, 4, 5, NA, NA, NA, 6, 7, NA), 100)
)

# enhanced implementation
roll_mean_k_efficient <- function(x, k){
  reduced_vec = na.omit(x)
  unique_means = rollapply(reduced_vec, width=k, mean)
  
  start = which(!is.na(x))[k] + 1
  repeater = which(is.na(x)) + 1
  repeater_cut = repeater[(repeater > start-1) & (repeater <= length(x))]
  
  final <- as.numeric(rep(NA, length(x)))
  index = start:length(x)
  final[setdiff(index, repeater_cut)] <- unique_means
  final[(start):length(final)] <- na.locf(final)
  final
}

# old implementation
roll_mean_previous_k <- function(x, k){
  res                      <- NA
  lagged_vector            <- dplyr::lag(x)
  lagged_vector_without_na <- lagged_vector[!is.na(lagged_vector)]
  previous_k_values        <- tail(lagged_vector_without_na, k)
  if (length(previous_k_values) >= k) res <- mean(previous_k_values)
  res
}

# wrapper function for the benchmarking below
roll_mean_benchmark = function(){
  res = tmp.df %>% 
    mutate(
      y = slider::slide_dbl(x, roll_mean_previous_k, k = 3, .before = Inf)
    ) 
  return(res)
}

# some benchmarking
benchmark(roll_mean_k_efficient(x = x, k=3), 
          roll_mean_benchmark(), 
          columns=c('test','elapsed','replications'),
          replications = 100)

此外，我扩展了您的示例向量x，以通过rbenchmark 包中的benchmark 函数获得一些更可靠的基准测试结果。在我的情况下，运行代码后打印的基准运行时是：

                                 test elapsed replications
2               roll_mean_benchmark()   4.463          100
1 roll_mean_k_efficient(x = x, k = 3)   0.039          100

【讨论】：

感谢您的回答。但是，我在获得所需的输出时遇到了一些问题，并且无法确定问题所在。如果需要，您可以查看我发布的模拟结果以复制警告和不同的结果。但绝对感谢您的宝贵时间！
哦，真的吗？该函数在原始问题中给出的示例数据集上没有按预期工作吗？在这种情况下，在我的机器上 - 返回您问题的确切结果。你能把你的模拟测试用例提供给我，然后我可以看看哪里出了问题。
对不起，我应该说它返回以下警告：“警告消息：在 final[setdiff(index, repeater_cut)]

【解决方案4】：

不使用zoo。在tidyverse 时尚中，您也可以使用purrr::map 进行操作


tmp.df %>% mutate(y = map(seq_along(x), ~ ifelse(length(tail(na.omit(tmp.df$x[1:(.x -1)]), 3)) ==3, 
                                                 mean(tail(na.omit(tmp.df$x[1:(.x -1)]), 3)), 
                                                 NA)))

    x  y
1  NA NA
2   1 NA
3   2 NA
4  NA NA
5   3 NA
6   4  2
7   5  3
8  NA  4
9  NA  4
10 NA  4
11  6  4
12  7  5
13 NA  6

【讨论】：