『続・わかりやすいパターン認識―教師なし学習入門―
クラスタリング
クラスタリングはデータをいくつかのグループに分ける手法で教師なし学習に分類される。類似する文書をグループ化し整理するテキストクラスタリングや DTW 距離を用いた時系列クラスタリングなどがある。
代表的なクラスタリング手法として以下がある。
- K-Means Clustering
- Mean-Shift Clustering
- Density-Based Spatial Clustering of Applications with Noise (DBSCAN)
- Expectation–Maximization (EM) Clustering using Gaussian Mixture Models (GMM)
- Agglomerative Hierarchical Clustering
k-Means
k-Meansは, 既知のクラスタ数を c として, 幾つかの距離 (e.g. ユーグリッド距離)に基づいてクラスタリングを行う。
d次元特徴空間に個々のパターン x1,x2…,xk が分布しており, 各パターンは異なる c種のいずれかのクラスに所属していることとする。
パターン xk に対して, c 次元ベクトル δk を以下に定義する。
δk = (0,0,…,0,1,0…,0)
δk のi番目のみ1として他は全て 0 とすることで, パターン xk をクラスタ ωi に割り当てたこととし, クラスタ ωi に属する全てのパターンを1個のプロトタイプ Pi に代表させる。
この処理は, 特徴空間上のパターンの分布を c 個のベクトル, P1, P2,…Pc で量子化 (Quantization) したとみなせるのでベクトル量子化という。
クラスタリングの妥当性を考える上で量子化誤差 e [1]を使って評価を行うことができる。e は特徴空間上でクラスタごとに集中してまとまっているほど小さくなる。
計算手順は以下となる。
- 初期状態としてd次元特徴空間上に各クラスタに対応する c個のプロトタイプ [2]を設定する。
- 全てのパターンに対して所属クラスタの割り当てを行う。各Piを固定してeをδkについて最小化する。
- プロトタイプの更新を行う。各δkを固定してeをPiについて最小化する。
- クラスタの再割当てが発生しなければ収束と判定する。そうでない場合は 2 に戻る。
k-means は教師なし学習として紹介されることが多いが, STEP3でプロトタイプの更新を行っている点は自ら教師信号を生成しながら学習を行っていると考える見方もある。
また, k-means は混合正規分布 (normal mixture distribution) のパラメータ推定の特別な場合と捉えることができる。
欠点としては下記が挙げられる。
- クラスタ数をあらかじめ設定しないとならない。従って, 最適なクラスタ数を探索する必要がある。
- 初期値の設定が不完全な場合, 大域的最適解ではなくて局所的最適解に収束してしまう。
上記の初期値の選択問題の改良を行った手法としては k-means++ があり, 実装では R の {LICORS} 等がある。
凸クラスタリング
上記の問題点の解決策のひとつとして, クラスタ数を推定できる 凸クラスタリング (Convex Clustering) [3]を用いる方法もある。
凸クラスタリングは混合正規分布のパラメータ推定を基礎にしている。
k-meansはハードな割り当てという各パターンにひとつのクラスタを割り当てていく手法であるが, 混合正規分布のパラメータ推定ならびに凸クラスタリングは複数のクラスタに確率的に割り当てられるためソフトな割り当てと呼ばれる。
詳細なアルゴリズムは本書を参考に。
凸クラスタリングは混合正規分布のパラメータ推定に加えて, 下記2つの条件を付加する。
- 正規分布の共分散行列は一定として, かつクラスタで共通とする。
- 各クラスタのプロトタイプはパターン集合の中から選ぶ。
従って, パターンが十分に多くて個々のクラスタが密であれば真のプロトタイプはほぼ一致する。
この手法ではクラスタのサイズを制御する分散パラメータの設定に恣意性が残るが, 与えられた分散パラメータに対して常に最適解が得られるメリットがある。
一方で, 現実には上記の制約条件が成り立たないことがあるので, この点でEMアルゴリズムに比べて不利である。この EMアルゴリズムのベイズ版としては, 変分ベイズ (Variational Bayesian) がある。
クラスタ数 c が未知のまま推定する手法としてはノンパラメトリックベイズモデルが有効で, 11章の話へと繋がっていく。
R で k-means
Rの Code は GitHub に置いた。今回用いた Wineデータセットは3つの異なる品種のワインを化学的な分析により, 13の特徴量で線形分離されたデータセット。
Rでは, stats::kmeans() で簡単に k-means が使える。cluster数に 3を設定。
wine.km <- kmeans(wine.nonlabeled, centers = 3)
summary(as.factor(wine.km$cluster))
table(wine.class, wine.km$cluster)
tableで分類結果を確認する。
wine.class 1 2 3
1 59 0 0
2 3 3 65
3 0 48 0
可視化にあたり, d次元特徴空間に対して主成分分析で次元削減を行い寄与度の高い2つを x, y に設定し plot した。
最適なクラスの決定の指標としてハーティガンルールを使う方法がある。Hartigan が 10 を超える場合に, k+1 を使うことに意味があると解釈する。[5]
今回は 5 と判定された。
Clusters Hartigan AddCluster
1 2 69.523332 TRUE
2 3 52.151051 TRUE
3 4 15.207285 TRUE
4 5 11.604607 TRUE
5 6 9.896997 FALSE
6 7 9.771937 FALSE
7 8 10.969714 TRUE
8 9 7.562533 FALSE
9 10 8.752139 FALSE
10 11 8.474434 FALSE
11 12 6.283500 FALSE
12 13 5.780901 FALSE
13 14 5.707151 FALSE
14 15 4.512510 FALSE
15 16 4.813462 FALSE
16 17 4.508877 FALSE
17 18 5.790036 FALSE
18 19 5.384525 FALSE
19 20 2.358511 FALSE
量子化誤差の収束過程
R の kmeans() で得られるオブジェクトからは量子化誤差が取得できなかったので, Python で量子化誤差を計算する k-means を書いて GitHub に置いた。現在の量子化誤差と修正したプロトタイプにおける量子化誤差の差が閾値未満の場合に収束と判定した。ちなみに Python で k-means を使うだけなら scikit-learn など多くの実装が既にある。
R と同様に Wine データセットを用いた。
$ python kmeans.py
Iter : 0 Quantization Error : 601.122433764 Diff : 871.573005395
Iter : 1 Quantization Error : 327.925733963 Diff : 273.196699801
Iter : 2 Quantization Error : 270.482767355 Diff : 57.4429666081
Iter : 3 Quantization Error : 261.264711188 Diff : 9.21805616624
Iter : 4 Quantization Error : 260.016662684 Diff : 1.24804850477
Iter : 5 Quantization Error : 260.016662684 Diff : 0.0
初期状態 (Iter=0)
収束した状態 (Iter=5) で Rとほぼ同じ結果になった。この時の量子化誤差は約260。
収束までの過程を plot してみた。
Agglomerative Hierarchical Clustering
凝集型階層的クラスタリング (Agglomerative Hierarchical Clustering; AHC) は, N個のデータで最も距離が近い2つのクラスタをマージし, 距離行列の更新する手順を単一のクラスタとなるまで繰り返す単純なアルゴリズム。
大規模データにはあまり向いていないが, 以下の利点がある。
- クラスタ数を事前に仮定しない。
- dendogram を切断することで任意の数のクラスタが得られる。
- K-means と比較して解釈性に優れる。
Iris を例に階層的クラスタリングを R で実行してみる。
library(dplyr)
library(forcats)
library(dendextend)
# change working directory
frame_files <- lapply(sys.frames(), function(x) x$ofile)
frame_files <- Filter(Negate(is.null), frame_files)
setwd(dirname(frame_files[[length(frame_files)]]))
feat <- iris[,-5]
label <- iris[,5]
# Agglomerative Hierarchical Clustering
feat_dist <- dist(feat, method = "euclidean")
hc <- hclust(feat_dist, method = "ward.D2")
dend <- hc %>%
as.dendrogram() %>%
color_branches(., k = 3, groupLabels = levels(label))
plot(dend, xlab = "", ylab = "", sub = "")
# cutting tree
cluster <- cutree(hc, k = 3) %>%
as.factor() %>%
fct_recode("setosa" = "1", "versicolor" = "2", "virginica" = "3")
# confusion matrix
table(cluster, label)
# label
# cluster setosa versicolor virginica
# setosa 50 0 0
# versicolor 0 50 14
# virginica 0 0 36
おわりに
1985年頃から現在まで株価が記録されている企業約900社の株価の変動を k-means でクラスタリングしてみたが, 業種ごとの差があまり見えず面白い結果にはならなかった。そもそも, k-means では相関の高いパターンを多く含むとクラスタリングが難しくなるのかもしれない。今回は外れ値として かわでんと NTTの2社が浮かび上がった。
かわでんは2000年に民事再生法適用を申請したため株価が急激に変動した時期があり, NTTはバブル期に上場したため 318万円の高値をつけたこともあるが, 現在(2015-10-02) は 4,213円となっているため, 外れ値となったのかもしれない。
[1] ベクトル量子化の際に伴う誤差を量子化誤差という。
[2] セントロイド (centroid) とも言う。
[3] 凸クラスタリングは, クラスタ数cを大きな値に設定しても, 最適化の過程で不要な分布はその事前確率が次第に0に近づいていく除去されるという考えに基づいている。
[4] The 5 Clustering Algorithms Data Scientists Need to Know
[5] みんなのR -データ分析と統計解析の新しい教科書-
[6] Hartigan-Wong のアルゴリズムを確認する
[7] パッケージユーザーのための機械学習(7)：K-meansクラスタリング
[8] クラスタリングの不可能性定理について – Qiita
[9] クラスタ数を自動推定するX-means法を調べてみた
[10] Hierarchical Clustering
[11] Hierarchical cluster analysis on famous data sets – enhanced with the dendextend package
