18

逐次累積計算

各行で計算される値は、前の行で計算された結果に依存する時系列計算を行う必要があります。の便利さを利用したいと考えていますdata.table。実際の問題は水文モデルです。つまり、各時間ステップで降雨量を加算し、現在の水量の関数として流出と蒸発を差し引く、累積的な水収支計算です。データセットには、さまざまな盆地とシナリオ (グループ) が含まれています。ここでは、問題の簡単な図を使用します。

計算の単純化された例は、時間ステップ (行) ごとに次のようになりますi。

 v[i] <- a[i] + b[i] * v[i-1]

aとbはパラメータ値のベクトルで、vは結果のベクトルです。最初の行 ( i == 1) では、 の初期値はvと見なされv0 = 0ます。

最初の試み

私の最初の考えは、で使用することshift()でしたdata.table。望ましい結果を含む最小限の例は次のとおりv.ansです。

library(data.table)        # version 1.9.7
DT <- data.table(a = 1:4, 
                 b = 0.1,
                 v.ans = c(1, 2.1, 3.21, 4.321) )
DT
#    a   b v.ans
# 1: 1 0.1 1.000
# 2: 2 0.1 2.100
# 3: 3 0.1 3.210
# 4: 4 0.1 4.321

DT[, v := NA]   # initialize v
DT[, v := a + b * ifelse(is.na(shift(v)), 0, shift(v))][]
#    a   b v.ans v
# 1: 1 0.1 1.000 1
# 2: 2 0.1 2.100 2
# 3: 3 0.1 3.210 3
# 4: 4 0.1 4.321 4

これは機能しません。shift(v)元の列のコピーがv1 行分シフトされるためです。への割り当ての影響を受けませんv。

cumsum() と cumprod() を使用して方程式を構築することも検討しましたが、それもうまくいきません。

ブルートフォースアプローチ

したがって、便宜上、関数内で for ループを使用します。

vcalc <- function(a, b, v0 = 0) {
  v <- rep(NA, length(a))      # initialize v
  for (i in 1:length(a)) {
    v[i] <- a[i] + b[i] * ifelse(i==1, v0, v[i-1])
  }
  return(v)
}

この累積関数は、data.table で正常に機能します。

DT[, v := vcalc(a, b, 0)][]
#    a   b v.ans     v
# 1: 1 0.1 1.000 1.000
# 2: 2 0.1 2.100 2.100
# 3: 3 0.1 3.210 3.210
# 4: 4 0.1 4.321 4.321
identical(DT$v, DT$v.ans)
# [1] TRUE

私の質問

data.table私の質問は、 for ループや関数定義を使用せずに、この計算をより簡潔で効率的な方法で記述できないかということです。set()おそらく使用していますか？

それとも、より良いアプローチがありますか？

編集：より良いループ

以下の David の Rcpp ソリューションifelse()は、forループから を削除するきっかけになりました。

vcalc2 <- function(a, b, v0 = 0) {
  v <- rep(NA, length(a))
  for (i in 1:length(a)) {
    v0 <- v[i] <- a[i] + b[i] * v0
  }
  return(v)
}

vcalc2()よりも 60% 高速ですvcalc()。

4

2 に答える 2

7

「data.table-way」を使用せず、まだ for ループを使用しているため、100% 探しているものではない可能性があります。ただし、このアプローチはより高速である必要があります (data.table と data.table-way を使用してコードを高速化すると仮定します)。Rcpp を利用してHydroFun、他の関数と同様に R で使用できる という短い関数を作成します (最初に関数をソースするだけです)。私の直感では、data.table の方法 (存在する場合) は、閉形式の解を計算できないため、かなり複雑であることがわかります (ただし、この点では間違っている可能性があります...)。

私のアプローチは次のようになります。

Rcpp 関数は次のようになります (ファイル内: hydrofun.cpp):

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector HydroFun(NumericVector a, NumericVector b, double v0 = 0.0) {
  // get the size of the vectors
  int vecSize = a.length();

  // initialize a numeric vector "v" (for the result)
  NumericVector v(vecSize);

   // compute v_0
  v[0] = a[0] + b[0] * v0;

  // loop through the vector and compute the new value
  for (int i = 1; i < vecSize; ++i) {
    v[i] = a[i] + b[i] * v[i - 1];
  }
  return v;
}

R で関数をソースして使用するには、次のようにします。

Rcpp::sourceCpp("hydrofun.cpp")

library(data.table)
DT <- data.table(a = 1:4, 
                 b = 0.1,
                 v.ans = c(1, 2.1, 3.21, 4.321))

DT[, v_ans2 := HydroFun(a, b, 0)]
DT
# a   b v.ans v_ans2
# 1: 1 0.1 1.000  1.000
# 2: 2 0.1 2.100  2.100
# 3: 3 0.1 3.210  3.210
# 4: 4 0.1 4.321  4.321

これにより、探している結果が得られます（少なくとも価値の観点から）。

速度を比較すると、約 65 倍の速度アップがわかりました。

library(microbenchmark)
n <- 10000
dt <- data.table(a = 1:n,
                 b = rnorm(n))

microbenchmark(dt[, v1 := vcalc(a, b, 0)],
               dt[, v2 := HydroFun(a, b, 0)])
# Unit: microseconds
# expr                                min        lq       mean    median         uq       max neval
# dt[, `:=`(v1, vcalc(a, b, 0))]    28369.672 30203.398 31883.9872 31651.566 32646.8780 68727.433   100
# dt[, `:=`(v2, HydroFun(a, b, 0))]   381.307   421.697   512.2957   512.717   560.8585  1496.297   100

identical(dt$v1, dt$v2)
# [1] TRUE

それは何らかの形であなたを助けますか？

于 2016-11-03T23:50:28.633 に答える