いま手元に 20万件くらいの時系列があって、それらを適当にクラスタリングしたい。どうしたもんかなあ、と調べていたら {TSclust} というまさになパッケージがあることを知った。

このパッケージでは時系列の類似度を測るためのさまざまな手法 (=クラスタリングのための距離) を定義している。うちいくつかの手法を確認し、動的時間伸縮法 / DTW (Dynamic Time Warping) を試してみることにした。

DTWの概要

時系列相関 (CCF) の場合は 片方を 並行移動させているだけなので 2つの系列の周期が異なる場合は 相関はでにくい。

DTW では 2つの時系列の各点の距離を総当りで比較した上で、系列同士の距離が最短となるパスを見つける。これが DTW 距離 になる。そのため、2つの系列の周期性が違っても / 長さが違っても DTW 距離を定義することができる。

アルゴリズム

英語版 Wikipedia がわかりやすい。{TSculst} とは別で DTW を計算できる {dtw} パッケージもあるが、アルゴリズムがシンプルなのでここでは 直接 実装してみた。

dtw_distance <- function(ts_a, ts_b, d = function(x, y) abs(x-y),
                         window = max(length(ts_a), length(ts_b))) {
  ts_a_len <- length(ts_a)
  ts_b_len <- length(ts_b)
  
  # コスト行列 (ts_a と ts_b のある2点間の距離を保存)
  cost <- matrix(NA, nrow = ts_a_len, ncol = ts_b_len)
  
  # 距離行列 (ts_a と ts_b の最短距離を保存)
  dist <- matrix(NA, nrow = ts_a_len, ncol = ts_b_len)
  
  cost[1, 1] <- d(ts_a[1], ts_b[1])
  dist[1, 1] <- cost[1, 1]

  for (i in 2:ts_a_len) {
    cost[i, 1] <- d(ts_a[i], ts_b[1])
    dist[i, 1] <- dist[i-1, 1] + cost[i, 1]
  }
  
  for (j in 2:ts_b_len) {
    cost[1, j] <- d(ts_a[1], ts_b[j])
    dist[1, j] <- dist[1, j-1] + cost[1, j]
  }
  
  for (i in 2:ts_a_len) {
    # 最短距離を探索する範囲 (ウィンドウサイズ = ラグ)
    window.start <- max(2, i - window)
    window.end <- min(ts_b_len, i + window)
    for (j in window.start:window.end) {
      # dtw::symmetric1 と同じパターン
      choices <- c(dist[i-1, j], dist[i, j-1], dist[i-1, j-1])
      cost[i, j] <- d(ts_a[i], ts_b[j])
      dist[i, j] <- min(choices) + cost[i, j]
    }
  } 
  return(dist[nrow(dist), ncol(dist)])
}

ts_a <- AirPassengers[31:45]
ts_b <- AirPassengers[41:55]

dtw_distance(ts_a, ts_b)
# 286

{dtw} パッケージの計算結果と一致することを確認。

library(dtw)
d <- dtw::dtw(ts_a, ts_b, step.pattern = symmetric1)
d$distance
# 286

可視化

{animation} を含むコードは gist においた (重いです)。

Dynamic Time Warping

• 2つの時系列の破線で結ばれた点同士を順に比較する
• 比較した2点間の距離 (コスト) を計算する。セルの色が緑色なら、比較した2点間の距離は小さい。(左下)
• その際、左側、左下、下にある距離行列のセル + 上で計算した2点間の距離 (コスト) を足して、最も小さい値をその時点の DTW 距離にする (そこまでの2点間の距離 (コスト) の和が最小になるようなパスが自動的に見つかる)。セルの色が青色なら、そのセルまでの DTW 距離は小さい。(右下)
• 距離行列の右上のセルに到達したとき、そのセルの値が 系列同士の DTW 距離になっている

{dtw} の dtwPlotTwoWay 関数を使うと DTW で 2 つの系列の各点がどのようにマッピングされたのかプロットできる。

dtwPlotTwoWay(d, ts_a, ts_b)

{dtw} の使い方

dtw では以下のようなオプションが指定でき、より柔軟に DTW 距離が求められる。

• dist.method : 2点間の距離(コスト)を求める関数
• step.pattern : DTW距離行列にコストを追加する際の関係式。
• window.type : window の種類
• open.begin, open.end : 最初 / 最後の要素のうちマッチしないものを捨てるかどうか

step.pattern なんかはいろいろあって奥が深そう。

# step.pattern のうち、上のサンプルと同じ動きのもの
symmetric1
# Step pattern recursion:
# g[i,j] = min(
#      g[i-1,j-1] +     d[i  ,j  ] ,
#      g[i  ,j-1] +     d[i  ,j  ] ,
#      g[i-1,j  ] +     d[i  ,j  ] ,
#   )
# 
#  Normalization hint: NA

# べつのパターン
asymmetric
# Step pattern recursion:
# g[i,j] = min(
#      g[i-1,j  ] +     d[i  ,j  ] ,
#      g[i-1,j-1] +     d[i  ,j  ] ,
#      g[i-1,j-2] +     d[i  ,j  ] ,
#   )
# 
#  Normalization hint: N