EMアルゴリズム答え合わせ

自分の実装した「Numpyで混合ガウス分布のEMアルゴリズムを実装した」のコードを中谷さんの「EM アルゴリズム実装(勉強用) - Mi manca qualche giovedi`?」と照らしあわせて答え合わせしてみる。

まず、EMアルゴリズムってなんなのかって話を簡潔に。観測できている確率変数Xの他に観測できてない確率変数Zもある状況。それを表現するために自分で何か確率モデルを仮定する。その確率モデルのパラメータθを、観測されたデータXから考えてもっとも良さげなものを選びたい。コレが目的。数式で言えばp(X|θ)を最大にするθを求めたい。

但し残念なことにp(X | θ)の式は簡単には最大化できない(できるならEMアルゴリズム使わなくていい)とする。そして幸運なことに完全データの対数尤度 ln p(X, Z | θ)の最大化は簡単だと仮定する。PRMLの下巻p.166あたりの議論と持橋さんの「自然言語処理のための 変分ベイズ法」を参考に、僕が分かる範囲で簡単に書くと:

• p(X | θ)は簡単には最大化できないが、Jensenの不等式によってそれがある式L(q(Z), θ)より大きいことが言える
• θを止めてq(Z)についてLを最大化することは、LからZに関係のないln p(X|θ)をくくりだしてみると残りがKL(q||p)なので、qをpにすることで可能(E-step)
• q(Z)を止めてθについてLを最大化することは、Lからθに無関係な分母をくくりだすと…えーとんーとQ関数ってどういう意味のものなんだろうか。完全データの対数尤度にp(Z | X, θ^)を掛けてzで和を取ったものなんだが。p(Z | X, θ^)がθに無関係な係数なので…うーん？
• メモ：EMアルゴリズムのMステップで、Q関数のargmax_thetaが、どうして完全データの対数尤度を微分して0ってことになるのかがわからない。

EM アルゴリズムについて説明するのが目的ではないからこれくらいにしておいて(コピペ)ソースコードの比較に移ろう。

データの準備

cls1 = c_[randn(100), randn(100)].dot(array([[10, 0], [0, 1]])).dot(array([[1, -1], [1, 1]]))
cls2 = c_[randn(100), randn(100)].dot(array([[10, 0], [0, 1]])).dot(array([[1, 1], [-1, 1]]))
data = r_[cls1, cls2]

中谷さんはOld Faithfulのデータを使っているけど、僕はk-meansでうまくできないデータがEMでうまくいくことを確認したかったので自分でクロスした分布を作っている。

パラメータの初期化

# クラス数
K <- 2;

# 平均、共分散、混合率の初期値(正規乱数)
mu <- matrix(rnorm(K*ncol(xx)), nrow=K);
mix <- numeric(K)+1/K;
sig <- list();
for(k in 1:K) sig[[k]] <- diag(ncol(xx));

# inital parameter
K = 2
D = 2
mus = [array([1, 1]), array([-1, -1])]
sigmas = [eye(2), eye(2)]
pis = [0.5, 0.5]
assert len(mus) == K and all(len(mu) == D for mu in mus)
assert len(sigmas) == K and all(s.shape == (D, D) for s in sigmas)
assert len(pis) == K

中谷さんのはちゃんとパラメータを変えられるように作ってあるな。僕のは、ハードコーディングだった。この機会に修正しておこう。実はk-meansを実装するときにそっちでは2個以上のクラスターに対応したんだった。

# inital parameter
K = 2
D = 2
_degrees = array(range(K)) * 2 * pi / K
mus = array([array([sin(th), cos(th)]) for th in _degrees])
sigmas = [eye(D) for k in range(K)]
pis = [1.0 / K] * K
assert len(mus) == K and all(len(mu) == D for mu in mus)
assert len(sigmas) == K and all(s.shape == (D, D) for s in sigmas)
assert len(pis) == K

ガウス分布の密度関数

# 多次元正規分布密度関数(パッケージ使えって？)
dmnorm <- function(x,mu,sig) {
	D <- length(mu);
	1/((2 * pi)^D * sqrt(det(sig))) * exp(- t(x-mu) %*% solve(sig) %*% (x-mu) / 2)[1];
}

def dens_gauss(x, mu, sigma):
    """
    calculate gauss distribution's density
    """
    Z = sqrt(2 * pi) ** x.size * sqrt(norm(sigma))
    v = x - mu
    return exp(-0.5 * v.dot(inv(sigma)).dot(v)) / Z

まあ同じだよね、と言おうとしたけど、中谷さんの方は2 * piに平方根がかかってないように見えるなぁ。

そして僕の方はnormじゃなくてdetを使うべきだったね。

In [1080]: norm(eye(2))
Out[1080]: 1.4142135623730951

In [1081]: det(eye(2))
Out[1081]: 1.0

Eステップ

# EM アルゴリズムの E ステップ
Estep <- function(xx, mu, sig, mix) {
	K <- nrow(mu);
	t(apply(xx, 1, function(x){
		numer <- sapply(1:K, function(k) {
			mix[k] * dmnorm(x, mu[k,], sig[[k]])
		});
		numer / sum(numer);
	}))
}

def resp(x, pis, mus, sigmas):
    """
    given: vector x, list of mu, list of sigma
    return: array of responsibility
    eq.9.13
    """
    v = array([p * dens_gauss(x, mu, sigma) 
               for p, mu, sigma in zip(pis, mus, sigmas)])
    s = v.sum()
    return v / s

(...snip...)

# E-step
# \gamma(z_nk)
# responsibility
resp_for_each_data = array([resp(x, mus, sigmas) for x in data])
assert len(resp_for_each_data) == N

# eq 9.18, Nk
num_for_each_cluster = array([col.sum() for col in resp_for_each_data.transpose()]) 
assert len(num_for_each_cluster) == K

ほぼ同じ。中谷さんはNkをMステップの中で求めてるのね。
そこだけ切り出してくるとこんな感じ。column方向に足す関数があるのか。Numpyにも探せば有りそうだなぁ。追記、sum(1)でok

N_k <- colSums(gamma_nk);

Mステップ

π

new_mix <- N_k / N;

# eq.9.22 update \pi_k
def update_pis(num_for_each_cluster):
    return num_for_each_cluster / N

まあそれ以外の書きようがないよね。

μ

new_mu <- (t(gamma_nk) %*% xx) / N_k;

# update parameters, M-step
# eq 9.17, update \mu_k
def update_mus(resp_for_each_data, data, num_for_each_cluster):
    return [
        resp_for_each_data[:, k].dot(data)
        / num_for_each_cluster[k]
        for k in range(K)]

あ、そうか、同じ長さ同士の配列で割り算をすればエレメントごとの割り算になるのか。なるほど。書きなおそう。

resp_for_each_data.T.dot(data) / num_for_each_cluster

あれ、結果が異なる。

(resp_for_each_data.T.dot(data).T / num_for_each_cluster).T

これなら一致するけどなんでだ？？

In [1078]: array([[1, 1], [1, 1]]) / array([1.0, 2.0])
Out[1078]: 
array([[ 1. ,  0.5],
       [ 1. ,  0.5]])

あ、そうなのか。「1次元目の長さが同じであればエレメントワイズになる」という解釈が間違いで、shapeが違うならブロードキャストになるのね。だから転置してから割って、結果を転置して戻せば機体と同じものになる、と。でもそんなコード嫌なので、んー、どうするのがエレガントかなぁ。良い書き方がありそうだけどとりあえず今は素直にzipで書いておくか。

array([v / s for v, s in zip(resp_for_each_data.T.dot(data), num_for_each_cluster)])

Σ

	new_sig <- list();
	for(k in 1:K) {
		sig <- matrix(numeric(D^2), D);
		for(n in 1:N) {
			x <- xx[n,] - new_mu[k,];
			sig <- sig + gamma_nk[n, k] * (x %*% t(x));
		}
		new_sig[[k]] <- sig / N_k[k]
	}

# eq. 9.19, update \Sigma_k
def update_sigmas(resp_for_each_data, data, mus, num_for_each_cluster):
    return [
        array([resp_for_each_data[n, k] 
               * outer(data[n] - mus[k], data[n] - mus[k])
               for n in range(N)]).sum(0)
        / num_for_each_cluster[k] for k in range(K)]

僕がsum(0)、つまり0軸目方向のsumで合算しているところを中谷さんはforで足しあわせている。そのかわり僕は2回書いちゃってるdata[n] - mus[k]を中谷さんはxに代入していて式が読みやすいかも。

まとめ

同じ内容を見ないで実装してからあとで照らし合わせると結構勉強になる。あとRとNumpyは似ている。