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)
## [1] 286
#[1] 286

library(dtw)
## Loading required package: proxy
## 
## Attaching package: 'proxy'
## The following objects are masked from 'package:stats':
## 
##     as.dist, dist
## The following object is masked from 'package:base':
## 
##     as.matrix
## Loaded dtw v1.21-3. See ?dtw for help, citation("dtw") for use in publication.
d <- dtw::dtw(ts_a, ts_b, step.pattern = symmetric1)
d$distance
## [1] 286
# 286

Dynamic Time Warping

2つの時系列の破線で結ばれた点同士を順に比較する

比較した2点間の距離 (コスト) を計算する。セルの色が緑色なら、比較した2点間の距離は小さい。(左下)

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

dtwPlotTwoWay(d, ts_a, ts_b)