読者です 読者をやめる 読者になる 読者になる

StatsFragments

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

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

R 時系列分析

概要

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

{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) 
}

そこそこうまくいったものを上位から順に並べる。結果として、(最遠隣法では) DTW が最も正分類率が高かった。

※ 各系列について、クラスタリング後のラベルで色分け / 左上の数値が正分類率。

DTW

f:id:sinhrks:20141116192156p:plain

CID

f:id:sinhrks:20141116192218p:plain

Wavelet変換

f:id:sinhrks:20141116192228p:plain

ユークリッド距離

f:id:sinhrks:20141116192239p:plain

まとめ

{TSclust} パッケージで 階層的クラスタリングをするとき、DTW を使うとユークリッド距離よりもうまくトレンドを拾って分類できている (ように見える)。

クラスタリングする時系列の数が少ない場合は 階層的クラスタリングは使えそう。ただし、ちゃんとやる場合は {TSclust} 側の手法だけでなく hclust の手法との組み合わせについても あてはまりをみるべき。

自分は最終的に より系列数が多い実データにあてたいので、DTW 距離を使って k-means 法でクラスタリングしたい。stats::kmeans では 自作の距離関数を利用できないが、ちょっと調べたところ {flexclust} というパッケージを使えばよさげ。

cluster analysis - How to specify distance metric while for kmeans in R? - Stack Overflow

つづく、、、。

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