Rcppで不等間隔の時系列データに対して移動平均を計算する

一年ぶりの更新ですね…。

不等間隔の時系列データに対して一定のwindow幅で移動平均や移動分散などを計算したいが、Rでは時間がかかる。という相談が slack の r-wakalang にありました。

推測ですが、Rで普通に（?）やるとwindowの範囲に入っている値を検索するのに、毎回全てのデータから検索してしまうところがボトルネックになっていると思います。そこでRcppを使って値の検索の部分を効率化してみました。

今回は移動平均を計算していますが、紹介する関数を雛形にすれば、windowに対して様々な処理を適用できるはずです。

まずはデータを準備します。

# windowの始点
from <- seq.POSIXt(as.POSIXlt("2017-01-01 00:00:00"), as.POSIXlt("2017-01-31 00:00:00"), by = 60*60)
# windowの終点
to <- from + 3*60*60

# 不等間隔な計測データを適当に再現する
library(dplyr)
time <- seq.POSIXt(min(from), max(to), by = 10) %>% sample(size = 10000) %>% sort
input <- data.frame(
  time,                    #計測時刻
  data = runif(length(time), 0, 5) #計測データ
)

こんな感じのデータ。

> input
                     time         data
1     2017-01-01 00:00:00 4.573989e+00
2     2017-01-01 00:02:40 1.899510e+00
3     2017-01-01 00:06:20 3.148162e+00
4     2017-01-01 00:10:20 1.928675e+00
5     2017-01-01 00:11:30 2.711351e+00
6     2017-01-01 00:13:30 2.464259e+00
7     2017-01-01 00:15:50 7.405586e-01
8     2017-01-01 00:16:40 2.433681e+00
9     2017-01-01 00:19:30 4.152487e+00
...略

Rバージョン

まずは、Rバージョン。工夫すればもっと効率化できるかもしれないですが、とりあえず最初の思いつきで書いたやつ。

# Rで移動平均（デフォルトのwindow幅は3時間）
moving_mean_r <- function(input, from, time_window=3*60*60){
  to <- from + time_window
  mean <- mapply(function(f, t){input %>% filter(time>=f, time<t) %>% '[['("data") %>% mean},from,to)
  data.frame(from, to, mean)
}

上のコードがあまりにも遅いので書き直しました。（ここでdplyr::filter使ってはいけなかったんや…。）

moving_mean_r2 <- function(input, from, time_window=3*60*60){
  to <- from + time_window
  mean <- mapply(function(f, t){input[input$time>=f & input$time<t,"data"] %>% mean;}
                 ,from,to)
  data.frame(from, to, mean)
}

Rcppバージョン

今回は標準C++のアルゴリズムを使っています。

#include <Rcpp.h>
#include <algorithm>

using namespace Rcpp;
//[[Rcpp::plugins(cpp11)]]



// [[Rcpp::export]]
DataFrame moving_mean_rcpp(
    DataFrame input,               //（計測時刻time, 計測値data）のデータフレーム
    newDatetimeVector from,        //移動平均を計算するタイミングの始点ベクトル
    double time_window = 3*60*60)  //移動平均を計算するwindow幅（秒）
{ 
  
  // 計測時刻と計測値をベクトルとして取り出す
  newDatetimeVector time = input["time"]; // 今回は time は昇順にソートされているのが前提です。
  NumericVector     data = input["data"];
  
  // 移動平均を計算するタイミングの終点ベクトル
  newDatetimeVector to = from + time_window;
  
  // 移動平均を計算する数
  R_xlen_t N = from.length();
  
  // 移動平均を格納するベクトル
  NumericVector value(N);
  
  // ベクトル要素の位置をあらわすオブジェクト
  newDatetimeVector::iterator begin = time.begin();
  newDatetimeVector::iterator end   = time.end();
  newDatetimeVector::iterator p1    = begin;
  newDatetimeVector::iterator p2    = begin;
  
  // window i についてループ
  for(R_xlen_t i = 0; i < N; ++i){
    // Rcout << "i=" << i << "\n";
    
    double f = from[i];         //windowの始点の時刻
    double t = f + time_window; //windowの終点の時刻
    
    // windowの終点が最初の計測時刻以前の時はNA、または
    // windowの始点が最後の計測時刻のより後の時はNA
    if(t <= *begin || f > *(end-1)){ 
      value[i]  = NA_REAL;
      continue;//次のループへ
    }
    
    // ベクトル time の位置 p1 以降の要素xから
    // 時刻がwindowの始点f「以降」である「最初の要素」の位置を p1 とする
    p1 = std::find_if(p1, end, [&f](double x){return f<=x;});
    // p1 = std::lower_bound(p1, end, f); //上と同義
    
    // ベクトル time の位置 p1 以降の要素xから
    // 時刻がwindowの終点t「より前」である「最後の要素」の位置を p2 とする
    // （下では、時刻がwindowの終点t「以降」である「最初の要素」の１つ前の位置、にすることで実現している’）
    p2 = std::find_if(p1, end, [&t](double x){return t<=x;}) - 1 ;
    // p2 = std::lower_bound(p1, end, t) - 1 ;//上と同義
    
    // 要素の位置p1,p2を、要素番号i1, i2に変換する
    R_xlen_t i1 = p1 - begin;
    R_xlen_t i2 = p2 - begin; 
    
    
    // 要素番号の確認
    // C++は要素番号が0から始まるのでRに合わせるために1を足している
    // Rcout << "i1 = " << i1+1 << " i2 = " << i2+1 << "\n";
    

    // 該当する範囲のデータについて平均を計算する
    if(i1>i2) {
      value[i] = NA_REAL; // window内にデータがない場合
    } else {
      value[i] = mean(data[seq(i1, i2)]);
    }
    // ↑を変更することで様々なwindow関数を作成できる
    
  }
  
  // 平均を計算した時間の範囲と、平均の値をデータフレームとして出力する
  DataFrame out =
    DataFrame::create(
      Named("from", from),
      Named("to", to),
      Named("mean", value));
  
  return out;
}

RバージョンとRcppバージョンの結果の一致を確認

> # Rバージョン
> res_r <- moving_mean_r(input, from)
> 
> # Rcppバージョン
> res_rcpp <- moving_mean_rcpp(input, from)
> 
> # 結果の比較（数値計算の誤差も考慮して）
> isTRUE(all.equal(res_r$mean, res_rcpp$mean))
[1] TRUE

パフォーマンス比較

> microbenchmark::microbenchmark(
+   moving_mean_r(input, from),
+   moving_mean_r2(input, from),
+   moving_mean_rcpp(input, from)
+ )
Unit: microseconds
                          expr        min          lq        mean      median          uq         max neval
    moving_mean_r(input, from) 904298.672 951065.9675 985553.3977 961292.6805 1015783.925 1179497.694   100
   moving_mean_r2(input, from) 191614.227 208638.5365 226343.9220 216555.8075  230387.330  317652.639   100
 moving_mean_rcpp(input, from)    441.276    488.5045    569.5601    510.9075     546.498    2998.015   100

RとRcppでかなり差が出ました。