StatsFragments

Python, R, Rust, 統計, 機械学習とか

{flexclust} + DTW で 時系列を k-means クラスタリングする

R 時系列分析

概要

下の記事のつづき。下の記事では DTW (Dynamic Time Warping) 距離を使って階層的クラスタリングを行った。続けて、 DTW 距離を使って 非階層的クラスタリング (k-means法) を試してみる。

stats::kmeans では任意の距離関数を利用することはできないため、任意の距離関数が利用できる {flexclust} というパッケージを使うことにした。

補足 DTW についてはこちら。

動的時間伸縮法 / DTW (Dynamic Time Warping) を可視化する - StatsFragments

インストール

install.packages('flexclust')
library(flexclust)

{flexclust} の使い方

サンプルデータは前回記事と同一。 まずは既定 (ユークリッド距離) で k-means してみる。サンプルデータは行持ちでわたす必要があるので転置する。

クラスタリングした結果は TSclust::cluster.evaluation を使って正分類率を算出して評価する。そのため、 {TSclust} もロードする。

library(TSclust)

# クラスタ数は 3
res = kcca(t(data), 3)
cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), res@cluster)
# [1] 0.7705387

任意の距離関数で k-means する場合は、kcca の family オプションで、 kccaFamily インスタンス化した距離関数を渡してやればよい。ここでは、既定の方法と同じ flexclust::distEuclidean を距離関数として渡してみる。

res = kcca(t(data), 3, family = kccaFamily(dist = distEuclidean))
cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), res@cluster)
# [1] 0.7705387

補足 ここではたまたま 1回目 / 2回目の正分類率が同じだったが、乱数 / 中心点の初期値によっては正分類率が異なることもありうる。

{flexclust} に任意の距離関数を渡す

DTW を使って k-means するためには、上記の distEuclidean を参考にして距離関数を定義してやればよさそうだ。どのような定義であればよいのか調べるため、flexclust::distEuclidean の実装をみてみる。引数と返り値は、

  • x: 入力データ (次元は データ数 x データの長さ)
  • centers : 現在のクラスタの中心点 (次元は クラスタ数 x データの長さ)
  • z: 各入力データと各クラスタの距離 (次元は データ数 x クラスタ数)
distEucledean
# function (x, centers) 
# {
#     if (ncol(x) != ncol(centers)) 
#         stop(sQuote("x"), " and ", sQuote("centers"), " must have the same number of columns")
#     z <- matrix(0, nrow = nrow(x), ncol = nrow(centers))
#     for (k in 1:nrow(centers)) {
#         z[, k] <- sqrt(colSums((t(x) - centers[k, ])^2))
#     }
#     z
# }
# <environment: namespace:flexclust>

# x として渡ってくるもの (データ数 x データの長さ の行列)
# サンプルの場合、実際は 60行目 / 20列目までデータあり
         [,1]       [,2]       [,3]       [,4]     [,5]      [,6]     [,7]     [,8]      [,9]
U1 -0.6327521  2.8410958  3.4529034 3.76949256 4.306223 5.2242329 1.305026 3.493067 8.0427642
U2 -0.9945734  3.4368117 -0.9907990 2.60606620 1.453351 0.8575501 2.541472 1.869034 0.6731607
U3  0.7014436  0.9778474  0.5219079 3.13321606 3.984465 3.8435608 3.707534 4.780511 7.7889760

# center として渡ってくるもの (クラスタ数 x 系列の長さ の行列)
# サンプルの場合、実際は 20列目までデータあり
           [,1]       [,2]      [,3]      [,4]       [,5]       [,6]       [,7]      [,8]
D3   0.62595969  1.1235099  1.592674 -0.949247 -0.8132157 -2.0273653 -4.4877420 -5.137003
N17 -0.05983237  0.2078591  1.261035  1.826098  1.4504397  2.3426289  3.7686407  3.463603
N4  -1.52153133 -2.0982198 -1.930059 -4.379938 -2.9796999 -0.9948817  0.4933508 -3.084777

# z として返すべきもの (データ数 x クラスタ数 の行列)
# サンプルの場合、実際は 60行目 までデータあり
         [,1]     [,2]     [,3]
[1,] 69.75362 16.51187 44.71353
[2,] 71.45904 18.28978 45.20475
[3,] 57.58485  9.22195 32.93481
[4,] 74.50569 21.34163 48.72438
[5,] 74.53808 20.40498 49.45025
[6,] 66.04604 14.45167 40.60733

形式がわかったので、 TSclust::diss.DTWARP を利用して同様のロジックを作成する。

