混合ガウス分布とEMアルゴリズムを解説して実装してみる

混合ガウス分布とは

EMアルゴリズムとは

更新式の導出

実装

実際に混合ガウス分布に従うデータを生成して、得られたデータ列をEMアルゴリズムによってクラスタリングしてみます。

あ、EMアルゴリズムでクラスタリングするっていうのは、EMアルゴリズムによって推定した複数の分布の中からそのデータが一番出やすい分布に分類するという感じです。

$\gamma_i^j(t)$は$j$ステップ目においてデータ$x_t$が$i$番目の確率分布から現れたデータである確率を表しているので、この値が最大になるインデックス$i$に分類してしまえば良さそうです。

Python3で実装しました。結果の後ろにソースコードも載せてあるので参考にしてみてください。

前回のk-means法の時と似たようなデータをクラスタリングします。

これを4つのグループにクラスタリングします。うまくいくといいな〜〜

結果

イテレーションを5回行うごとにその時点でのパラメータを用いてクラスタリングしたものを載せていきます。

黒丸は推定した各グループの中心を表しています。等高線は各分布のガウス分布を表しています。2次元の確率分布も等高線を使えばわかりやすいですね。

 5回終了時。割と初期化した平均点が正解に近いところだったため、最初から比較的うまくクラスタリングされています。

初期化をランダムにしている部分があるので、同じコードを試しても結果は異なります。今回は運が良いですね。

 10回終了時。少し落ち着いてきました。

 15回終了時。だいたいもう終わってます。

20回で収束しました！綺麗にクラスタリングされています。わーい

対数尤度関数の値はこんな感じ。上昇して行って最後には増加が止まるので収束していることがわかります。

iter:1 likelihood:-3842.720847712288
iter:2 likelihood:-3649.8295196239096
iter:3 likelihood:-3483.322755557905
iter:4 likelihood:-3345.8369806950527
iter:5 likelihood:-3271.919667255609
iter:6 likelihood:-3242.2223846206166
iter:7 likelihood:-3228.6747589443034
iter:8 likelihood:-3221.5026106295704
iter:9 likelihood:-3217.1428785467497
iter:10 likelihood:-3213.7649598658277
iter:11 likelihood:-3209.998035905365
iter:12 likelihood:-3203.6673419243098
iter:13 likelihood:-3191.6960964562663
iter:14 likelihood:-3176.4622164239267
iter:15 likelihood:-3169.8894090988674
iter:16 likelihood:-3169.68998823923
iter:17 likelihood:-3169.6892371688186
iter:18 likelihood:-3169.6891464002156
iter:19 likelihood:-3169.6891259683107
iter:20 likelihood:-3169.6891213337417

混合係数の推定値はこんな感じになりました。

ピンク	赤	緑	青
0.09997503	0.1997618	0.29994096	0.40032222

生成元の混合係数は次のように設定したので、ほぼ正しく推定できてることがわかります。

ソースコード

import numpy as np
from genarateData import genarateData, plot
import math
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt

Colors = ["red", "green", "m", "blue"]
N = 1000
iter_num = 100
Means = [[2, 2], [-2, -2], [2, -2], [-2, 2]]
Covs = [[[0.5, 0], [0, 0.25]], [[0.3, 0], [0, 0.6]],
        [[0.5, 0], [0, 0.25]], [[0.3, 0], [0, 0.6]]]
Pi = [0.1, 0.2, 0.3, 0.4]
K = len(Pi)
d = len(Means[0])


def gaussian(x, y, mu, sigma):
    x = np.matrix([x, y])
    mu = np.matrix(mu)
    sig = np.matrix(sigma)
    assert sig.shape == (d, d), sig.shape
    if np.linalg.det(sig) * (2 * np.pi)**sig.ndim < 0:
        print("error!")
    a = np.sqrt(np.linalg.det(sig) * (2 * np.pi)**sig.ndim)
    b = np.linalg.det(-0.5 * (x - mu) * sig.I * (x - mu).T)
    return np.exp(b) / a


