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

六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

パッケージユーザーのための機械学習(7):K-meansクラスタリング

本シリーズ記事のカテゴリからPythonが消えて久しい今日この頃ですが、皆様いかがお過ごしでしょうか*1。とかいう前口上はどうでも良くて、とっとと今回のお題に入りましょう。今回はクラスタリングのド定番、K-means(k平均)クラスタリングです。


K-meansはどの本にも載っている手法なので参考文献には何を挙げても良い気がしますが、とりあえずこれまたド定番ということで相変わらずはじパタを挙げておきます。


はじめてのパターン認識

はじめてのパターン認識


ちなみにPRMLでも下巻のp.140-146で触れてて*2、K-meansアルゴリズムEMアルゴリズムとの類似点についてもコメントがあります。この辺はまた後で触れることにしましょう。


まずRでどんなものか見てみる


ということで、これまたド定番ですがいつものGitHubに置いてある多変量データを使って試してみましょう。これまたこれまたいつも通りdとかいう名前で置いといて、cvカラムをfactor型に直しておきます。


K-meansクラスタリング自体はkmeans(){stats}でやれるので、特にパッケージのインストールなどは要りません。めっちゃ簡単に実行できます。

> head(d[sample(3000,3000,replace=F),])
     a1 a2 a3 a4 a5 a6 a7  cv
1970  1  1  1  1  1  1  1 Yes
1929  0  1  0  1  1  1  1 Yes
434   0  0  0  1  1  0  1  No
2413  0  0  0  1  1  1  0 Yes
1114  0  0  1  1  0  0  1  No
284   1  1  0  1  0  0  1  No

# データは相変わらずこんな感じ

d$cv<-as.factor(d$cv)
# cvカラムをfactor型に直す

> d.km<-kmeans(d[,-8],centers=2)
# 特徴ベクトルのレコードと、centers引数でクラスタ数Kを与えておく

> summary(as.factor(d.km$cluster))
   1    2 
1393 1607 
# d.km$clusterにクラスタのインデックスが格納されている

> table(d$cv,d.km$cluster)
     
         1    2
  No    57 1443
  Yes 1336  164

# 試しにcvカラムのラベルと比較してみたら、
# 正答率(本当は意味ないんだけど)が92.6%になった
# ※ラベルは毎回変動するのでどれかしら塊が合っていればOK

> d.km4<-kmeans(d[,-8],centers=4)
# ちょっとイタズラでクラスタ数4でやってみた

> table(d$cv,d.km4$cluster)
     
        1   2   3   4
  No   27 408 826 239
  Yes 875 282 131 212

# めちゃくちゃな結果になった(1,3クラスタは合ってるっぽいけど)


それなりに綺麗な結果になったんじゃないかと思います。こんな感じで、事前にクラスタ数が分かっていればそれなりに妥当そうな結果になるし、分かっていなければそれなりにめちゃくちゃな結果になるわけです。


K-meansクラスタリングとは何ぞや


はじパタの解説と、Wikipedia英語版の図説が非常に分かりやすいので、基本的にはそちらを見てもらった方が早いでしょう。はじパタp.155には以下のように定義が書かれています。

K-平均法(K-means法)の目的は、d次元のN個のデータ{\cal D} = \{ \mathbf{x}_1, \cdots , \mathbf{x}_N \}, \mathbf{x}_i \in {\cal M}^dを、データ間の類似性(距離)を尺度に、あらかじめ定めたK個のクラスタに分類することである。


前回の階層的クラスタリングと比べてみると、

  1. データ間距離を尺度に決めている点で同じ
  2. クラスタ個数Kを事前に固定する点が異なる


と言ったところでしょうか。階層的クラスタリングは深さ方向(クラスタ個数が増える方向)にどんどん行ってしまう*3のでデータ容量が大きくなると計算量が爆発して簡単にコケてしまいますが、K-meansはクラスタ個数が固定されているので大容量データでも動く点がメリットと言って良いと思います*4。そのアルゴリズムは同じくはじパタp.155によればこんな感じです。

問題設定


クラスタの代表ベクトルの集合を {\cal M} = \{ \mathbf{\mu}_1,\cdots,\mathbf{\mu}_k \}とする。k番目の代表ベクトルが支配するクラスタ(ボロノイ領域になる)を M(\mathbf{\mu}_k)とする。i番目のデータが、k番目の代表ベクトルが支配するクラスタ M(\mathbf{\mu}_k)に帰属するか否かを表す帰属変数q_{ik}を、


