logo标签列表

可调窗口和步长的滚动窗口函数，适用于不规则间隔观测数据。

rtime-seriessliding-window
9
假设有一个2列数据框，其中一列是时间或距离，按顺序递增，另一列是观测值，可能会有一些NA。如何有效地使用滑动窗口函数，在持续时间为X（例如5秒）的窗口中获取某些统计信息，例如平均值，在Y秒（例如2.5秒）内滑动窗口，重复此过程...窗口中的观测数量基于时间列，因此每个窗口中的观测数量和滑动窗口的观测数量都可能不同该函数应接受任何窗口大小，最多到观测数，并具有步长。

以下是样本数据（请参见“ 编辑：”以获取更大的样本集）

set.seed(42)
dat <- data.frame(time = seq(1:20)+runif(20,0,1))
dat <- data.frame(dat, measure=c(diff(dat$time),NA_real_))
dat$measure[sample(1:19,2)] <- NA_real_
head(dat)
      time   measure
1 1.914806 1.0222694
2 2.937075 0.3490641
3 3.286140        NA
4 4.830448 0.8112979
5 5.641746 0.8773504
6 6.519096 1.2174924

期望输出针对5秒窗口、2.5秒步长、第一个窗口从-2.5到2.5的特定情况，na.rm=FALSE：

 [1] 1.0222694
 [2]        NA
 [3]        NA
 [4] 1.0126639
 [5] 0.9965048
 [6] 0.9514456
 [7] 1.0518228
 [8]        NA
 [9]        NA
[10]        NA

解释：在期望的输出中，第一个窗口查找-2.5到2.5之间的时间。该窗口内有一次测量观察值，并且它不是NA，因此我们得到观察值1.0222694。下一个窗口是从0到5，窗口内有一个NA，所以我们得到NA。对于从2.5到7.5的窗口也是如此。下一个窗口是从5到10。窗口内有5个观察值，没有一个是NA。因此，我们得到这5个观察值的平均值（即mean(dat[dat$time>5&dat$time<10,'measure'])）。
我尝试过的方法：以下是我针对步长为窗口持续时间的1/2的特定情况所尝试的方法：
windo <- 5  # duration in seconds of window

# partition into groups depending on which window(s) an observation falls in
# When step size >= window/2 and < window, need two grouping vectors
leaf1 <- round(ceiling(dat$time/(windo/2))+0.5)
leaf2 <- round(ceiling(dat$time/(windo/2))-0.5) 

l1 <- tapply(dat$measure, leaf1, mean)
l2 <- tapply(dat$measure, leaf2, mean)

as.vector(rbind(l2,l1))

不灵活，不优雅，不高效。如果步长不是1/2窗口大小，那么这种方法就行不通。

对于这种问题，您有什么通用解决方案的想法吗？任何解决方案都可以接受。速度越快越好，但我更喜欢使用基本R、数据表、Rcpp和/或并行计算的解决方案。在我的真实数据集中，有几百万个观察结果存储在数据框列表中（最大数据框约为400,000个观察结果）。

以下是附加信息： 更大的样本集

编辑：根据要求，这里有一个更大、更现实的示例数据集，其中包含许多NA和最小时间跨度（~0.03）。确切地说，数据框列表包含像上面那个小数据框一样的小数据框，以及像下面这样的更大的数据框：

set.seed(42)
dat <- data.frame(time = seq(1:50000)+runif(50000, 0.025, 1))
dat <- data.frame(dat, measure=c(diff(dat$time),NA_real_))
dat$measure[sample(1:50000,1000)] <- NA_real_
dat$measure[c(350:450,3000:3300, 20000:28100)] <- NA_real_
dat <- dat[-c(1000:2000, 30000:35000),] 

# a list with a realistic number of observations:
dat <- lapply(1:300,function(x) dat)
1
你看过RcppRoll及其相关内容吗？我在这个问题中编写了一个很酷的窗口平均函数，是否与你所需的类似？ - bright-star
@TrevorAlexander 感谢您指引我使用 RcppRoll；我会去看一下。至于您编写的函数，据我所知，窗口是基于观测数量而不是时间持续期间，这不是我想要的。 - Jota
是的，我认为你需要像你在问题中所提到的那样编写代码，将时间间隔分成离散的索引。 - bright-star
1
我们需要一个更大的真实样本集：其中包含适量的NA值，并且时间维度上的最小间距得到了充分的体现。 - IRTFM
5个回答
8