def likelihood(X, Y, pi, mu, sigma):
    L = sum([np.log(sum([pi[i] * gaussian(X[t], Y[t], mu[i], sigma[i]) for i in range(K)]))
             for t in range(N)])
    return L


def initialize(X, Y, classes):
    pi = [1 / classes] * classes
    mu = [np.random.uniform(min(X), max(X), d) for _ in range(classes)]
    sigma = [np.eye(d) for _ in range(classes)]
    return pi, mu, sigma


def expectation(X, Y, pi, mu, sigma):
    r = []
    for t in range(len(X)):
        denom = sum([pi[l] * gaussian(X[t], Y[t], mu[l], sigma[l])
                     for l in range(K)])
        numer = np.array([pi[i] * gaussian(X[t], Y[t], mu[i], sigma[i])
                          for i in range(K)])
        r.append(numer / denom)
    assert np.array(r).shape == (N, K)
    return r


def maximization(X, Y, r):
    N_k = np.sum(r, axis=0)
    assert len(N_k) == K
    pi = N_k / N
    assert len(pi) == K
    mu = []
    sigma = []
    for i in range(K):
        mu_i = np.zeros(d)
        for t in range(N):
            mu_i += r[t][i] * np.array([X[t], Y[t]])
        mu_i /= N_k[i]
        mu.append(mu_i)
        sigma_i = np.zeros((d, d))
        for t in range(N):
            temp = np.array([X[t], Y[t]]) - mu[i]
            sigma_i += r[t][i] * \
                np.matrix(temp).reshape(d, 1) * np.matrix(temp).reshape(1, d)
        sigma_i /= N_k[i]
        sigma.append(sigma_i)
    assert np.array(mu).shape == (K, d), mu.shape
    assert np.array(sigma).shape == (K, d, d), sigma.shape
    return pi, mu, sigma


def show_param(pi, mu, sigma):
    print("pi:{}".format(pi))
    print("mu:")
    for i in range(K):
        print(mu[i])
    print("sigma:")
    for i in range(K):
        print(sigma[i])


def EM_plot(X, Y, r, iter, mu, sigma):
    plt.title = "iter{}".format(iter)
    # クラスタリング
    for t in range(N):
        i = np.argmax([r[t][i] for i in range(K)])
        plt.scatter(X[t], Y[t],
                    c=Colors[i], marker="o", s=5)
    # 平均点と等高線の表示
    xlist = np.linspace(min(X), max(X), 100)
    ylist = np.linspace(min(Y), max(Y), 100)
    x, y = np.meshgrid(xlist, ylist)
    for i in range(K):
        plt.scatter(mu[i][0], mu[i][1], c="k", marker="o")
        z = bivariate_normal(x, y, np.sqrt(sigma[i][0][0]), np.sqrt(
            sigma[i][1][1]), mu[i][0], mu[i][1], sigma[i][0][1])
        cs = plt.contour(x, y, z, 3, colors='k', linewidths=1)

    plt.show()


def EM(X, Y, classes, iter_num):
    pi, mu, sigma = initialize(X, Y, classes)
    for j in range(iter_num):
        prev = likelihood(X, Y, pi, mu, sigma)
        r = expectation(X, Y, pi, mu, sigma)
        pi, mu, sigma = maximization(X, Y, r)
        now = likelihood(X, Y, pi, mu, sigma)
        print("iter:{} likelihood:{}".format(j + 1, now))
        if now - prev < 0.00001:
            print("converged!")
            show_param(pi, mu, sigma)
            EM_plot(X, Y, r, j + 1, mu, sigma)
            return pi, mu, sigma
        if (j + 1) % 5 == 0:
            EM_plot(X, Y, r,  j + 1, mu, sigma)
    print("not converged")
    return pi, mu, sigma


if __name__ == "__main__":
    X, Y = genarateData(Means, Covs, Pi, N)
    plot(X, Y, Pi)
    EM(X, Y, K, iter_num)

まとめ