26

是否有某种方法可以使用 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参数有一些技巧可以使它起作用，但我无法弄清楚。在此先感谢您的帮助。

4

5 回答 5

5

从 v1.9.8 版本开始（CRAN 2016 年 11 月 25 日），已获得执行可在此处使用的非 equi 连接的能力。

OP已要求

对于不规则时间序列中的每个元素，我想计算一个具有 N 天窗口的滚动函数。也就是说，该窗口应包括当前观察前 N 天的所有观察。时间序列也可能包含重复项。

请注意，OP 已要求在当前观察前 N 天包括所有观察结果。这与请求当天前 N 天的所有观察结果不同。

对于后者，我希望有一个值1/11/2011，即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")]
于 2018-12-26T13:54:43.637 回答
3

1）rollapply没有检查速度，但如果没有日期超过出现次数，那么它必须是最后 5 * max.dup 条目包含最后 5 天，因此传递给下面显示max.dup的单行函数将执行此操作：fnrollapplyr

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 天前将a这些行加入到每一行，然后按行分组，取加入行的中位数。bab

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")
于 2015-11-24T13:10:36.667 回答
2

我建议使用经过优化的runner包来执行本主题中要求的操作。根据文档中的日期转到 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

PS 一应该知道，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
于 2019-10-27T11:34:22.307 回答
1

大多数答案都建议插入 NA 以使时间序列有规律。但是，在长时间序列的情况下，这可能会很慢。此外，它不适用于不能与 NA 一起使用的功能。

rollapply (zoo package) 的 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 的参数“部分”不起作用。

于 2017-06-07T08:23:07.607 回答
0

这是我对这个问题的修补。如果那种得到你想要的（我不知道它在速度方面是否令人满意），我可以把它写成更详细的答案（即使它基于@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
于 2015-11-24T11:08:10.953 回答