【clustering】R と Python で k-means を試してみた【PCA】

| 0件のコメント

続・わかりやすいパターン認識―教師なし学習入門―を読んでいる。今回は第10章 クラスタリングの中の k-meansについて。

クラスタリング

クラスタリングはデータをいくつかのグループに分ける手法。ラベル付与を自動化して省力化できる。また, 類似する文書をグループ化し整理するテキストクラスタリングがある。
ビジネス的な利用としては, 顧客の分類などマーケティング活動で使われる。

k-means

k-meansは, 既知のクラスタ数を cとして, 幾つかの距離 (ユーグリッド距離 etc)に基づいてクラスタリングを行う。

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は特徴空間上でクラスタごとに集中してまとまっているほど小さくなる。

計算手順は下記となる。(間違っていたらすみません)

  • STEP1 : 初期状態としてd次元特徴空間上に各クラスタに対応する c個のプロトタイプ [2]を設定する。
  • STEP2 : 全てのパターンに対して所属クラスタの割り当てを行う。各Piを固定してeをδkについて最小化する。
  • STEP3 : プロトタイプの更新を行う。各δkを固定してeをPiについて最小化する。
  • STEP4 : クラスタの再割当てが発生しなければ収束と判定する。そうでない場合は STEP2に戻る。

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 に置いた。[4]
今回用いた 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した。

kmeans-wine

最適なクラスの決定の指標としてハーティンガンルールを使う方法がある。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

hartigans-rule

量子化誤差の収束過程

k-meansを使うだけなら, scikit-learnなどに実装があるが, Rのkmeans()で得られる objectからは量子化誤差を取得できなかったので, Pythonで k-meansを書いて GitHub に置いた。 現在の量子化誤差と修正したプロトタイプにおける量子化誤差の差が閾値未満の場合に収束と判定した。

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)

wine-kmeans-init

収束した状態 (Iter=5) で Rとほぼ同じ結果になった。この時の量子化誤差は約260。

wine-kmeans-convergence

収束までの過程を plot してみた。

Quantization-Error

おわりに

1985年頃から現在まで株価が記録されている企業約900社の株価の変動を k-means でクラスタリングしてみたが, 業種ごとの差があまり見えず面白い結果にはならなかった。そもそも, k-means では相関の高いパターンを多く含むとクラスタリングが難しくなるのかもしれない。今回は外れ値として かわでんと NTTの2社が浮かび上がった。

chart-6648

chart-9432

かわでんは2000年に民事再生法適用を申請したため株価が急激に変動した時期があり, NTTはバブル期に上場したため 318万円の高値をつけたこともあるが, 現在(2015/10/02)は 4,213円となっているため, 外れ値となったのかもしれない。他の視点でも分析してみて, 面白い結果になったらまとめようかなと思う。


[1] ベクトル量子化の際に伴う誤差を量子化誤差という。
[2] セントロイド (centroid)とも言う。
[3] 凸クラスタリングは, クラスタ数cを大きな値に設定しても, 最適化の過程で不要な分布はその事前確率が次第に0に近づいていく除去されるという考えに基づいている。
[4] R-Study
[5] みんなのR -データ分析と統計解析の新しい教科書-
[6] Hartigan-Wong のアルゴリズムを確認する
[7] パッケージユーザーのための機械学習(7):K-meansクラスタリング
[8] クラスタリングの不可能性定理について - Qiita

コメントを残す

必須欄は * がついています