Hiroshi TAKEMOTO (take@pwv.co.jp)

第9章-混合ガウス分布のEMアルゴリズムをSageで試す

混合ガウス分布の最尤推定問題をEMアルゴリズムで解いてみます。

ここで紹介するのは、PRMLの図9.8、

です。反復を繰り返して収束する様子がよく分かります。

テストデータ

テストデータには、Old Faithful間欠泉データをX-Y平面に-2.5から2.5の範囲に正規化したものを使用します。

{{{id=0| # 混合ガウス分布の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) /// }}} {{{id=2| # ベクトルのガウス分布を定義 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]; /// }}} {{{id=12| # 対数尤度 def ln_p(): sum = 0.0 for n in range(N): sum += ln(pi_N[n].sum()) return sum /// }}} {{{id=20| # (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)]) /// }}} {{{id=6| # 定数の初期化 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} /// }}} {{{id=18| # アクション関数 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 \} $$

{{{id=4| # 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に近い結果を得ることができました。

{{{id=21| # 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アルゴリズムもμの初期値に依存するのではないかと思います。

{{{id=22| # μの初期値に結果が大きく依存する # 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) /// }}} {{{id=29| /// }}}