{TSclust} ではじめる時系列クラスタリング

概要

こちらで書いた動的時間伸縮法 / DTW (Dynamic Time Warping) を使って時系列をクラスタリングしてみる。ここからはパッケージ {TSclust} を使う

<a href="http://sinhrks.hatenablog.com/entry/2014/11/14/232603">動的時間伸縮法 / DTW (Dynamic Time Warping) を可視化する - StatsFragments</a>

{TSclust} のインストール

install.packages('TSclust')
library(TSclust)

サンプルデータの準備

{TSclust} では時系列間の距離を計算する方法をいくつか定義している。クラスタリングの際にどの定義 (距離) を使えばよいかは時系列を何によって分類したいのかによる。{TSclust} に実装されているものをいくつかあげると、

diss.ACF : ACF
diss.CID : Complexity Correlations (よくわからん) で補正したユークリッド距離
diss.COR : ピアソン相関 (ラグは考慮しない)
diss.EUCL : ユークリッド距離
diss.DTWARP : DTW 距離
diss.DWT : Wavelet変換してユークリッド距離

自分は、"ある時点からのトレンドがどういったパターン (上昇, 下降, etc...) に属するかで分類したい" ので、今回はユークリッド距離か DTW 距離が適切かなと思う。

まず、自分が実際に処理したいものに近そうなサンプルデータを準備した。条件として考慮したのは、

背後にランダムだが一定のトレンドを持つ
かつランダムな AR 構造を持つ

もの。上記のようなデータを以下の3グループ, 各 20 系列作成した。

上昇トレンド : 変数名で group1。ラベルは U で始まる。
停滞トレンド : 変数名で group2。ラベルは N で始まる。
下降トレンド : 変数名で group3。ラベルは D で始まる。

set.seed(1)

# 各グループの系列数
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
}

group1 = generate_ts(TREND, label = 'U')
group2 = generate_ts(0, label = 'N')
group3 = generate_ts(-TREND, label = 'D')

data <- cbind(group1, group2, group3)
data <- as.data.frame(data)

作成した系列 (グループで色分け)

f:id:sinhrks:20141116191700p:plain

階層的クラスタリングの実行と評価

まずは DTW を使って階層的クラスタリングを行い、デンドログラムを描く。手順は、

TSclust::diss で距離行列を求める。
あとは普通の hclust 同様。

# DTW 距離で距離行列を作成
d <- diss(data, "DTWARP")
# hclust は既定で実行 = 最遠隣法
h <- hclust(d)
par(cex=0.6)
plot(h, hang = -1)

結果をみると、左のノードから順に、上昇グループ ("U"), 下降グループ ("D"), 停滞グループ ("N") でけっこう綺麗に分かれている気がする。

f:id:sinhrks:20141116191734p:plain

また、TSclust::cluster.evaluation を使ってクラスタの正分類率が算出できる。DTW を使った場合の正分類率は 86.8 % になった。

補足ここでの正分類率はサンプル作成の際に利用したベースのトレンドを正とした。上のグラフで分かる通り AR 成分によって各グループが混ざってしまっているため、どんな方法を使っても正分類率 100 % になることはないと思うが、ひとつの基準ということで。

# クラスタ数は 3 とする
clusters <- cutree(h, 3)
# 正分類率の算出
cluster.evaluation(rep(c(1, 2, 3), rep(N, 3)), clusters)
# 0.8675466

手法の比較

{TSclust} のいくつかの手法を比較してみる。{TSclust} の実行中に何かエラーが出た手法は除いた。

plot_group <- function(plot.data, cluster = NULL, method = '') {
  library(ggplot2)
  library(tidyr)
  plot.data$index <- 1:SPAN
  plot.data <- gather(plot.data, variable, value, -index)
  if (is.null(cluster)) {
    method = 'Original'
    plot.data$colour <- substr(plot.data$variable, 1, 1)
    plot.data$colour <- factor(plot.data$colour, levels = c('U', 'N', 'D'))
  } else {
    cluster <- as.factor(cluster)
    plot.data$colour <- rep(cluster, rep(SPAN, length(cluster)))
  }

  p <- ggplot(plot.data, aes(x = index, y = value, group = variable, colour = colour)) + 
    geom_line() +
    xlab('') + ylab('') + ggtitle(method)
  if (!is.null(cluster)) {
    cluster.true <- rep(c(1, 2, 3), rep(N, 3))
    rate <- cluster.evaluation(cluster.true, cluster)
    rate <- paste0(round(rate * 100, 1), '%')
    p <- p + annotate(geom = 'text', x = 5, y = 20, label = rate, size = 10)
  }
  print(p)
}

methods <- c("ACF", "AR.LPC.CEPS", "AR.PIC", "CID", "COR", "CORT",
             "DTWARP", "DWT", "EUCL", "INT.PER", "PACF", "PDC",
             "PER", "SPEC.LLR", "SPEC.GLK", "SPEC.ISD")
for (method in methods) {
  print(method)
  d <- diss(data, method)
  h <- hclust(d)
  clusters <- cutree(h, 3)
  plot_group(data, cluster = clusters, method = method) 
}