如何通过顺序操作加速 R 循环

Question

我有一个模型，它有多个条件并返回一个值，该值取决于下一个预测。假设给定 A 和 B 的时间序列，模型返回 C 变量的值，该值又用于估计 D 的值。在沿着新 A 和 B 的下一次迭代中，模型还使用估计的 D 作为输入：

df = data.frame(A = sample(-5:5, 10000, replace = TRUE),
                 B = sample(-5:5, 10000, replace = TRUE),
                 C = 0,
                D=0)

for(i in 1:nrow(df)){
  
    if (df$A[i]< 0 & df$B[i]>0){     
      df$C[i]<-df$B[i]
    
      } else if(df$A[i]==0  & df$B[i]==0 ){ 
      df$C[i]<-0
    
      }  else {
      df$C[i]<-df$A[i]+df$B[i]-df$D[i]  
        }
    
    df$D[i+1]<-ifelse(df$D[i]<=-df$C[i],0,df$D[i]+df$C[i]) # this is a cumulative sum-reset function
    
}

尽管代码运行良好，但由于我有数十万个观察结果，所以速度非常慢。我将不胜感激任何可以加快速度的建议。

Answer 1

由于每一行都依赖于前一行的结果，因此很难以可以利用 R 向量化的方式编写。在这种情况下，我们在 Rcpp 中编写代码会获得巨大的优势。

library(Rcpp)

cppFunction('

DataFrame f_Rcpp(DataFrame df) {

  NumericVector A = df["A"];
  NumericVector B = df["B"];
  NumericVector C = df["C"];
  NumericVector D = df["D"];

  for(int i = 0; i < (df.nrows() - 1); ++i) {
    
    if (A[i] < 0 && B[i] > 0) {     
      C[i] = B[i];
      
    } else if(A[i] == 0 && B[i] == 0 ) { 
      C[i] = 0;
      
    }  else {
      C[i] = A[i] + B[i] - D[i];
    }
    
    if(D[i] <= -C[i]) {
    D[i+1] = 0;
    } else {
    D[i+1] = D[i] + C[i]; 
    }
  }
  return(df);
}
            
')

如果我们将您自己的代码包装为一个函数以便我们可以对其进行比较，我们会看到我们的 Rcpp 函数给出了相同的结果：

f_R <- function(df) {
  for(i in 1:(nrow(df) - 1)) {
    
    if (df$A[i] < 0 & df$B[i] > 0) {     
      df$C[i] <- df$B[i]
      
    } else if(df$A[i] == 0 & df$B[i] == 0 ){ 
      df$C[i] <- 0
      
    }  else {
      df$C[i] <- df$A[i] + df$B[i] - df$D[i]  
    }
    
    df$D[i+1] <- ifelse(df$D[i] <= -df$C[i], 0, df$D[i] + df$C[i]) 
    
  }
  return(df)
}

res1 <- f_R(df)
res2 <- f_Rcpp(df)

identical(res1, res2)
#> [1] TRUE

但是看看当我们进行基准测试时会发生什么：

microbenchmark::microbenchmark(f_R(df), f_Rcpp(df), times = 10)
#> Unit: microseconds
#>       expr         min        lq         mean      median          uq         max neval cld
#>    f_R(df) 1746032.401 1793779.0 1794274.9209 1802222.051 1810686.801 1815285.001    10   b
#> f_Rcpp(df)     567.701     585.9     610.1607     601.851     642.801     650.101    10  a 

Rcpp 函数在不到一毫秒的时间内处理所有 10,000 行，而在基本 R 中则需要近 2 秒。Rcpp 版本的速度几乎快了 3,000 倍。

编辑

要使用您自己的数据进行此操作，请尝试：

cppFunction('

DataFrame f_Rcpp(DataFrame df, NumericVector v) {
  NumericVector A = df["Tav"];
  NumericVector B = df["dprcp"];
  NumericVector C = df["dSWE"];
  NumericVector D = df["simSWE"];
  NumericVector E = df["dSWElag"];

  for(int i = 5; i < (df.nrows() - 1); ++i) {
    if (A[i] < -1 && B[i] > 0) {     
      C[i] = B[i];
    } else if(A[i] < -1 && B[i] == 0 ) { 
      C[i] = 0;
    }  else {
      C[i] = v[i];
    }
    
    if(D[i-1] <= -C[i]) {
      D[i] = 0;
    } else {
      D[i] = D[i-1] + C[i]; 
    }
    E[i + 1] = C[i];
  }

  df["dSWE"] = C;
  df["simSWE"] = D;
  df["dSWElag"] = E;
  
  return(df);
}        
')

你可以这样称呼：

preds <- predict(svm_model,station)

station2 <- f_Rcpp(station, preds)