r - 当矢量化不可行时，在数据框中迭代行的 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 循环，无法很好地扩展到更大的数据集:

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)))
}

x
#> # A tibble: 10 x 4
#>        t     a        b        c
#>    <int> <dbl>    <dbl>    <dbl>
#>  1     1  100  7.50e- 1 9.00e- 1
#>  2     2  110  8.98e- 1 6.57e- 1
#>  3     3  120  7.31e- 1 1.37e- 4
#>  4     4  130  2.08e- 1 6.96e- 1
#>  5     5  140  6.70e- 1 8.82e- 1
#>  6     6  154  1.03e+ 2 1.34e+ 4
#>  7     7  169. 1.75e+ 4 2.45e+ 6
#>  8     8  186. 3.26e+ 6 5.02e+ 8
#>  9     9  205. 6.68e+ 8 1.13e+11
#> 10    10  225. 1.51e+11 2.80e+13

Answer 1

我认为对于这种本质上的迭代过程来说，真的很难击败 for环形。 @Shree 提出的方法取决于 NA 是连续的并从已知位置开始。
这是我对您的循环的轻微改进，我认为它更具可读性，速度提高了大约 2.5 倍，并且可能比将矢量化操作与循环相结合的方法更好地扩展。通过完全离开 tidyverse 并采用真正一次对每一行起作用的 rowwise 循环，我们在两个方面都获得了一些效率:

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 解决方案。