教師なし学習シリーズもいよいよ佳境に入ってきましたねー、と言いつつ前回記事から既に2ヶ月半ぐらい経ってますが。。。ここからは主に混合モデルを取り上げていく予定です。今回もはじパタpp.165-174をベースにやっていきます。
- 作者: 平井有三
- 出版社/メーカー: 森北出版
- 発売日: 2012/07/31
- メディア: 単行本（ソフトカバー）
- 購入: 1人 クリック: 7回
- この商品を含むブログ (4件) を見る
もちろん細かいところはPRMLなどでチェックすると良いかと思います。PRMLなら下巻の第9章が丸々そのまま混合モデルとEMアルゴリズムの解説に充てられているので、はじパタで説明が足りないなと思ったところについては参照してみてください。
まずRでどんなものか見てみる
いつも通りですが、GitHubに置いてある多変量データで試してみましょう。これまたいつも通りdとかいう名前でインポートしておきます。
混合モデルとEMアルゴリズムに基づくクラスタリングは、Rでは{mclust}パッケージで実装されています。インストールしてから以下のように実行してみましょう。
> require("mclust") > d <- read.table("~/experiments/R_work/conflict_sample.txt", header=T, quote="\"") > plot(mclustBIC(d[,-8])) > mhc<-hc(modelName="EEV",data=d[,-8]) Error in hc(modelName = "EEV", data = d[, -8]) : invalid model name for hierarchical clustering # "EEV"はhc()関数に実装されていない > mhc<-hc(modelName="EII",data=d[,-8]) > d.cl<-hclass(mhc,2) > table(d[,8],d.cl) d.cl 1 2 No 1419 81 Yes 310 1190 # 試しに真のラベルと比較してみたら正答率87%弱という結果に # 意外と悪いな。。。
ちなみにこれはplot(mclustBIC(d[,-8])して得られるプロットです。mclustBIC()関数は混合正規分布モデルの正規分布の「形状」に対して、例えば球形or楕円球形、等分散or可変分散、等体積or異なる体積、同じ形状＆向きor異なる形状＆向き（楕円球形の場合）といったモデルのそれぞれのBICを算出してくれます。そこで、目指すクラスタ数において最もBICの高いものを選べば良い、というわけです。
ただし、hc()関数パッケージに実装されているモデルには限りがあるので、?hcとして実装されているモデルを確認しておいた方が良いでしょう。上の例でも実装されているモデルの中でBICが最も高かった"EII"（球形＆等体積）を選んでいます。
混合モデル＆EMアルゴリズムとは何ぞや
ぶっちゃけ、今回の話はものすごく難しいのでかなり端折ります*1。数式が全部並ぶと辛い人も多い*2と思うので、要点だけ以下に書いておきます。詳しい式ははじパタpp.165-174およびPRML下巻第9章を熟読してください。
ということで、式を端折った要点は以下の通りです。
- まず混合正規分布関数とその線形和を立式する（はじパタ10.4.1節）
- 隠れ変数を置き、これと実際に観測される変数との同時分布から式変形して最終的に隠れ変数の事後確率を立式する（10.4.2節）
- 完全データ（観測データと隠れ変数とを混ぜたもの）の対数尤度を立式する（10.4.3節）
- 隠れ変数に関する期待値は、式変形の都合上隠れ変数の事後確率と等しくなる（Eステップ）。そこで対数尤度関数の隠れ変数の箇所を、隠れ変数の事後確率に置き換えたQ関数（Mステップ）を立てる（10.4.4節）
あー、式書かないと簡単ですわ（笑）。ポイントとしては、1番目に立式されたモデルのパラメータは基本的に隠れ変数の事後確率から算出することが可能な一方で、隠れ変数の事後確率自体もそのパラメータから求められ、他方で完全データの対数尤度はこれらとは別に算出されるという点が重要です。
そこで、1番目の式のパラメータに適当な初期値を与えて適当に隠れ変数を算出し（Eステップ）、そこからパラメータを逆算し（Mステップ）、対数尤度関数に代入し、対数尤度が収束するまでこれを繰り返す、というやり方で逐次最大化することで最尤推定量を求めることができるわけです。これをEMアルゴリズムと呼びます。
なお、EMアルゴリズム自体はもっと一般的なもので、Kullback-Leibler divergenceとか皆さん大好きそうな要素*3の散りばめられた大変楽しいものなので、興味のある人は是非はじパタ10.4.7節やPRML下巻9.4節などを参照してみてください。
そうそう、1番目に立式されたモデルのところに混合正規分布関数が出てきますが、これを変えることで冒頭にコメントした球形or楕円球形みたいなモデルの仮定を与えることができます。モデルの妥当性はBICで求めれば良いというわけです。
真の分布の分かっているデータをクラスタリングしてみる
そんなわけで、これまたいつも通り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つ配置されたデータです。以上の点を踏まえて、xors / xorc双方のデータに対して{mclust}パッケージを用いてEMアルゴリズムによる混合正規分布モデルクラスタリングを行ってみましょう。
> plot(mclustBIC(xors[,-3])) > xors.mhc<-hc(modelName="EII",data=xors[,-3]) > xors.cl<-hclass(xors.mhc,4) > table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xors.cl) xors.cl 1 2 3 4 1 25 0 0 0 2 0 24 1 0 3 0 0 0 25 4 1 0 24 0 # 正答率98%となかなか > plot(xors[,-3],pch=19,col=c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),cex=3,main="True") > plot(xors[,-3],pch=19,col=xors.cl,cex=3,main="Predicted") > plot(mclustBIC(xorc[,-3])) > xorc.mhc<-hc(modelName="EII",data=xorc[,-3]) > xorc.cl<-hclass(xorc.mhc,4) > table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xorc.cl) xorc.cl 1 2 3 4 1 14 5 5 1 2 0 12 12 1 3 0 1 10 14 4 0 19 6 0 # 正答率59%。。。だいぶ低いっすね > plot(xorc[,-3],pch=19,col=c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),cex=3,main="True") > plot(xorc[,-3],pch=19,col=xorc.cl,cex=3,main="Predicted")
シンプルパターンと、
複雑パターンの結果。単純に正解ラベルとの比較でいうと、シンプルパターンでは正答率はなかなかなんですが、複雑パターンはちょっとどうしようもないですね（汗）。そこで、ちょっと変えてみました。
> xorc.mhc2<-hcVII(xorc[,-3],minclus=4) # "VII"（球形＆異なる体積）にモデル変更 > xorc.cl2<-hclass(xorc.mhc2,4) > table(c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),xorc.cl2) xorc.cl2 1 2 3 4 1 14 2 4 5 2 0 19 2 4 3 0 2 1 22 4 1 7 17 0 # 正答率72%。。。ちょっと改善したかも > par(mfrow=c(1,2)) > plot(xorc[,-3],pch=19,col=c(rep(1,25),rep(2,25),rep(3,25),rep(4,25)),cex=3,main="True") > plot(xorc[,-3],pch=19,col=xorc.cl2,cex=3,main="Predicted")
見た感じでは前よりもマシになり、なおかつ真の分布に少し近づいた気がします。mclustBIC()で描いたプロットから見ると、G = 3のところではVIIモデルのBICが高かったので試しにVIIに変えてみたというだけなんですが、随分印象が違って見えますね。
ということで、混合モデルの場合は分布の選び方によってかなりクラスタリング性能に差が出るのかな？とも思いました。もっと高次元化した場合はどうなるのか分かりませんが、注意しておきたいところです。
最後に
次回は混合ディリクレ過程のつもりですが、多分「何ぞや」のところもまともに書けないと思われます（泣）。そして意外と良い資料が少なくて、書籍となると和書ではほぼ皆無、洋書でも稀少ということらしく。。。どなたか良い資料ご存じないですか？
