Hatena::ブログ(Diary)

西尾泰和のはてなダイアリー

2011-11-09

numpyのndarrayで内積と外積を計算

内積(v^t v)や行列の積はndarrayのdotメソッドでできるが、外積(v v^t)はどうやってやるんだろう…それらしきメソッドがないなぁ…と思ったらouterって関数があった。めでたしめでたし。

In [833]: array([1, 2, 3])
Out[833]: array([1, 2, 3])

In [834]: v = _

In [835]: v.dot(v)
Out[835]: 14

In [836]: outer(v, v)
Out[836]: 
array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])

numpyで対角行列を作る

与えられた成分が対角成分に入っているような対角行列を作りたい場合、diagを使えばよい。

In [837]: diag([1, 2, 3])
Out[837]: 
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])

NumPy+Matplotlibで散布図の頂点に色をつける

散布図のマーカーに"+"を指定して「あれー、色を指定しても反映されないなー」と悩んでいたが、指定した色はマーカーの塗りに使われるのであった。"+"だと塗りがないから色が変わらないというオチ。edgecolorを指定すればいいんだろうけどマーカーを"o"に変えたらなかなかかっこ良かったのでこっちにする。

f:id:nishiohirokazu:20111109101621p:image

scatterの引数cに各頂点の色を指定すればよい。下のコードではresp_for_each_dataが合計すると1になる2要素のベクトルになっているので、その比率で赤と緑を混ぜ合わせてcolorsを作っている。

    red = array([1, 0, 0])
    green = array([0, 1, 0])
    colors = [red * r[0] + green * r[1] for r in resp_for_each_data]
    scatter(data[:, 0], data[:, 1], c=colors, alpha=0.5, marker="o")

numpyでk-means法を実装した

f:id:nishiohirokazu:20111109120733p:image

NumPyすごいな。学習部分は実質2行だ。

# E-step
nearest_cluster = array([argmin([norm(x - mu) for mu in mus]) for x in data])

# M-step
mus = [average(data[nearest_cluster == k], axis=0) for k in range(K)]


「パターン認識と機械学習」(PRML)読書会 #11 + R で K-means - Mi manca qualche giovedi`?

# K-means の1ステップ
(mu <- t(sapply(1:K,function(k)colMeans(xx[max.col(-sapply(1:K,function(i)apply(xx,1,function(x)sum((mu[i,]-x)^2))))==k,]))));

むむっ

mu = [average(xx[array([argmin([norm(x - m) for m in mu]) for x in xx]) == k], axis=0) for k in range(K)]

よし。(何)


クラスタ数を増やしたら枝が1本ずつ取られるかとおもいきやこんなことになった。

f:id:nishiohirokazu:20111109124003p:image

numpy+matplotlibで散布図の上にバツ印をつける

「バツ印をつける」と言った場合、多くの場合言語化されていない暗黙の要求仕様がある。「バツ印はグラフの軸に影響されて横長になったりしないでほしい」とか。

f:id:nishiohirokazu:20111109144009p:image

前回 「NumPy+Matplotlibで散布図の上に平均と分散を表示する」 ではCircleをPatchCollectionに入れて云々、とやったが、今回はそうではなくバツ印でクラスタの中心だけを表現したい。Circleを置いたのでは軸に影響されて横長になったりして期待と違う。

f:id:nishiohirokazu:20111111000335p:image

どうしたらいいかなぁ、としばらく考えてから、scatterでいいじゃないかと気づいた。散布図のデータ点は拡縮されていないからね。白で太めに書いてからその中に各クラスターの色で十字を書いている。とても手軽に実現できた。

    clf()
    colors = color_for_cluster[nearest_cluster]
    scatter(data[:, 0], data[:, 1], c=colors, alpha=0.5, marker="o")
    scatter(mus[:, 0], mus[:, 1], s=50, linewidths=3, edgecolor=[1, 1, 1], marker="+")
    scatter(mus[:, 0], mus[:, 1], s=50, edgecolor=color_for_cluster, marker="+")

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は似ている。