六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

単純なK-meansと{TSclust}のDTWによる時系列クラスタリングとではどう違うのか実験してみた

これは単なる備忘録です。詳細を知りたいという方は、この記事の元ネタになった以下のid:sinhrksさんの記事をお読みください。



ここでの問題意識は非常にシンプルで「そもそも時系列クラスタリングをかなり膨大な行数のデータに対して実行する際にどれほど厳密にやるべきか?」というお話です。というのも、たかだか数百行ぐらいのデータならリンク先にもあるように動的時間伸縮法(Dynamic Time Warping: DTW)で距離行列を求めてhclust関数で動く階層的クラスタリングでやるのが正しいと思うのですが、これが1万行とか10万行とかになるとそもそもdiss関数やhclust関数が重過ぎて動かないわけです。


そういうケースで、例えば原理的には正しくないのを承知の上で単純なK-means(リンク先にあるような{flexclust} + DTWでやるのではなく)でやっても、それなりの精度が出せるのであれば「仮に」使う目的で代用しても良いのかどうか?というのが気になった、というのが今回の記事のモチベーションです。言い方は悪いですが「K-meansで何とかなるのなら大規模データにはK-meansで代用してしまう」が是か非かを知りたい、ということです。


ということで、面倒なので一気に本題に行きたいと思います。


上昇トレンド・平均回帰・下降トレンド・季節調整あり時系列を用意する


やることは以下のRスクリプトを回すだけです。完全にリンク先の記事を踏襲していますが、季節調整っぽいものを加えるために適当にsinカーブを異なる周期で加えたものを2つ追加しています。

> # 各グループの系列数
> N <- 20
> # 系列の長さ
> SPAN <- 24
> # トレンドが上昇/ 下降する時の平均値
> TREND <- 0.5
> 
> generate_ts <- function(m, label) {
+     library(dplyr)
+     # ランダムな AR 成分を追加
+     .add.ar <- function(x) {
+         x + arima.sim(n = SPAN, list(ar = runif(2, -0.5, 0.5)))
+     }
+     # 平均が偏った 乱数を cumsum してトレンドとする
+     d <- matrix(rnorm(SPAN * N, mean = m, sd = 1), ncol = N) %>%
+         data.frame() %>%
+         cumsum()
+     d <- apply(d, 2, .add.ar) %>%
+         data.frame()
+     colnames(d) <- paste0(label, seq(1, N))
+     d
+ }
> 
> set.seed(101)
> group1 <- generate_ts(TREND, label = 'U')
> set.seed(102)
> group2 <- generate_ts(0, label = 'N')
> set.seed(103)
> group3 <- generate_ts(-TREND, label = 'D')
> set.seed(104)
> group4 <- generate_ts(0, label = 'S1_')
> group4 <- group4 + 5*sin(seq(0, 4*pi, length.out = SPAN))
> set.seed(105)
> group5 <- generate_ts(0, label = 'S2_')
> group5 <- group5 + 5*sin(seq(0, 8*pi, length.out = SPAN))
> 
> data <- cbind(group1, group2, group3, group4, group5)
> data <- as.data.frame(data)
> matplot(data, type = 'l', lty = 1, col = c(rep(1,20), rep(2,20), rep(3,20), rep(4,20), rep(5,20)))
> true.cluster <- c(rep(1,20), rep(2,20), rep(3,20), rep(4,20), rep(5,20))

f:id:TJO:20190118002610p:plain

こんな感じの時系列になると思います。


単純なK-meansでやってみる


横持ちのデータを縦持ちに直せば、普通にkmeans関数で回せます。

> model.km <- kmeans(t(data), centers = 5)
> model.km$cluster
   U1    U2    U3    U4    U5    U6    U7    U8    U9   U10   U11   U12   U13 
    1     3     1     3     3     4     1     3     3     3     4     1     1 
  U14   U15   U16   U17   U18   U19   U20    N1    N2    N3    N4    N5    N6 
    1     3     1     1     1     3     4     1     1     4     5     5     5 
   N7    N8    N9   N10   N11   N12   N13   N14   N15   N16   N17   N18   N19 
    1     1     5     2     5     5     4     5     5     4     1     4     5 
  N20    D1    D2    D3    D4    D5    D6    D7    D8    D9   D10   D11   D12 
    4     2     2     2     5     2     2     2     2     2     2     5     5 
  D13   D14   D15   D16   D17   D18   D19   D20  S1_1  S1_2  S1_3  S1_4  S1_5 
    2     2     2     2     5     2     5     5     4     5     4     4     1 
 S1_6  S1_7  S1_8  S1_9 S1_10 S1_11 S1_12 S1_13 S1_14 S1_15 S1_16 S1_17 S1_18 
    5     5     5     5     4     5     5     5     1     5     5     5     5 
