R - 计算前 k 个非 NA 值的滚动平均值
R - Calculate rolling mean of previous k non-NA values
我正在尝试计算 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 应该是：
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
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，但除非我弄错了，否则它们似乎在固定滚动 window 以及处理 NA 值的选项（例如忽略它们）。就我而言，我想“扩展”我的 window 以包含 k 非 NA 值。
欢迎任何thoughts/suggestions。
编辑 - 模拟结果
感谢所有贡献者。首先，我没有提到我的数据集确实更大，而且 运行 通常，所以任何性能改进都是最受欢迎的。因此，我 运行 下面的模拟来检查执行时间，然后再决定接受哪个答案。请注意，某些答案需要对 return 所需的输出进行小的调整，但如果您觉得您的解决方案被歪曲（因此效率低于预期），请随时告诉我，我会相应地进行编辑.我在下面的回答中使用了 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 解决方案似乎更快。非常欢迎任何最终建议。
非常感谢您的宝贵时间！！
由于我不知道在任何标准库中有现成的计算输出的方法，我想出了下面的实现 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
rollapplyr. 关于问题中关于 rollmean 的评论，zoo 也有 rollappy 和 rollapplyr（右对齐），它们允许每个组件的不同宽度（和偏移量）通过指定向量（就像我们在这里所做的那样）或宽度列表来输入——请参阅 ?rollapply 了解更多信息。我们在下面使用了一个相对简单的宽度向量，还展示了一些改进的宽度向量，它们 运行 更快。
操作 创建一个 Mean 函数，它接受一个向量，删除最后一个元素和所有 NA，然后取剩下的最后 k 个元素，将其扩展为 k 个 NA 元素需要。最后取其平均值。我们使用 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))
将 rollapplyr 的第二个参数替换为 w，此处显示了 w。这里的想法是使用 findInterval 找到元素 k 非 NA 的背面，它提供了更严格的界限。这个（加上使用 na_rm）运行 在尝试使用相同的 1300 行时，以增加 2 行代码为代价，几乎是问题中代码的两倍。
tt <- length(x) - rev(cumsum(rev(!is.na(x))))
w <- seq_along(tt) - findInterval(tt - k - 1, tt)
Code. 使用问题中的数据，下面的代码（未使用上述改进）运行 比问题中的代码稍快（不是很多）基于我的基准测试的问题，它只有两行代码。
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))
给予：
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
不使用 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
使用 runner 它将类似于 mean 的 3 元素 tail window 的非 na 值。您可以使用 slider
获得相同的结果
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
我正在尝试计算 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 应该是：
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
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，但除非我弄错了，否则它们似乎在固定滚动 window 以及处理 NA 值的选项（例如忽略它们）。就我而言，我想“扩展”我的 window 以包含 k 非 NA 值。
欢迎任何thoughts/suggestions。
编辑 - 模拟结果
感谢所有贡献者。首先，我没有提到我的数据集确实更大，而且 运行 通常，所以任何性能改进都是最受欢迎的。因此，我 运行 下面的模拟来检查执行时间，然后再决定接受哪个答案。请注意，某些答案需要对 return 所需的输出进行小的调整，但如果您觉得您的解决方案被歪曲（因此效率低于预期），请随时告诉我，我会相应地进行编辑.我在下面的回答中使用了 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 解决方案似乎更快。非常欢迎任何最终建议。
非常感谢您的宝贵时间！！
由于我不知道在任何标准库中有现成的计算输出的方法，我想出了下面的实现 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
rollapplyr. 关于问题中关于 rollmean 的评论，zoo 也有 rollappy 和 rollapplyr（右对齐），它们允许每个组件的不同宽度（和偏移量）通过指定向量（就像我们在这里所做的那样）或宽度列表来输入——请参阅 ?rollapply 了解更多信息。我们在下面使用了一个相对简单的宽度向量，还展示了一些改进的宽度向量，它们 运行 更快。
操作 创建一个 Mean 函数，它接受一个向量，删除最后一个元素和所有 NA，然后取剩下的最后 k 个元素，将其扩展为 k 个 NA 元素需要。最后取其平均值。我们使用 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))将 rollapplyr 的第二个参数替换为 w，此处显示了 w。这里的想法是使用 findInterval 找到元素 k 非 NA 的背面，它提供了更严格的界限。这个（加上使用 na_rm）运行 在尝试使用相同的 1300 行时，以增加 2 行代码为代价，几乎是问题中代码的两倍。
tt <- length(x) - rev(cumsum(rev(!is.na(x)))) w <- seq_along(tt) - findInterval(tt - k - 1, tt)
Code. 使用问题中的数据，下面的代码（未使用上述改进）运行 比问题中的代码稍快（不是很多）基于我的基准测试的问题，它只有两行代码。
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))
给予：
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
不使用 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
使用 runner 它将类似于 mean 的 3 元素 tail window 的非 na 值。您可以使用 slider
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