sage/PRML - 混合ガウス分布のEMアルゴリズム

FrontPage

2011/06/29からのアクセス回数 9195

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://sage.math.canterbury.ac.nz/home/pub/108/

また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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

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

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

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

テストデータ †

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

sageへの入力:

# 混合ガウス分布の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)

sageへの入力:

# ベクトルのガウス分布を定義
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).transpose()/2;
    a = (2*pi)**(-d/2) * sigma_abs_sqrt**-0.5
    return a * e**val[0];

sageへの入力:

# 対数尤度
def ln_p():
    sum = 0.0
    for n in range(N):
        sum += ln(pi_N[n].sum())
    return sum

sageへの入力:

# (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)])

sageへの入力:

# 定数の初期化
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}

sageへの入力:

# アクション関数
def action_f(mu_k, sigma_k, gamma):
    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アルゴリズムは、以下の通りです。

1. 平均\(\mu_k\)、 分散\(\Sigma_k\)、混合係数\(\pi_k\)を初期化する
2. 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)} $$
3. 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 \} $$

sageへの入力:

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

sageへの入力:

# 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)]
_EM(pi_k, mu_k, sigma_k)

μの初期値に依存 †

平均の初期値が図9.8と少しずれていたので、\(\mu_1 = (-1.5, 1.0), \mu_2 = (1.5, -1.0)\) としたところ、途中まで上手く計算していたのですが、L=20回目がまったく異なる結果になってしまい ました。

これはもしかすると、図9.7

で説明されている特異性の例かも知れない。

sageへの入力:

# μの初期値に結果が大きく依存する
# 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)]
_EM(pi_k, mu_k, sigma_k)

コメント †

皆様のご意見、ご希望をお待ちしております。

お名前:

(Input image string)