第9章-混合ガウス分布のEMアルゴリズムをSageで試す
混合ガウス分布の最尤推定問題をEMアルゴリズムで解いてみます。
ここで紹介するのは、PRMLの図9.8、
テストデータ
テストデータには、Old Faithful間欠泉データをX-Y平面に-2.5から2.5の範囲に正規化したものを使用します。
# 混合ガウス分布のEMアルゴリズム
from numpy import loadtxt
data = loadtxt(DATA+"faithful.txt")
# 1カラムはx, 3カラムにはyがセットされている
# これを-2.5, 2.5に正規化する
x = data[:,0]- 3.5
y = (data[:,1]-70)*5/60
# プロット
data_plt = list_plot(zip(x,y))
data_plt.show(aspect_ratio=1, figsize=(3), xmin=-2.5, xmax=2.5, ymin=-2.5, ymax=2.5)
|
# ベクトルのガウス分布を定義
def _gauss(v, mu, sigma):
d = len(v);
sigma_inv = sigma.inverse();
sigma_abs_sqrt = sigma.det().sqrt();
# sage 4.6.2ではtranspose()の代わりにcolumn()を使用する
val = -(v - mu) * sigma_inv * (v - mu).column()/2;
a = (2*pi)**(-d/2) * sigma_abs_sqrt**-0.5
return a * e**val[0];
|
|
# 対数尤度
def ln_p():
sum = 0.0
for n in range(N):
sum += ln(pi_N[n].sum())
return sum
|
|
# (x - u)(x - u)^Tを計算
def covMatrix(x, u):
d = x - u
return matrix(RDF, [[d[i]*d[j] for j in range(D)] for i in range(D)])
|
|
# 定数の初期化
D = 2
K = 2
N = len(x)
beta = 0.5
X = matrix(RDF, zip(x, y))
ident = identity_matrix(2)
# アクションのインデックス
action_idx = {0:0, 1:1, 5:5, 20:20}
|
|
# アクション関数
def action_f(mu_k, sigma_k, gamma):
x,y = var('x y')
pt_plt = Graphics()
for n in range(N):
pt_plt += point(X[n], rgbcolor=(gamma[n][1], 0, gamma[n][0]))
r_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu_k[1], sigma_k[1]), [x, -2.5, 2.5], [y, -2.5, 2.5], contours = 1, cmap=['red'], fill=False)
b_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu_k[0], sigma_k[0]), [x, -2.5, 2.5], [y, -2.5, 2.5], contours = 1, cmap=['blue'], fill=False)
(r_cnt_plt + b_cnt_plt + pt_plt).show(aspect_ratio=1, figsize=(3), xmin=-2.5, xmax=2.5, ymin=-2.5, ymax=2.5)
|
|
EMアルゴリズム
混合ガウス分布のEMアルゴリズムは、以下の通りです。
- 平均$\mu_k$、 分散$\Sigma_k$、混合係数$\pi_k$を初期化する
- Eステップ：式(9.23)を使って負担率を計算する $$ \gamma(z_{nk}) = \frac{\pi_k \mathcal{N}(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^N \pi_j \mathcal{N}(x_n|\mu_j, \Sigma_j)} $$
- Mステップ：上記の負担率を使って式(9.24), (9.25), (9.26), (9.27)のパラメータ値を再計算する $$ \mu_k^{new} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma(z_{nk})x_n $$ $$ \Sigma_k^{new} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma(z_{nk}) (x_n - \mu_k^{new}) (x_n - \mu_k^{new})^T $$ $$ \pi_k^{new} = \frac{N_k}{N} $$ $$ N_k = \sum_{n=1}^{N} \gamma(z_{nk}) $$
- 対数尤度：式(9.28)を計算し、対数尤度が収束するまでステップ２に戻る $$ \ln p(X|\mu, \Sigma, \pi) = \sum_{n=1}^{N} ln \left \{ \sum_{k=1}^{K} \pi_k \mathcal{N}(x_n|\mu_k, \Sigma_k) \right \} $$
# EMアルゴリズム
def _EM(pi_k, mu_k, sigma_k):
# 対数尤度の初期値
o_lnP = ln_p()
diff = 1
for l in range(21):
# Eステップ
pi_N = matrix(RDF, [[pi_k[k]*_gauss(X[n], mu_k[k], sigma_k[k]) for k in range(K)] for n in range(N)])
gamma = matrix(RDF, N, K)
for n in range(N):
gamma.set_row(n, pi_N[n]/(pi_N[n].sum()))
# 最初のEステップの状態
if l == 0:
action_f(mu_k, sigma_k, gamma)
# Mステップ
N_k = [gamma.column(k).sum() for k in range(K)]
for k in range(K):
mu_k[k] = 0
for n in range(N):
mu_k[k] += gamma[n][k]*X[n]
mu_k[k] /= N_k[k]
sigma_k[k] = 0
for n in range(N):
sigma_k[k] = sigma_k[k] + gamma[n][k]*covMatrix(X[n], mu_k[k])
sigma_k[k] /= N_k[k]
pi_k[k] = N_k[k]/N
# 新しい対数尤度
lnP = ln_p()
diff = abs(lnP - o_lnP)
o_lnP = lnP
# 図化
if action_idx.has_key(l):
action_f(mu_k, sigma_k, gamma)
|
|
計算結果
平均の初期値を$\mu_1 = (-1.5, 0.5), \mu_2 = (1.5, -0.5)$、 分散の初期値を0.5として、 計算した結果を以下に示します。
図9.8に近い結果を得ることができました。
# mu=(-1.5, 0.5), (1.5, -0.5)で計算
pi_k = [0.5, 0.5]
mu_k = [vector([-1.5, 0.5]), vector([1.5, -0.5])]
sigma_k = [matrix(RDF, D, D, beta*ident), matrix(RDF, D, D, beta*ident)]
pi_N = matrix(RDF, [[pi_k[k]*_gauss(X[n], mu_k[k], sigma_k[k]) for k in range(K)] for n in range(N)])
_EM(pi_k, mu_k, sigma_k)
    
|
μの初期値に依存
平均の初期値が図9.8と少しずれていたので、$\mu_1 = (-1.5, 1.0), \mu_2 = (1.5, -1.0)$ としたところ、途中まで上手く計算していたのですが、L=20回目がまったく異なる結果になってしまい ました。
混合密度ネットワークの時のように混合ガウス分布のEMアルゴリズムもμの初期値に依存するのでは ないかと思います。
# μの初期値に結果が大きく依存する
# mu=(-1.5, 1.0), (1.5, -1.0)で計算
pi_k = [0.5, 0.5]
mu_k = [vector([-1.5, 1.0]), vector([1.5, -1.0])]
sigma_k = [matrix(RDF, D, D, beta*ident), matrix(RDF, D, D, beta*ident)]
pi_N = matrix(RDF, [[pi_k[k]*_gauss(X[n], mu_k[k], sigma_k[k]) for k in range(K)] for n in range(N)])
_EM(pi_k, mu_k, sigma_k)
    
|
|
|