distDTWARP <- function (x, centers) {
  if (ncol(x) != ncol(centers)) 
    stop(sQuote("x"), " and ", sQuote("centers"), " must have the same number of columns")
  z <- matrix(0, nrow = nrow(x), ncol = nrow(centers))
  for (k in 1:nrow(centers)) {
    z[, k] <- apply(x, 1, function(x) diss.DTWARP(x, centers[k, ]))
  }
  z
}

いざ実行。

res_dtw = kcca(t(data), 3, family = kccaFamily(dist = distDTWARP))

、、、。

、、、。

、、、反応ないな、、、。

動いていることを信じて処理時間を計ってみる。

res_dtw = kcca(t(data), 3, family = kccaFamily(dist = distDTWARP))

#     ユーザ    システム        経過  
#  1142.323     44.924   1306.208 

cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), res_dtw@cluster)
# [1] 0.7653292

なんだこれは、、、。遅いし あんま結果もよくない。

{flexclust} に任意のクラスタ中心点更新関数を渡す

調べてみると、どうも kccaFamily に原因がありそうだ。 kccaFamily には 距離関数 dist のほかに、クラスタ中心点の更新に使う関数 cent が渡せる。ここで、dist には関数を指定 / cent には何も指定しない場合、kccaFamily では以下の関数 centOptim を使ってクラスタ中心点を更新する。DTW は 普通の距離関数に比べると計算コストが高いため、この optim での最適化に非常に時間がかかっていたようだ。

centOptim <- function (x, dist) {
    foo <- function(p) sum(dist(x, matrix(p, nrow = 1)))
    optim(colMeans(x), foo)$par
}

クラスタ中心点の更新処理については ユークリッド距離の最小化 or クラスタの中央値で行うようにしてみる。感覚的には、これらの方法で更新されたクラスタの中心点を使った場合も、クラスタ内の各要素 / 中心点の DTW 距離は小さくなるはずだ。

# ユークリッド距離最小となるように中心点更新
cent <- function(x) centOptim(x, dist = distEuclidean)
res_dtw_ceuc <- kcca(t(data), 3, family = kccaFamily(dist = distDTWARP, cent = cent))

#    ユーザ    システム        経過  
#    3.374      0.118      3.942 

cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), res_dtw_ceuc@cluster)
# [1] 0.7638889

# クラスタの中央値で中心点更新
cent <- function(x) apply(x, 2, median)
res_dtw_cmed <- kcca(t(data), 3, family = kccaFamily(dist = distDTWARP, cent = cent))

#    ユーザ    システム        経過  
#     5.926      0.209      7.039 

cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), res_dtw_cmed@cluster)
# [1] 0.7638889

試してみると、処理時間はかなり短くなった。

比較

k-means 法でのクラスタリングは初期値の影響を受けるため、100回ほどループさせて比較してみた。 比較したのは 以下の3手法

横軸に正分類率 / 縦軸に頻度として頻度分布を描いてみた。ベストな値は DTW (Euclidean) で出ていたが、顕著な差があるようには見えない。ループ回数を増やせば /　元データがちょっと違えば逆転しそうである。

ただ、DTW (Euclidean) の結果はばらつきが少ないようには見える。

f:id:sinhrks:20141121073027p:plain

まとめ

  • {flexclust} を使って DTW 距離で k-means してみた。その際、クラスタ中心点の更新を既定 ( DTW距離 + optim ) で行うのは処理速度の点で現実的でない。クラスタ中心点の更新を ユークリッド距離 or 中央値で行うことで処理速度は改善する。しかし、上記の方法/データではユークリッド距離での k-meansと比較して 劇的な差はみられなかった。

  • 確かめられていないのは、

補足

こちらの本に DTW での k-meansというトピック (時系列ではなく手書き文字についてだが) があるようだ。高い、、、。

Developing Concepts in Applied Intelligence (Studies in Computational Intelligence)

Developing Concepts in Applied Intelligence (Studies in Computational Intelligence)

  • 作者: Kishan G. Mehrotra,Chilukuri Mohan,Jae C. Oh,Pramod K. Varshney,Moonis Ali
  • 出版社/メーカー: Springer
  • 発売日: 2011/06/10
  • メディア: ハードカバー
  • この商品を含むブログを見る

2015/06/03 追記: 直近でブクマいただいていたので補足を。DTW は ユークリッド距離と比べると計算量が多いため、 k-means のように距離計算を複数回 繰り返すアルゴリズムとはあまり相性がよくない。今回、自分はトレンドによって時系列を分類したかったのだが、最終的には上昇/下降など 解釈可能なパターンをいくつか教師データとして用意して 最近傍法を使って分類した (距離計算の回数を減らし、かつ 解釈できるパターンに分類する必要があったため)。