以下是使用Rcpp的尝试。该函数假定数据已按时间排序。建议进行更多测试并进行调整。

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
NumericVector rollAverage(const NumericVector & times, 
                          NumericVector & vals, 
                          double start,
                          const double winlen, 
                          const double winshift) {
  int n = ceil((max(times) - start) / winshift);
  NumericVector winvals;
  NumericVector means(n);
  int ind1(0), ind2(0);
  for(int i=0; i < n; i++) {
    if (times[0] < (start+winlen)) {
      while((times[ind1] <= start) & 
                (times[ind1+1] <= (start+winlen)) & 
                (ind1 < (times.size() - 1))) {
        ind1++;
      }    

      while((times[ind2+1] <= (start+winlen)) & (ind2 < (times.size() - 1))) {
        ind2++;
      }  

      if (times[ind1] >= start) {
        winvals = vals[seq(ind1, ind2)];
        means[i] = mean(winvals);
      } else {
        means[i] = NA_REAL;
      }
      } else {
        means[i] = NA_REAL;
    }

    start += winshift;    
  }

   return means;
}

测试它:

set.seed(42)
dat <- data.frame(time = seq(1:20)+runif(20,0,1))
dat <- data.frame(dat, measure=c(diff(dat$time),NA_real_))
dat$measure[sample(1:19,2)] <- NA_real_

rollAverage(dat$time, dat$measure, -2.5, 5.0, 2.5)
#[1] 1.0222694        NA        NA 1.0126639 0.9965048 0.9514456 1.0518228        NA        NA        NA

使用data.table处理数据框的列表：

set.seed(42)
dat <- data.frame(time = seq(1:50000)+runif(50000, 0.025, 1))
dat <- data.frame(dat, measure=c(diff(dat$time),NA_real_))
dat$measure[sample(1:50000,1000)] <- NA_real_
dat$measure[c(350:450,3000:3300, 20000:28100)] <- NA_real_
dat <- dat[-c(1000:2000, 30000:35000),] 

# a list with a realistic number of observations:
dat <- lapply(1:300,function(x) dat)

library(data.table)
dat <- lapply(dat, setDT)
for (ind in seq_along(dat)) dat[[ind]][, i := ind]
#possibly there is a way to avoid these copies?

dat <- rbindlist(dat)

system.time(res <- dat[, rollAverage(time, measure, -2.5, 5.0, 2.5), by=i])
#user  system elapsed 
#1.51    0.02    1.54 
print(res)
#           i        V1
#      1:   1 1.0217126
#      2:   1 0.9334415
#      3:   1 0.9609050
#      4:   1 1.0123473
#      5:   1 0.9965922
#     ---              
#6000596: 300 1.1121296
#6000597: 300 0.9984581
#6000598: 300 1.0093060
#6000599: 300        NA
#6000600: 300        NA
是的，抱歉。我忘记删除这行代码了。我已经进行了编辑，但现在无法测试（稍后会尝试）。希望它仍然能正常工作。 - Roland
我现在在Windows计算机上运行它，编译器抱怨'vals'是一个常量。所以我也改变了它。由于函数的更改和不同的CPU速度，计时是不同的。 - Roland
它运作得很好！速度快，易于使用。缺点是您需要硬编码您想要使用的函数（例如，在本例中为“mean”）。据我所知，在第一次之前完全出现窗口时存在问题（即参见“testdf <- data.frame(time=10:40, measure=30:0) rollAverage2(testdf$time, testdf$measure, 0, 5, 1)”）。 - Jota
可能有一种方法可以将一个R函数传递给它，当然它需要一些输入检查，正如你注意到的一些边缘情况需要修复（我已经修复了你发现的那个）。剩下的就交给你了。 - Roland
通常来说，一个功能越专业化，它就越有效率。如果您将一个R函数传递给此函数，则会因性能下降而付出代价。 - Roland
2

这里有一个函数，可以为您的小型数据框提供相同的结果。它并不特别快：在第二个 dat 示例中的较大数据集上运行需要几秒钟。

rolling_summary <- function(DF, time_col, fun, window_size, step_size, min_window=min(DF[, time_col])) {
    # time_col is name of time column
    # fun is function to apply to the subsetted data frames
    # min_window is the start time of the earliest window

    times <- DF[, time_col]

    # window_starts is a vector of the windows' minimum times
    window_starts <- seq(from=min_window, to=max(times), by=step_size)

    # The i-th element of window_rows is a vector that tells us the row numbers of
    # the data-frame rows that are present in window i 
    window_rows <- lapply(window_starts, function(x) { which(times>=x & times<x+window_size) })

    window_summaries <- sapply(window_rows, function(w_r) fun(DF[w_r, ]))
    data.frame(start_time=window_starts, end_time=window_starts+window_size, summary=window_summaries)
}