q_{ik} = \left\{ \begin{array}{l} 1 (\mathbf{x}_i \in M(\mathbf{\mu}_k)) \\ 0 (else) \end{array} \right. (10.2)


で定義する。K-平均法の評価関数を、


 J(q_{ik},\mathbf{\mu}_k) = \displaystyle{\sum^{N}_{i=1}} \displaystyle{\sum^{K}_{k=1}} q_{ik} | | \mathbf{x}_i - \mathbf{\mu}_k | |^2(10.3)


と定義する。この評価関数には、q_{ik}\mathbf{\mu}_kに関する最適化が含まれている。\mathbf{\mu}_kに関する最適化は、


 \frac{\partial J(q_{ik},\mathbf{\mu}_k)}{\partial \mathbf{\mu}_k} = 2 \displaystyle{\sum^{N}_{i=1}} q_{ik} (\mathbf{x}_i - \mathbf{\mu}_k) = 0

 \Rightarrow \mathbf{\mu}_k = \frac{\displaystyle{\sum^{N}_{i=1}}q_{ik}\mathbf{x}_i}{\displaystyle{\sum^{N}_{i=1}q_{ik}}}(10.4)


に従って行う。すなわち、各クラスタの代表ベクトルは、そのクラスタに帰属するデータの平均ベクトルとなる。しかし、 q_{ik}と同時に最適化するのは難しいので、次のような手順で逐次最適化する。


アルゴリズム


初期化:N個のデータをランダムにK個のクラスタに振り分け、それぞれのクラスタの平均ベクトル(重心、セントロイド)を求め、\mathbf{\mu}_k (k=1,\cdots,K)とする。


(1) q_{ik}に関する最適化:\mathbf{\mu}_kを固定し、帰属変数q_{ik}を、

q_{ik} = \left\{ \begin{array}{l} 1 (k = arg min_j || \mathbf{x}_i - \mathbf{\mu}_j ||^2) \\ 0 (else) \end{array} \right.

に従って決める。


(2) \mathbf{\mu}_kの最適化:q_{ik}を固定し、式(10.4)に従ってセントロイド\mathbf{\mu}_kを決める。


(3) 繰り返し:上記(1)と(2)を収束するまで(状態変化がなくなるまで)繰り返す。


何となく雰囲気は分かると思うんですが、Wikipedia英語版にこのプロセスを分かりやすく図説にしたものがあるので貼っておきます。

f:id:TJO:20140228124150p:plain

1) k initial "means" (in this case k=3) are randomly generated within the data domain (shown in color).


f:id:TJO:20140228124215p:plain

2) k clusters are created by associating every observation with the nearest mean. The partitions here represent the Voronoi diagram generated by the means.


f:id:TJO:20140228124232p:plain

3) The centroid of each of the k clusters becomes the new mean.


f:id:TJO:20140228124247p:plain

4) Steps 2 and 3 are repeated until convergence has been reached.


(wikipedia:en:K-means_clustering)


こんな感じに、2ステップある逐次最適化*5を繰り返すことで最適なクラスタ振り分けを決める、というのがK-meansクラスタリングです。ちなみにこのK-meansにおける「2ステップを交互に繰り返す」というのは今後紹介する予定のEMアルゴリズムと実は全く同じ方法論で、PRMLのp.140にもその旨がコメントされているんですねー。覚えておいて損はないと思います。


注意点が3つあって、1つはこの手の逐次最適化手法は初期値に依存するので、最適解に近い値が欲しい場合は初期値を変えながらヒューリスティックに繰り返す必要があるという点*6。もう1つは、代表ベクトルはあくまでも平均ベクトルであって、実際のデータベクトルとは一般に一致しないという点。ちなみにこれはK-medoidsクラスタリングという、データ間の非類似度行列が定義できる場合に距離の1乗(K-meansは距離の2乗)で誤差評価するような関連手法を使うことで解決できます*7


最後に、K-meansはデータ構造的にデータベクトルの平均ベクトルが定義できない場合は使えないという点*8。これは結構厄介らしく、実は僕も時々仕事でぶつかることがあって頭を悩ませることが多く。。。このシリーズの後の方で良い解決策にたどりつけたらいいなぁ、と(汗)。


真の分布の分かっているデータをクラスタリングしてみる


そんなわけで、これまで通りXORパターンを使ってK-meansをやってみましょう。いつものようにGitHubからシンプルパターン複雑パターンのデータを取ってきて、それぞれxors, xorcという名前でインポートしておきます。これら2つのデータの真の分布は、

# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> xors<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(xors)<-c("x","y","label")

