Numpyを使って混合ガウス分布のEMアルゴリズムを実装

Pythonの練習ということで，Numpyを使って混合ガウス分布のEMアルゴリズムによる最尤推定を実装してみた．そもそもPythonを書いた経験があまり無いうえに，全く知らないNumpyを使って行列演算や確率計算をしようということで，手探りでかなり苦戦してしまったが，何とか形にはなったと思う．ということで，次の勉強に活かすためにもここでコードを振り返ってみる．

注意：以下のコードはテストデータでしか確かめてないので多分どこかバグってる．あと確率値に対数を取ってないので，値が限りなく小さくなってゼロ除算になることがある．

	import math
	import numpy as np

	def mv_norm(x, mu, sigma):
	norm1 = 1 / (math.pow(2 * math.pi, len(x)/2.0) * math.pow(np.linalg.det(sigma), 1.0/2.0))
	x_mu = np.matrix(x-mu)
	norm2 = np.exp(-0.5 * x_mu * sigma.I * x_mu.T)
	return float(norm1 * norm2)

	# test data
	data = np.array([[-1,-1],[-1,0],[0,1],[1,1],[1,2]])

	K = 2
	N = len(data)

	mu = [np.array([0, 0]), np.array([1, 0])]
	sigma = [np.eye(2), np.eye(2)]
	pi_k = [0.5, 0.5]

	L = []
	mu_iter = []
	sigma_iter = []
	pi_k_iter = []
	diff = 1

	while diff > 0.1 :
	# E-step
	likelihood = np.zeros((N,K))
	gamma_nk = np.zeros((N,K))
	for k in range(K):
	likelihood[:,k] = [mv_norm(d, mu[k], np.array(sigma[k]))*pi_k[k] for d in data]

	for n in range(N):
	gamma_nk[n,:] = likelihood[n,:] / sum(likelihood[n,:])

	# M-step
	N_k = np.array([sum(gamma_nk[:,k]) for k in range(K)])
	pi_k = N_k/sum(N_k)
	mu = np.dot(gamma_nk.T, data)/N_k
	for k in range(K):
	sig = 0
	for n in range(N):
	x_mu = data[n,:] - mu[:,k]
	sig += gamma_nk[n,k] * np.outer(x_mu, x_mu.T)
	sigma[k] = np.array(sig/N_k[k])

	# iteration
	mu_iter.append(mu)
	sigma_iter.append(sigma)
	pi_k_iter.append(pi_k)

	l = sum(map(np.log,[sum(likelihood[n,:]) for n in range(N)]))/N
	print l
	if L:
	diff = math.fabs(L[-1] - l)
	L.append(l)
	else:
	L.append(l)

view raw gmm_em.py hosted with ❤ by GitHub

実装

とっかかり

まずscipyあたりで多変量正規分布が無いか探す
- 普通の正規分布scipy.stats.normはあるが多変量は無いっぽい
- 自分で関数を書く
逆行列とか転置とか書き方が複数あって混乱する
- 今回はとりあえず動いたらいいので適当に使う
- np.matrixとnp.arrayが混ざってるの良くないっぽい

EMアルゴリズムの作り始め

Numpyの行列の表し方を一通り見る
- Rっぽい
- 0 originか1 originかで混乱する
更新するパラメータは基本的に入れ物となる変数を作っておいてそれに代入していく
- 最初はappendでリストに突っ込もうかと考えたけど，ndarrayとかになる場合もあるので面倒でやめた
Numpyの簡潔な書き方がわからない
- Rのapplyみたいに列指向で計算する方法がわからない
- 最初はリスト内包使ったりmap,zipつかったりしたけど，行列計算ぽく書けばそれなりに動くことに気付く
- それでもだめならfor文で直感的に計算する

終わりかけ

共分散行列の更新で意図した結果にならずにだいぶ悩む
- とりあえず人のを真似て書いてみる
対数尤度の収束規準をどうするかで悩む
- 今回は適当に更新したときの対数尤度の差を見るようにした
- 収束度合いの検出を改善するか，単調増加する対数尤度をプロットして適当なところで切るみたいなことをしたほうがいいかも

一通り書き終えて

まずはNumpyのスタイルを一通り身につけないといけない
- Numpyの流儀を学ぶ
- パッケージや方法が混在しているので，書き方を統一することが必要
実装の途中で数式の添字とかで混乱した
- EMアルゴリズムの資料としてPRMLともう一つ別のスライドを参考にしたのが原因
- 実装前にしっかりと頭のなかを整理しとくことが必要
全体的に対数尤度で計算する
- いわゆるlogsumexpでアンダフローに対処する

参考

Numpy/Scipy資料まとめ

英語

日本語

書籍