S1_19 S1_20  S2_1  S2_2  S2_3  S2_4  S2_5  S2_6  S2_7  S2_8  S2_9 S2_10 S2_11 
    4     2     5     4     4     4     4     1     4     4     2     4     4 
S2_12 S2_13 S2_14 S2_15 S2_16 S2_17 S2_18 S2_19 S2_20 
    4     2     4     4     5     5     5     5     4 
> table(true.cluster, model.km$cluster)
            
true.cluster  1  2  3  4  5
           1  9  0  8  3  0
           2  5  1  0  5  9
           3  0 14  0  0  6
           4  2  1  0  5 12
           5  1  2  0 12  5
> matplot(data, type = 'l', lty = 1, col = model.km$cluster)

f:id:TJO:20190118152052p:plain

パッと見ではそれほど悪くないように見えるんですが、confusion matrix的なものを出して「それぞれの行の最大値」を正解の場合の数とみなして計算すると、正答率は56%に留まります


{TSclust}のDTWで階層的クラスタリングをやってみる


これは完全にリンク先を踏襲しています。

> # DTW 距離で距離行列を作成
> d <- diss(data, "DTWARP")
> # hclust は既定で実行 = 最遠隣法
> h <- hclust(d)
> par(cex=0.6)
> plot(h, hang = -1)
> h.cluster <- cutree(h, 5)
> h.cluster
   U1    U2    U3    U4    U5    U6    U7    U8    U9   U10   U11   U12   U13 
    1     1     1     1     1     2     3     1     1     1     2     3     1 
  U14   U15   U16   U17   U18   U19   U20    N1    N2    N3    N4    N5    N6 
    1     1     1     3     1     1     2     3     1     2     4     4     4 
   N7    N8    N9   N10   N11   N12   N13   N14   N15   N16   N17   N18   N19 
    3     3     2     4     4     4     2     4     4     2     3     2     4 
  N20    D1    D2    D3    D4    D5    D6    D7    D8    D9   D10   D11   D12 
    2     5     4     4     4     4     4     4     5     5     5     4     2 
  D13   D14   D15   D16   D17   D18   D19   D20  S1_1  S1_2  S1_3  S1_4  S1_5 
    5     5     4     5     4     5     4     4     2     4     2     2     3 
 S1_6  S1_7  S1_8  S1_9 S1_10 S1_11 S1_12 S1_13 S1_14 S1_15 S1_16 S1_17 S1_18 
    4     4     4     4     2     4     4     4     3     4     4     2     4 
S1_19 S1_20  S2_1  S2_2  S2_3  S2_4  S2_5  S2_6  S2_7  S2_8  S2_9 S2_10 S2_11 
    2     4     4     2     2     4     2     3     2     4     4     3     4 
S2_12 S2_13 S2_14 S2_15 S2_16 S2_17 S2_18 S2_19 S2_20 
    4     4     2     2     4     4     4     4     2 
> table(true.cluster, h.cluster)
            h.cluster
true.cluster  1  2  3  4  5
           1 14  3  3  0  0
           2  1  6  4  9  0
           3  0  1  0 11  8
           4  0  6  2 12  0
           5  0  7  2 11  0
> matplot(data, type = 'l', lty = 1, col = h.cluster)
> cluster.evaluation(rep(c(1, 2, 3, 4, 5), rep(N, 5)), h.cluster)
[1] 0.4774603

f:id:TJO:20190118003043p:plain
f:id:TJO:20190118152337p:plain

これも見た目にはそれほど悪くないような気がしますが、{TSclust}のcluster.evaluationを使うと正答率は47.7%。先ほどの理屈でconfusion matrixから計算しても57%という結果になりました。


現時点での感想


ただのK-meansを使っても{TSclust} + DTWでやってもそれほど変わらないかもしれない、という感想を持ちました。とは言え、この感想が妥当だとは正直思えない(理論的には当然このやり方は正しくない)ので、もう少し色々調べてみようと思います。。。