本シリーズ記事のカテゴリからPythonが消えて久しい今日この頃ですが、皆様いかがお過ごしでしょうか*1。とかいう前口上はどうでも良くて、とっとと今回のお題に入りましょう。今回はクラスタリングのド定番、K-means（k平均）クラスタリングです。
K-meansはどの本にも載っている手法なので参考文献には何を挙げても良い気がしますが、とりあえずこれまたド定番ということで相変わらずはじパタを挙げておきます。
- 作者: 平井有三
- 出版社/メーカー: 森北出版
- 発売日: 2012/07/31
- メディア: 単行本（ソフトカバー）
- 購入: 1人 クリック: 7回
- この商品を含むブログ (4件) を見る
ちなみに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個のデータ
を、データ間の類似性（距離）を尺度に、あらかじめ定めたK個のクラスタに分類することである。
前回の階層的クラスタリングと比べてみると、
- データ間距離を尺度に決めている点で同じ
- クラスタ個数Kを事前に固定する点が異なる
と言ったところでしょうか。階層的クラスタリングは深さ方向（クラスタ個数が増える方向）にどんどん行ってしまう*3のでデータ容量が大きくなると計算量が爆発して簡単にコケてしまいますが、K-meansはクラスタ個数が固定されているので大容量データでも動く点がメリットと言って良いと思います*4。そのアルゴリズムは同じくはじパタp.155によればこんな感じです。
何となく雰囲気は分かると思うんですが、Wikipedia英語版にこのプロセスを分かりやすく図説にしたものがあるので貼っておきます。
1) k initial "means" (in this case k=3) are randomly generated within the data domain (shown in color).
2) k clusters are created by associating every observation with the nearest mean. The partitions here represent the Voronoi diagram generated by the means.
3) The centroid of each of the k clusters becomes the new mean.
4) Steps 2 and 3 are repeated until convergence has been reached.
こんな感じに、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つ配置されたデータです。プロットすると以下にようになるはず。
という点を踏まえて、実際に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))
K-meansのシンプルパターンと、
複雑パターンの結果。上の方の「正解」のプロットと見比べてもらえば分かるかと思いますが、大体問題なくクラスタリングできてる感じです。
ところで、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))
K-medoidsのシンプルパターンと、
複雑パターンの結果。パフォーマンスで言うと1%だけ良くなってますが（笑）、このデータではあまり変わりがないと見る方が妥当でしょう。
むしろ、このK-meansとK-medoidsの結果で注目したいのが「実データよりもさらに真の分布に近いクラスタリング結果になっている」点かな、と。どちらも特に複雑パターンでは、実データの方がぐちゃぐちゃになっているのに対してクラスタリング結果はきちんと真の分布を反映して4象限にほぼ均等に分かれています。
そういう意味で言うと、K-means / K-medoidsは「より自然な真の分布形状に近付こうとする」のかな？と僕には思われるのですが、実際問題どんなものなんでしょうか？ 教えて、偉い人！
最後に
最近になって色々また教師なし学習手法について知る機会があったので、その辺もこのシリーズで取り上げていければ良いかなぁと考えてます。例えばLatent Dirichlet Allocationとか、Affinity Propagationとか。。。RパッケージがあったりそれこそStanでやったりもできるみたいですが、まずは勉強せねば。