rolling_summary(DF=dat,
                time_col="time",
                fun=function(DF) mean(DF$measure),
                window_size=5,
                step_size=2.5,
                min_window=-2.5)
非常好。从“Rprof”输出的解释中，我认为lapply（window_starts，function（x）which（times> = x＆times <x + window_size））是最慢的一行，但我还没想出如何改进它。我正在尝试使用data.table来提高性能，但到目前为止，我只让事情变得更慢了。 - Jota
2

下面是一些函数，它们会在您的第一个示例中产生相同的输出：

partition <- function(x, window, step = 0){
    a = x[x < step]    
    b = x[x >= step]
    ia = rep(0, length(a))
    ib = cut(b, seq(step, max(b) + window, by = window))    
    c(ia, ib)
}

roll <- function(df, window, step = 0, fun, ...){
    tapply(df$measure, partition(df$time, window, step), fun, ...)
}

roll_steps <- function(df, window, steps, fun, ...){
    X = lapply(steps, roll, df = df, window = window, fun = fun, ...)
    names(X) = steps
    X
}

您的第一个示例的输出：

> roll_steps(dat, 5, c(0, 2.5), mean)
$`0`
        1         2         3         4         5 
       NA 1.0126639 0.9514456        NA        NA 

$`2.5`
        0         1         2         3         4 
1.0222694        NA 0.9965048 1.0518228        NA

您可以轻松地忽略缺失的值，方法如下：
> roll_steps(dat, 5, c(0, 2.5), mean, na.rm = TRUE)
$`0`
        1         2         3         4         5 
0.7275438 1.0126639 0.9514456 0.9351326       NaN 

$`2.5`
        0         1         2         3         4 
1.0222694 0.8138012 0.9965048 1.0518228 0.6122983

这也可以用于数据框的列表：

> x = lapply(dat2, roll_steps, 5, c(0, 2.5), mean)
2
好的，那这样怎么样？
library(data.table)
dat <- data.table(dat)
setkey(dat, time)

# function to compute a given stat over a time window on a given data.table
window_summary <- function(start_tm, window_len, stat_fn, my_dt) {
  pos_vec <- my_dt[, which(time>=start_tm & time<=start_tm+window_len)]
  return(stat_fn(my_dt$measure[pos_vec]))
}

# a vector of window start times
start_vec <- seq(from=-2.5, to=dat$time[nrow(dat)], by=2.5)

# sapply'ing the function above over vector of start times 
# (in this case, getting mean over 5 second windows)
result <- sapply(start_vec, window_summary, 
                 window_len=5, stat_fn=mean, my_dt=dat)

在我的电脑上，它处理了你的大型数据集的前20000行，用时13.06781秒；全部行用时51.58614秒。

我猜这比詹姆斯的解决方案慢，但或许看到另一种方法也会有所帮助。 - arvi1000
2
这里有另一种使用纯粹的data.table方法及其between函数的尝试。
与上述答案（除@Rolands的答案）进行了比较，并且似乎是最优化的。尚未测试是否存在错误，但如果您喜欢，我可以扩展答案。
使用上面提到的dat数据。
library(data.table)
Rollfunc <- function(dat, time, measure, wind = 5, slide = 2.5, FUN = mean, ...){
  temp <- seq.int(-slide, max(dat$time), by = slide)
  temp <- cbind(temp, temp + wind)
  setDT(dat)[, apply(temp, 1, function(x) FUN(measure[between(time, x[1], x[2])], ...))]
}

Rollfunc(dat, time, measure, 5, 2.5)

## [1] 1.0222694        NA        NA 1.0126639 0.9965048 0.9514456 1.0518228        NA        NA
## [10]        NA

你还可以指定函数及其参数，例如：
Rollfunc(dat, time, measure, 5, 2.5, max, na.rm = TRUE)

也可以工作

编辑: 我对@Roland进行了一些基准测试，他的方法明显胜出（远远超过其他方法），因此我会选择Rcpp方法

它的胜率有多高？我很好奇，因为data.table通常具有非常强大的性能。如果除了“使其变成c”之外还可以获得相当大的性能提升，那么我认为Hadley Wickam（等人）非常希望将其推广并使R在这方面胜出。 - EngrStudent
1
@EngrStudent 请忽略这个答案，因为这是我还不太熟悉data.table时的旧答案。如果你在任何与data.table相关的地方看到apply(..., 1, ...)，你可以放心地给我点个踩。我认为今天我会通过类似这样的方式来解决这个问题，但是三年过去了，我懒得修改这个答案。 - David Arenburg
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接