计算移动平均线
Posted
技术标签:
【中文标题】计算移动平均线【英文标题】：Calculating moving average 【发布时间】：2010-10-19 03:10:49 【问题描述】：我正在尝试使用 R 来计算矩阵中一系列值的移动平均值。 R 中似乎没有built-in function 可以让我。是否有任何套餐提供一个？还是我需要自己写？
【问题讨论】：
【参考方案1】：或者你可以简单地使用过滤器计算它，这是我使用的函数：
ma <- function(x, n = 5)filter(x, rep(1 / n, n), sides = 2)
如果你使用dplyr，请注意在上面的函数中指定stats::filter。
【讨论】：
我应该指出，“sides=2”可能是许多人不想忽视的用例中的一个重要选项。如果您只希望移动平均线中的尾随信息，则应使用边数=1。 几年后，dplyr 现在有了过滤功能，如果你加载了这个包，请使用stats::filter
sides = 2 相当于 zoo::rollmean 或 RcppRoll::roll_mean 的 align="center"。 sides = 1 相当于“右”对齐。我看不到“左”对齐或使用“部分”数据（2 个或更多值）计算的方法？
stats::filter 给出一个时间序列对象。将结果传递给as.vector 以获取向量。【参考方案2】：
zoo 包中的滚动均值/最大值/中位数 (rollmean)
TTR 中的移动平均值
马在forecast
【讨论】：
R 中不包含给定时间戳的未来值的移动平均线是什么？我检查了forecast::ma，它包含所有社区，不对。
尝试使用stats::filter 函数。在那里，您可以设置 sides = 1 仅用于过去的值。例如。 stats::filter(x, rep(1,5), sides = 1)/5 表示超过 5 个值的平均值。【参考方案3】：
使用cumsum 应该足够有效。假设你有一个向量 x 并且你想要 n 个数字的运行总和
cx <- c(0,cumsum(x))
rsum <- (cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n
正如@mzuther 在 cmets 中指出的那样，这假设数据中没有 NA。处理这些需要将每个窗口除以非 NA 值的数量。这是一种方法，结合了@Ricardo Cruz 的评论：
cx <- c(0, cumsum(ifelse(is.na(x), 0, x)))
cn <- c(0, cumsum(ifelse(is.na(x), 0, 1)))
rx <- cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]
rn <- cn[(n+1):length(cx)] - cn[1:(length(cx) - n)]
rsum <- rx / rn
这仍然存在一个问题，如果窗口中的所有值都是 NA，那么就会出现除以零的错误。
【讨论】：
此解决方案的一个缺点是它无法处理缺失：cumsum(c(1:3,NA,1:3))
@Ricardo Cruz：删除 NA 并相应地调整向量长度可能会更好。想想一个有很多 NA 的向量——零会将平均值拉向零，而移除 NA 会使平均值保持原样。当然，这完全取决于您的数据和您想要回答的问题。 :)
@mzuther，我在您的 cmets 之后更新了答案。感谢您的输入。我认为处理丢失数据的正确方法不是扩展窗口（通过删除 NA 值），而是通过正确的分母对每个窗口进行平均。
rn
【参考方案4】：
在 data.table 1.12.0 中添加了新的 frollmean 函数，以计算快速准确的滚动平均值，仔细处理 NA、NaN 和 +Inf、-Inf 值.
由于问题中没有可重复的示例，因此这里没有更多要解决的问题。
您可以在手册中找到有关?frollmean 的更多信息，也可以通过?frollmean 在线获取。
以下手册中的示例：
library(data.table)
d = as.data.table(list(1:6/2, 3:8/4))
# rollmean of single vector and single window
frollmean(d[, V1], 3)
# multiple columns at once
frollmean(d, 3)
# multiple windows at once
frollmean(d[, .(V1)], c(3, 4))
# multiple columns and multiple windows at once
frollmean(d, c(3, 4))
## three above are embarrassingly parallel using openmp
【讨论】：
【参考方案5】：caTools 包具有非常快速的滚动平均值/最小值/最大值/sd 和一些其他功能。我只使用过runmean 和runsd，它们是迄今为止提到的任何其他软件包中最快的。
【讨论】：
这太棒了！它是唯一能以一种简单、漂亮的方式完成此任务的函数。现在是 2018 年......【参考方案6】：您可以使用RcppRoll 获得用 C++ 编写的非常快速的移动平均线。只需调用roll_mean 函数。可以在here 找到文档。
否则，这个（较慢的）for 循环应该可以解决问题：
ma <- function(arr, n=15)
res = arr
for(i in n:length(arr))
res[i] = mean(arr[(i-n):i])
res
【讨论】：
你能否详细解释一下，这个算法是如何工作的？因为我无法理解这个想法 首先他用res = arr初始化了一个相同长度的向量。然后有一个循环，从n 或第 15 个元素开始迭代到数组的末尾。这意味着他取平均值的第一个子集是arr[1:15]，它填充了res[15]。现在，我更喜欢设置res = rep(NA, length(arr)) 而不是res = arr，所以res[1:14] 的每个元素都等于NA 而不是一个数字，我们不能完全取15 个元素的平均值。
我觉得应该是arr[(i-n+1):i]【参考方案7】：
其实RcppRoll很好。
cantdutchthis发布的代码必须在第四行更正到窗口被修复：
ma <- function(arr, n=15)
res = arr
for(i in n:length(arr))
res[i] = mean(arr[(i-n+1):i])
res
另一种处理缺失的方法是here。
第三种方法，改进cantdutchthis 代码以计算部分平均值，如下：
ma <- function(x, n=2,parcial=TRUE)
res = x #set the first values
if (parcial==TRUE)
for(i in 1:length(x))
t<-max(i-n+1,1)
res[i] = mean(x[t:i])
res
else
for(i in 1:length(x))
t<-max(i-n+1,1)
res[i] = mean(x[t:i])
res[-c(seq(1,n-1,1))] #remove the n-1 first,i.e., res[c(-3,-4,...)]
【讨论】：
【参考方案8】：这里的示例代码展示了如何使用 zoo 包中的 rollmean 函数计算居中移动平均线和追踪移动平均线。
library(tidyverse)
library(zoo)
some_data = tibble(day = 1:10)
# cma = centered moving average
# tma = trailing moving average
some_data = some_data %>%
mutate(cma = rollmean(day, k = 3, fill = NA)) %>%
mutate(tma = rollmean(day, k = 3, fill = NA, align = "right"))
some_data
#> # A tibble: 10 x 3
#> day cma tma
#> <int> <dbl> <dbl>
#> 1 1 NA NA
#> 2 2 2 NA
#> 3 3 3 2
#> 4 4 4 3
#> 5 5 5 4
#> 6 6 6 5
#> 7 7 7 6
#> 8 8 8 7
#> 9 9 9 8
#> 10 10 NA 9
【讨论】：
您可以对多个新列使用一个 mutate 调用，方法是用逗号分隔每个新列。【参考方案9】：为了补充cantdutchthis和Rodrigo Remedio的答案；
moving_fun <- function(x, w, FUN, ...)
# x: a double vector
# w: the length of the window, i.e., the section of the vector selected to apply FUN
# FUN: a function that takes a vector and return a summarize value, e.g., mean, sum, etc.
# Given a double type vector apply a FUN over a moving window from left to the right,
# when a window boundary is not a legal section, i.e. lower_bound and i (upper bound)
# are not contained in the length of the vector, return a NA_real_
if (w < 1)
stop("The length of the window 'w' must be greater than 0")
output <- x
for (i in 1:length(x))
# plus 1 because the index is inclusive with the upper_bound 'i'
lower_bound <- i - w + 1
if (lower_bound < 1)
output[i] <- NA_real_
else
output[i] <- FUN(x[lower_bound:i, ...])
output
# example
v <- seq(1:10)
# compute a MA(2)
moving_fun(v, 2, mean)
# compute moving sum of two periods
moving_fun(v, 2, sum)
【讨论】：
【参考方案10】：您可以通过以下方式计算窗口宽度为k 的向量x 的移动平均值：
apply(embed(x, k), 1, mean)
【讨论】：
data.frames 的扩展是：apply(df,rc,FUN=function(x) apply(embed(x, k),1,mean))。 rc 可以是一个或两个，分别代表行或列。【参考方案11】：
滑块包可以用于此。它的界面经过专门设计，感觉类似于 purrr。它接受任意函数，并且可以返回任何类型的输出。数据帧甚至可以逐行迭代。 pkgdown 站点是here。
library(slider)
x <- 1:3
# Mean of the current value + 1 value before it
# returned as a double vector
slide_dbl(x, ~mean(.x, na.rm = TRUE), .before = 1)
#> [1] 1.0 1.5 2.5
df <- data.frame(x = x, y = x)
# Slide row wise over data frames
slide(df, ~.x, .before = 1)
#> [[1]]
#> x y
#> 1 1 1
#>
#> [[2]]
#> x y
#> 1 1 1
#> 2 2 2
#>
#> [[3]]
#> x y
#> 1 2 2
#> 2 3 3
slider 和 data.table 的 frollapply() 的开销应该非常低（比 zoo 快得多）。 frollapply() 在这里的这个简单示例看起来要快一些，但请注意它只接受数字输入，并且输出必须是标量数值。滑块函数是完全通用的，您可以返回任何数据类型。
library(slider)
library(zoo)
library(data.table)
x <- 1:50000 + 0L
bench::mark(
slider = slide_int(x, function(x) 1L, .before = 5, .complete = TRUE),
zoo = rollapplyr(x, FUN = function(x) 1L, width = 6, fill = NA),
datatable = frollapply(x, n = 6, FUN = function(x) 1L),
iterations = 200
)
#> # A tibble: 3 x 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 slider 19.82ms 26.4ms 38.4 829.8KB 19.0
#> 2 zoo 177.92ms 211.1ms 4.71 17.9MB 24.8
#> 3 datatable 7.78ms 10.9ms 87.9 807.1KB 38.7
【讨论】：
【参考方案12】：编辑：非常高兴添加side 参数，例如移动平均线（或总和，或...） Date 向量的过去 7 天。
对于只想自己计算的人来说，无非是：
# x = vector with numeric data
# w = window length
y <- numeric(length = length(x))
for (i in seq_len(length(x)))
ind <- c((i - floor(w / 2)):(i + floor(w / 2)))
ind <- ind[ind %in% seq_len(length(x))]
y[i] <- mean(x[ind])
y
但是让它独立于mean() 会很有趣，所以你可以计算任何“移动”函数！
# our working horse:
moving_fn <- function(x, w, fun, ...)
# x = vector with numeric data
# w = window length
# fun = function to apply
# side = side to take, (c)entre, (l)eft or (r)ight
# ... = parameters passed on to 'fun'
y <- numeric(length(x))
for (i in seq_len(length(x)))
if (side %in% c("c", "centre", "center"))
ind <- c((i - floor(w / 2)):(i + floor(w / 2)))
else if (side %in% c("l", "left"))
ind <- c((i - floor(w) + 1):i)
else if (side %in% c("r", "right"))
ind <- c(i:(i + floor(w) - 1))
else
stop("'side' must be one of 'centre', 'left', 'right'", call. = FALSE)
ind <- ind[ind %in% seq_len(length(x))]
y[i] <- fun(x[ind], ...)
y
# and now any variation you can think of!
moving_average <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = mean, side = side, na.rm = na.rm)
moving_sum <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = sum, side = side, na.rm = na.rm)
moving_maximum <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = max, side = side, na.rm = na.rm)
moving_median <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = median, side = side, na.rm = na.rm)
moving_Q1 <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = quantile, side = side, na.rm = na.rm, 0.25)
moving_Q3 <- function(x, w = 5, side = "centre", na.rm = FALSE)
moving_fn(x = x, w = w, fun = quantile, side = side, na.rm = na.rm, 0.75)
【讨论】：
【参考方案13】：虽然有点慢，但您也可以使用 zoo::rollapply 对矩阵执行计算。
reqd_ma <- rollapply(x, FUN = mean, width = n)
其中x是数据集，FUN = mean是函数；您也可以将其更改为 min、max、sd 等，width 是滚动窗口。
【讨论】：
并不慢；。与基础 R 相比，它快得多。set.seed(123); x <- rnorm(1000); system.time(apply(embed(x, 5), 1, mean)); library(zoo); system.time(rollapply(x, 5, mean)) 在我的机器上它是如此之快以至于它返回一个 0 秒的时间。【参考方案14】：
可以使用runner 包来移动函数。在这种情况下，mean_run 函数。 cummean 的问题在于它不处理 NA 值，但 mean_run 处理。 runner 包还支持不规则的时间序列和窗口可以依赖日期：
library(runner)
set.seed(11)
x1 <- rnorm(15)
x2 <- sample(c(rep(NA,5), rnorm(15)), 15, replace = TRUE)
date <- Sys.Date() + cumsum(sample(1:3, 15, replace = TRUE))
mean_run(x1)
#> [1] -0.5910311 -0.2822184 -0.6936633 -0.8609108 -0.4530308 -0.5332176
#> [7] -0.2679571 -0.1563477 -0.1440561 -0.2300625 -0.2844599 -0.2897842
#> [13] -0.3858234 -0.3765192 -0.4280809
mean_run(x2, na_rm = TRUE)
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.12188853 -0.13873536
#> [7] -0.13873536 -0.14571604 -0.12596067 -0.11116961 -0.09881996 -0.08871569
#> [13] -0.05194292 -0.04699909 -0.05704202
mean_run(x2, na_rm = FALSE )
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.12188853 -0.13873536
#> [7] NA NA NA NA NA NA
#> [13] NA NA NA
mean_run(x2, na_rm = TRUE, k = 4)
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.10546063 -0.16299272
#> [7] -0.21203756 -0.39209010 -0.13274756 -0.05603811 -0.03894684 0.01103493
#> [13] 0.09609256 0.09738460 0.04740283
mean_run(x2, na_rm = TRUE, k = 4, idx = date)
#> [1] -0.187600111 -0.090220655 -0.004349696 0.168349653 -0.206571573 -0.494335093
#> [7] -0.222969541 -0.187600111 -0.087636571 0.009742884 0.009742884 0.012326968
#> [13] 0.182442234 0.125737145 0.059094786
还可以指定其他选项，例如lag，并且只滚动at 特定索引。更多内容请参见 package 和 function 文档。
【讨论】：
【参考方案15】：这是一个简单的函数，filter 演示了一种使用填充处理开始和结束 NA 的方法，并使用自定义权重计算加权平均值（filter 支持）：
wma <- function(x)
wts <- c(seq(0.5, 4, 0.5), seq(3.5, 0.5, -0.5))
nside <- (length(wts)-1)/2
# pad x with begin and end values for filter to avoid NAs
xp <- c(rep(first(x), nside), x, rep(last(x), nside))
z <- stats::filter(xp, wts/sum(wts), sides = 2) %>% as.vector
z[(nside+1):(nside+length(x))]
【讨论】：
【参考方案16】：vector_avg <- function(x)
sum_x = 0
for(i in 1:length(x))
if(!is.na(x[i]))
sum_x = sum_x + x[i]
return(sum_x/length(x))
【讨论】：
请添加描述以获取更多详细信息。 请将您的答案与问题联系起来，并包含一些显示问题已得到回答的输出。请参阅 How to Answer 以获取有关做出良好答案的指导。【参考方案17】：我将聚合与由 rep() 创建的向量一起使用。这具有使用 cbind() 一次在数据框中聚合多于 1 列的优点。以下是长度为 1000 的向量 (v) 的移动平均值 60 的示例：
v=1:1000*0.002+rnorm(1000)
mrng=rep(1:round(length(v)/60+0.5), length.out=length(v), each=60)
aggregate(v~mrng, FUN=mean, na.rm=T)
注意 rep 中的第一个参数是简单地根据向量的长度和要平均的数量为移动范围获取足够的唯一值；第二个参数保持长度等于向量长度，最后一个参数重复第一个参数的值与平均周期相同的次数。
总的来说，您可以使用多个函数（中值、最大值、最小值） - 例如所示的平均值。同样，可以使用带有 cbind 的公式在数据框中的多个（或所有）列上执行此操作。
【讨论】：
以上是关于计算移动平均线的主要内容，如果未能解决你的问题，请参考以下文章