具有基于时间的窗口的不规则时间序列上的优化滚动函数
Posted
技术标签:
【中文标题】具有基于时间的窗口的不规则时间序列上的优化滚动函数【英文标题】：optimized rolling functions on irregular time series with time-based window 【发布时间】：2013-04-04 07:56:43 【问题描述】：有没有办法使用 rollapply（来自zoo 包或类似的东西）优化函数（rollmean、rollmedian 等）来计算具有基于时间的窗口的滚动函数，而不是基于数字的滚动函数观察？我想要的很简单：对于不规则时间序列中的每个元素，我想计算一个具有 N 天窗口的滚动函数。也就是说，该窗口应包括当前观察前 N 天之前的所有观察。时间序列也可能包含重复项。
下面是一个例子。给定以下时间序列：
date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3
具有 5 天窗口的滚动中位数，向右对齐，应进行以下计算：
> c(
median(c(5)),
median(c(5,4)),
median(c(5,4,2)),
median(c(1)),
median(c(1,0)),
median(c(0,0)),
median(c(0,0,0)),
median(c(0,0,0,1)),
median(c(1,4)),
median(c(3))
)
[1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
我已经找到了一些解决方案，但它们通常很棘手，这通常意味着缓慢。我设法实现了自己的滚动函数计算。问题在于，对于很长的时间序列，中值（rollmedian）的优化版本会产生巨大的时间差异，因为它考虑了窗口之间的重叠。我想避免重新实现它。我怀疑rollapply参数有一些技巧可以使它起作用，但我无法弄清楚。提前感谢您的帮助。
【问题讨论】：
rollapply 无法做到这一点。您可以使用window 滚动您自己的函数（双关语）。
这个问答有什么帮助吗？ ***.com/questions/10465998/…
rollapply “作弊”，如果您使用 median 作为乐趣，请致电 rollmedian。比较：system.time(rollapply(runif(100000), 5, function(x) median(x))) 与 system.time(rollapply(runif(100000), 5, median))（前者慢 30 倍）。如果您想要与rollapply 在没有“作弊”的情况下达到的速度相当的速度，我可以提供一些解决方案。此外，rollmedian 还“作弊”，因为它需要奇怪的观察，所以很明显它只是定义了一个“中间”值的索引，与您尝试做的相比，这是微不足道的。
查看这个答案 (***.com/questions/20134823/…) 以了解具有基于时间的窗口的快速 Rcpp rollmean 函数。
是否可以通过填充 NA 来使时间序列有规律，然后对其应用固定大小的滚动窗口？
【参考方案1】：
从 v1.9.8 版（CRAN 2016 年 11 月 25 日）开始，data.table 已获得执行非等值连接的能力，可在此处使用。
OP 已请求
对于不规则时间序列中的每个元素，我想计算一个 具有 N 天窗口的滚动函数。也就是说，窗口应该 包括当前 N 天之前的所有观测值 观察。时间序列也可能包含重复项。
请注意，OP 已要求包括在当前观察之前最多 N 天的所有观察。这与请求当前 day 前 N 天的所有观察结果不同。
对于后者，我希望 1/11/2011 的 one 值，即 median(c(5, 4, 2)) = 4。
显然，OP 期望基于 观察 的滚动窗口限制为 N 天。因此，非等连接的连接条件也要考虑行号。
library(data.table)
n_days <- 5L
setDT(DT)[, rn := .I][
.(ur = rn, ud = date, ld = date - n_days),
on = .(rn <= ur, date <= ud, date >= ld),
median(as.double(value)), by = .EACHI]$V1
[1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
为了完整起见，基于天的滚动窗口的可能解决方案可能是：
setDT(DT)[.(ud = unique(date), ld = unique(date) - n_days), on = .(date <= ud, date >= ld),
median(as.double(value)), by = .EACHI]
date date V1 1: 2011-11-01 2011-10-27 4.0 2: 2011-11-08 2011-11-03 1.0 3: 2011-11-13 2011-11-08 0.5 4: 2011-11-14 2011-11-09 0.0 5: 2011-11-15 2011-11-10 0.0 6: 2011-11-18 2011-11-13 0.0 7: 2011-11-21 2011-11-16 2.5 8: 2011-12-05 2011-11-30 3.0
数据
library(data.table)
DT <- fread(" date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3")[
# coerce date from character string to integer date class
, date := as.IDate(date, "%d/%m/%Y")]
【讨论】：
【参考方案2】：1) rollapply 没有检查速度，但如果没有日期的出现次数超过max.dup，那么它必须是最后 5 * max.dup 条目包含最后 5 天，所以如下所示的单行函数fn 传递给rollapplyr 即可：
k <- 5
dates <- as.numeric(DF$date)
values <- DF$value
max.dup <- max(table(dates))
fn <- function(ix, d = dates[ix], v = values[ix], n = length(ix)) median(v[d >= d[n]-k])
rollapplyr(1:nrow(DF), max.dup * k, fn, partial = TRUE)
## [1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
2) sqldf 我们可以使用 SQL 自连接来执行此操作。我们将不超过 5 天的 b 行加入到每个 a 行，然后按 a 行分组，取加入的 b 行的中位数。
library(sqldf)
k <- 5
res <- fn$sqldf("select a.date, a.value, median(b.value) median
from DF a
left join DF b on b.date between a.date - $k and a.date and b.rowid <= a.rowid
group by a.rowid")
给予：
res$median
## [1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
注意：我们将此用于DF：
Lines <- "
date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3
"
DF <- read.table(text = Lines, header = TRUE)
DF$date <- as.Date(DF$date, format = "%d/%m/%Y")
【讨论】：
【参考方案3】：我推荐使用runner 包，该包经过优化以执行本主题中要求的操作。转至documentation 中的 Windows 取决于日期 部分，以获得进一步说明。
为了解决您的任务，可以使用runner 函数，它可以在运行的窗口中执行任何 R 函数。单线在这里：
df <- read.table(
text = "date value
2011-11-01 5
2011-11-01 4
2011-11-01 2
2011-11-08 1
2011-11-13 0
2011-11-14 0
2011-11-15 0
2011-11-18 1
2011-11-21 4
2011-12-05 3", header = TRUE, colClasses = c("Date", "integer"))
library(runner)
runner(df$value, k = 5, idx = df$date, f = median)
[1] 5.0 4.5 4.0 1.0 0.0 0.0 0.0 0.0 2.5 3.0
附：应该知道 5 天窗口是 [i-4, i-3, i-2, i-1, i] 而不是 (i-5):i（6 天窗口）。下面的插图可以更好地解释这个概念。
我在 5 天窗口上做了示例，但如果想按照 OP 的要求重现结果，可以指定 6 天窗口：
identical(
runner(df$value, k = 6, idx = df$date, f = median),
c(5.0, 4.5, 4.0, 1.0, 0.5, 0.0, 0.0, 0.0, 2.5, 3.0)
)
# [1] TRUE
【讨论】：
【参考方案4】：大多数答案建议插入 NA 以使时间序列有规律。 但是，在长时间序列的情况下，这可能会很慢。此外，它不适用于不能与 NA 一起使用的功能。
rollapply（动物园包）的width参数可以是一个列表（详见rollapply的帮助）。基于此，我编写了一个函数，该函数创建一个列表，用于 rollapply 作为宽度参数。如果移动窗口是时间而不是基于索引的，则该函数提取不规则动物园对象的索引。所以动物园对象的索引应该是实际时间。
# Create a zoo object where index represents time (e.g. in seconds)
d <- zoo(c(1,1,1,1,1,2,2,2,2,2,16,25,27,27,27,27,27,31),
c(1:5,11:15,16,25:30,31))
# Create function
createRollapplyWidth = function(zoodata, steps, window )
mintime = min(time(zoodata))
maxtime = max(time(zoodata))
spotstime = seq(from = mintime , to = maxtime, by = steps)
spotsindex = list()
for (i in 1:length(spotstime))
spotsindex[[i]] = as.numeric(which(spotstime[i] <= time(zoodata) & time(zoodata) < spotstime[i] + window))
rollapplywidth = list()
for (i in 1:length(spotsindex))
if (!is.na(median(spotsindex[[i]])) )
rollapplywidth[[round(median(spotsindex[[i]]))]] = spotsindex[[i]] - round(median(spotsindex[[i]]))
return(rollapplywidth)
# Create width parameter for rollapply using function
rollwidth = createRollapplyWidth(zoodata = d, steps = 5, window = 5)
# Use parameter in rollapply
result = rollapply(d, width = rollwidth , FUN = sum, na.rm = T)
result
限制：不是基于日期，而是基于时间，以秒为单位。 rollapply 的参数“partial”不起作用。
【讨论】：
【参考方案5】：这是我对问题的修补。如果那种得到你想要的（我不知道它在速度方面是否令人满意），我可以把它写成更详细的答案（即使它基于@rbatt的想法）。
library(zoo)
library(dplyr)
# create a long time series
start <- as.Date("1800-01-01")
end <- as.Date(Sys.Date())
df <- data.frame(V1 = seq.Date(start, end, by = "day"))
df$V2 <- sample(1:10, nrow(df), replace = T)
# make it an irregular time series by sampling 10000 rows
# including allowing for duplicates (replace = T)
df2 <- df %>%
sample_n(10000, replace = T)
# create 'complete' time series & join the data & compute the rolling median
df_rollmed <- data.frame(V1 = seq.Date(min(df$V1), max(df$V1), by = "day")) %>%
left_join(., df2) %>%
mutate(rollmed = rollapply(V2, 5, median, na.rm = T, align = "right", partial = T)) %>%
filter(!is.na(V2)) # throw out the NAs from the complete dataset
【讨论】：
以上是关于具有基于时间的窗口的不规则时间序列上的优化滚动函数的主要内容，如果未能解决你的问题，请参考以下文章