# 複雑なXORパターン
> p21<-cbind(rnorm(n=25,mean=1,sd=1),rnorm(n=25,mean=1,sd=1))
> p22<-cbind(rnorm(n=25,mean=-1,sd=1),rnorm(n=25,mean=1,sd=1))
> p23<-cbind(rnorm(n=25,mean=-1,sd=1),rnorm(n=25,mean=-1,sd=1))
> p24<-cbind(rnorm(n=25,mean=1,sd=1),rnorm(n=25,mean=-1,sd=1))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> xorc<-as.data.frame(cbind(rbind(p21,p23,p22,p24),t))
> names(xorc)<-c("x","y","label")


と表されるような、(1,1), (-1,1), (-1,-1), (1,-1)のそれぞれを中心とする標準偏差0.5ないし1.0の二次元正規分布が4つ配置されたデータです。プロットすると以下にようになるはず。


f:id:TJO:20140228144250p:plain

f:id:TJO:20140228144300p:plain


という点を踏まえて、実際にK-meansクラスタリングをやってみましょう。

> xors.km<-kmeans(xors,centers=4)
> xorc.km<-kmeans(xorc,centers=4)

> table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xors.km$cluster)
   
     1  2  3  4
  1  0 25  0  0
  2  0  0  0 25
  3 25  0  0  0
  4  0  0 25  0

# 正答率100%

> table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xorc.km$cluster)
   
     1  2  3  4
  1  1  3 20  1
  2 21  2  0  2
  3  2  1  0 22
  4  4 21  0  0

# 正答率84%
# ※これもラベルの順番がばらけているので要注意

> plot(xors[,-3],pch=19,cex=3,col=xors.km$cluster,xlim=c(-3,3),ylim=c(-3,3))
> plot(xorc[,-3],pch=19,cex=3,col=xorc.km$cluster,xlim=c(-3,3),ylim=c(-3,3))


f:id:TJO:20140228144312p:plain

K-meansのシンプルパターンと、

f:id:TJO:20140228144332p:plain

複雑パターンの結果。上の方の「正解」のプロットと見比べてもらえば分かるかと思いますが、大体問題なくクラスタリングできてる感じです。


ところで、Rではpam(){cluster}でK-medoidsクラスタリングを行うことができます。全く同じことを今度はK-medoidsでやってみましょう。

> require("cluster")
Loading required package: cluster
> xors.pam<-pam(xors,k=4)
> xorc.pam<-pam(xorc,k=4)

> table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xors.pam$clustering)
   
     1  2  3  4
  1 25  0  0  0
  2  0 25  0  0
  3  0  0 25  0
  4  0  0  0 25

# 正答率100%

> table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xorc.pam$clustering)
   
     1  2  3  4
  1 20  3  1  1
  2  0  2  2 21
  3  0  1 22  2
  4  0 22  0  3

# 正答率85%
# ※これまたラベルの順番がばらけているので要注意

> plot(xors[,-3],pch=19,cex=3,col=xors.pam$clustering,xlim=c(-3,3),ylim=c(-3,3))
> plot(xorc[,-3],pch=19,cex=3,col=xorc.pam$clustering,xlim=c(-3,3),ylim=c(-3,3))

f:id:TJO:20140228144659p:plain

K-medoidsのシンプルパターンと、

f:id:TJO:20140228144713p:plain

複雑パターンの結果。パフォーマンスで言うと1%だけ良くなってますが(笑)、このデータではあまり変わりがないと見る方が妥当でしょう。


むしろ、このK-meansとK-medoidsの結果で注目したいのが「実データよりもさらに真の分布に近いクラスタリング結果になっている」点かな、と。どちらも特に複雑パターンでは、実データの方がぐちゃぐちゃになっているのに対してクラスタリング結果はきちんと真の分布を反映して4象限にほぼ均等に分かれています。


そういう意味で言うと、K-means / K-medoidsは「より自然な真の分布形状に近付こうとする」のかな?と僕には思われるのですが、実際問題どんなものなんでしょうか? 教えて、偉い人!


最後に


最近になって色々また教師なし学習手法について知る機会があったので、その辺もこのシリーズで取り上げていければ良いかなぁと考えてます。例えばLatent Dirichlet Allocationとか、Affinity Propagationとか。。。RパッケージがあったりそれこそStanでやったりもできるみたいですが、まずは勉強せねば。

*1:だってPythonコード1文字たりとも書いてないし(笑)

*2:第9章「混合モデルとEM」の中

*3:厳密には深さ方向を逆に上に上ってくる

*4:もちろん物には限度があるので要注意

*5:ヒューリスティクスですな

*6:ただし僕自身の感覚でいうと例えばシンプレックス法みたいに極端に初期値に依存するということはない気がする

*7:説明がめんどくさいので詳細ははじパタp.156読んでくださいなー

*8:二値データなんかだと実はダメだったりする。Jaccard距離を使えとか色々な考え方があり。。。