マハラノビス距離を用いた非階層的クラスタリング

こちらのサイトやこちらのサイトを参考に，豊田・池原（2011）によるk-means法改良版を実行するスクリプトを作成してみました。

この方法では，最初にk-means法を用いて基本的なクラスタリングしたうえで，各データポイントとクラスタ間のマハラノビス距離を求め，より最適な状態にクラスタリングし直します。

計算方法はシンプルなのにパワフル・・・と言いたかったのですが，私の理解不足あるいは技量不足のせいか，思ったほど使いやすいものにはなりませんでした。

しかし，せっかく作ったので，記録として一応残しておきます。以下では，まずスクリプトについて簡単に説明したうえで，（このスクリプトの）問題点について述べておきます。なお，実行にはstatsパッケージとmvtnormパッケージが必要です。

作成したスクリプト

mkmeans <- function (df, centers, iter.max =10,... ){
	require(stats)
	require(mvtnorm)

	# 繰り返し数上限
	iter.max <- iter.max

	# クラスタの数
	n_cl <- centers

	# 尤度行列
	L_mat<-numeric()

	# 尤度の計算
	calc_likelihood <- function(df, mu, sig){
		n <- length(mu)
		apply(df, 1, 
			function(x) dmvnorm(x, mean = mu, sigma =sig))
	}

    # クラスタごとの尤度
    dmnorm.cl <- function (df, cluster){
		# 平均ベクトル　μ_k
    	df.center <- sapply(split(df, cluster), function(x) sapply(x, mean))
		
		# 分散共分散行列 Σ_k
		df.cvmat <- lapply(split(df, cluster), function(x) var(x))

		# 尤度算出
		L_i <- lapply(1:n_cl, 
			function(y) calc_likelihood(df, df.center[,y], df.cvmat[[y]]))

		# リストをマトリクスに変換
		matrix(unlist(L_i),ncol=n_cl)

   }

	# マハラノビス距離によるクラスタリング
	mahalanobis_clustering <- function(x, cluster){
			# 2. 各個体の尤度の算出
			L_mat <<- dmnorm.cl(x, cluster)

			# 3. 尤度が最大のクラスタに再割り当て
			apply(L_mat, 1, which.max)
	}
	
###################################################
#### コメント中のナンバリングは，豊田・池原（2011）の
#### 記述と対応させたもの
###################################################

	# 1. 初期クラスタリング
	kcl <- kmeans(df, centers=n_cl,iter.max=iter.max, ...)
	iter<-0
	
	# 2-3. 再クラスタリング
	cl.old <- kcl$cluster
	cl.new <- mahalanobis_clustering(df, cl.old)
	while (!(all(cl.old == cl.new))){
			iter <- iter +1
			cl.old <- cl.new
			cl.new <- mahalanobis_clustering(df, cl.old)
			### 途中でクラスタが結合してしまった場合，初期化し直す
			if (length(table(cl.new)) != n_cl) {
				cl.new <- kmeans(df, centers=n_cl,iter.max=iter.max, ...)$cluster
			}
			if (iter > iter.max){
				cat('Maximum number of iterations exceeded\n')
				break
			}
		}
	
	# 4. 個体数と所属割合の算出
	w_k <- table(cl.new) / length(cl.new)

	# 5. 各個体の尤度算出　※ この辺の計算が怪しい？
	L1 <- sapply(1:n_cl, function(x) w_k[x] * L_mat[,x])
	L2 <- apply(L1,1,sum)
	L  <- prod(L2)

	# 6. AICとBICの算出
	# 自由母数
	K <- length(table(cl.new))
	p <- ncol(df)
	l <- p * (p +3	) * K / 2

	AIC <-   -2 * log(L) + 2 * l 
	BIC <-   -2 * log(L) + l * log (nrow(df)) 
    structure(list( "cluster"=cl.new,
                    "AIC"=AIC,
                    "BIC"=BIC))    
}

スクリプトの説明

このスクリプトでは，最初のクラスタリングにkmeans()を使用します。なお，kmeans()ではクラスタリングのアルゴリズムをHartigan-Wong，Lloyd，Forgy，MacQueenから選択することができますが，このスクリプトでもオプションはそのままkmeans()に渡すようにしているので，kmeans()のアルゴリズムを指定すれば，その方法でクラスタリングされます。なお，豊田・池原（2011）では，初期クラスタリングの方法としてMacQueenが使用されています。

最初のクラスタリングが終わった後は，各レコードについてそれぞれのクラスタに対する尤度を求め，尤度が最大のクラスタに再分類する手続きを，結果が収束するか，反復数の最大になるまで繰り返します。

なお，論文中にはとくに書かれていないのですが，自分の手持ちのデータで色々試していたところ，データによっては途中で複数のクラスタが結合されてしまい，計算がストップしてしまうケースがありました。そこで，スクリプト中では，そのような場合には再度k-means法で初期化して再計算するようにしています。

結果が収束した場合には，クラスタリングの結果と，各レコードの尤度をもとにしたAICとBICを出力しておしまいです。この方法ではAICやBICを算出することができるので，それらを参考にしながら最適なクラスター数を決定できるのが強みとされています。

（このスクリプトの）問題点

なお，分類結果自体は豊田・池原（2011）で報告されているものとほぼ同じなのですが，このスクリプトで算出されるAICとBICの値は，論文のものとは随分桁の違うものになっています。対数尤度の算出方法を根本的に間違えている可能性もありますので，その様な場合にはどうかご指摘ください。

また，これはスクリプトの問題というよりはこの手法そのものの問題なのだと思いますが，初期解として用いるk-means法の結果によっては，変なところで結果が収束してしまう場合があります。irisデータや紙幣データでクラスタの数をいじりながら何度か実行していたところ，そこそこ高い割合で変な結果が出てきました。irisデータの場合ですが，うまく結果が収束すると次の表の上の様に分類されるのですが，表の下の様な結果が出てくることが，感覚的には10回に1回ぐらいの割合であります。

正しい分類結果

	1	2	3
setosa	0	50	0
versicolor	48	0	2
virginica	1	0	49

おかしな分類結果

	1	2	3
setosa	0	33	17
versicolor	46	0	4
virginica	50	0	0