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