当矢量化不可行时，在数据框中迭代行的 tidyverse 方法是什么？

Question

我想知道当第 n 行变量的值取决于第 n-1 行和/或第 n-2 行变量的值时，迭代数据帧行的最佳方法。理想情况下，我想以“tidyverse”的方式来做到这一点，也许使用 purrr::pmap()。

例如，假设我有这个数据框：

library(dplyr)

x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))

x
#> # A tibble: 10 x 4
#>        t     a      b         c
#>    <int> <dbl>  <dbl>     <dbl>
#>  1     1   100  0.750  0.900   
#>  2     2   110  0.898  0.657   
#>  3     3   120  0.731  0.000137
#>  4     4   130  0.208  0.696   
#>  5     5   140  0.670  0.882   
#>  6     6    NA NA     NA       
#>  7     7    NA NA     NA       
#>  8     8    NA NA     NA       
#>  9     9    NA NA     NA       
#> 10    10    NA NA     NA

我已经知道时间 (t) = 5 之前的值。除此之外，我希望使用以下公式预测值：

a = lag(a) * 1.1
b = a * lag(b)
c = b * lag(a, 2)

这段代码实现了所需的输出，但它是一个笨拙、可怕的 for 循环，无法很好地扩展到更大的数据集：

a = lag(a) * 1.1
b = a * lag(b)
c = b * lag(a, 2)

Answer 1

我认为对于这种本质上迭代的过程来说，确实很难击败循环for。@Shree 提出的方法取决于 NA 是连续的并且从已知点开始。

这是我对循环的轻微改进，我认为它更具可读性，速度提高了大约 2.5 倍，并且可能比将矢量化操作与循环相结合的方法更好地扩展。通过完全摆脱 tidyverse 并采用真正一次只处理每一行的行循环，我们在这两方面都获得了一些效率：

method_peter <- function(x){
  for(i in 2:nrow(x)){
    x[i, "a"] <- ifelse(is.na(x[i, "a"]), x[i - 1, "a"] * 1.1,       x[i, "a"])
    x[i, "b"] <- ifelse(is.na(x[i, "b"]), x[i, "a"] * x[i - 1, "b"], x[i, "b"])
    x[i, "c"] <- ifelse(is.na(x[i, "c"]), x[i, "b"] * x[i - 2, "a"], x[i, "c"])
  }
  return(x)
}

毫无疑问，可能会提高效率，当然这是用 C++ 重写它的理想选择:)。

这大约是您的方法的两倍，如下所示：

method_matt <- function(x){
  for(i in 1:nrow(x)) {
    x <- x %>%
      mutate(a = if_else(!is.na(a), a, lag(a, 1) * 1.1),
             b = if_else(!is.na(b), b, a * lag(b, 1)),
             c = if_else(!is.na(c), c, b * lag(a, 2)))
  }
  return(x)
}

set.seed(123)
x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))

stopifnot(identical(method_matt(x), method_peter(x)))

library(microbenchmark)
microbenchmark(
  method_matt(x),
  method_peter(x)
)

返回：

Unit: milliseconds
            expr     min       lq     mean   median      uq     max neval
  method_matt(x) 24.1975 25.50925 30.64438 26.33310 31.8681 74.5093   100
 method_peter(x) 10.0005 10.56050 13.33751 11.06495 13.5913 42.0568   100

@Shree 的方法又快得多，并且非常适合示例数据，但我不确定它是否足够灵活以适用于您的所有用例。

如果有的话，我希望看到一个 tidyverse 解决方案。

Answer 2

我认为没有任何简单的方法可以tidyverse进行行依赖性的计算。带有Reduce或 的东西gather + spread是可能的，但我不希望它们在可读性上得分。

无论如何，从好的方面来说，您的计算可以使用dplyr和zoo包进行矢量化 -

x %>% 
  mutate(
    a = ifelse(is.na(a), na.locf(a) * 1.1^(t-5), a),
    b = ifelse(is.na(b), na.locf(b) * c(rep(1, 5), cumprod(a[6:n()])), b),
    c = ifelse(is.na(c), b * lag(a, 2), c)
  )

 # A tibble: 10 x 4
 t     a        b        c
 <int> <dbl>    <dbl>    <dbl>
 1     1  100  1.85e- 1 9.43e- 1
 2     2  110  7.02e- 1 1.29e- 1
 3     3  120  5.73e- 1 8.33e- 1
 4     4  130  1.68e- 1 4.68e- 1
 5     5  140  9.44e- 1 5.50e- 1
 6     6  154  1.45e+ 2 1.89e+ 4
 7     7  169. 2.46e+ 4 3.45e+ 6
 8     8  186. 4.59e+ 6 7.07e+ 8
 9     9  205. 9.40e+ 8 1.59e+11
10    10  225. 2.12e+11 3.95e+13

数据 -

set.seed(2)
x <- tibble(t = c(1:10),
            a = c(seq(100, 140, 10), rep(NA_real_, 5)),
            b = c(runif(5), rep(NA_real_, 5)),
            c = c(runif(5), rep(NA_real_, 5)))