R：具有可调节窗口和步长的滚动窗口功能，用于不规则间隔的观察

Posted 2023-02-14

【中文标题】R：具有可调节窗口和步长的滚动窗口功能，用于不规则间隔的观察

【英文标题】：R: Rolling window function with adjustable window and step-size for irregularly spaced observations

【发布时间】：2014-09-04 21:25:58

【问题描述】：

假设有一个 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

Desired Output 针对 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

我尝试了什么：以下是我针对步长为窗口持续时间 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、data.table、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)

【问题讨论】：

你见过 RcppRoll 和它的朋友吗？我在this question 中做了一个很酷的窗口平均函数；这和你追求的相似吗？

@TrevorAlexander 感谢您将我指向RcppRoll；我会看看的。至于您编写的函数，据我所知，该窗口基于观察次数而不是时间持续时间，这不是我所追求的。

是的，我认为您需要像问题中那样的代码来将时间持续时间分类为离散索引。

我们需要一个更大的真实样本集：一个具有真实数量的 NA 并且表示沿时间维度的最小间距的样本集。

【参考方案1】：

这是对 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.frames 列表（使用 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

【讨论】：

是的，对不起。我忘了删除这条线。我已编辑，但现在无法测试（将在今天晚些时候尝试）。希望它仍然有效。

我现在在一台 win 机器上运行它，编译器抱怨 vals 是一个常量。所以，我也改变了这一点。由于功能的变化和 CPU 速度的不同，时序有所不同。

效果很好！它快速且易于使用。缺点是你需要硬编码你想使用的函数（例如mean在这种情况下。），afaik。当窗口在第一次之前完全出现时会出现问题（即参见testdf <- data.frame(time=10:40, measure=30:0) rollAverage2(testdf$time, testdf$measure, 0, 5, 1)。

可能有一种方法可以将 R 函数传递给它，当然它需要一些输入检查，并且您注意到需要修复一些边缘情况（我已经修复了您找到的那个）。剩下的就交给你了。

一般来说，功能越专业，效率越高。如果您将 R 函数传递给它，您将为此付出性能损失。

【参考方案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)

【讨论】：

+1 非常好。在我看来（根据我对Rprof 输出的解释）lapply(window_starts, function(x) which(times>=x & times<x+window_size)) 是最慢的线路，但我还没有弄清楚如何改进它。我正在尝试使用data.table 来提高性能，但到目前为止我只是让事情变慢了。

【参考方案3】：

以下是一些函数，它们将在您的第一个示例中提供相同的输出：

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 

这也可以用于 data.frames 列表：

> x = lapply(dat2, roll_steps, 5, c(0, 2.5), mean)

【参考方案4】：

好的，这个怎么样。

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)

在我的机器上，它在 13.06781 秒内处理大数据集的前 20,000 行； 51.58614 秒内的所有行

【讨论】：

（我想这比 James 的解决方案要慢，但无论如何可能有助于查看另一种方法）

【参考方案5】：

这是使用纯data.table 方法及其between 函数的另一种尝试。

已将Rprof 与上述答案（@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 请忽略这个答案，因为当我不太了解 data.table 时，这是一个非常古老的答案。如果您在data.table 附近的任何地方看到apply(..., 1, ...)，我允许您投反对票。我想今天我会通过做类似this 的事情来解决这个问题，但我懒得在 3 年后修改这个答案。