Hope is a Dream. Dream is a Hope.

非公開ブログは再開しました。

GMMとEMアルゴリズム | Python + Numpy

GMM混合ガウス分布MixtureGaussianModel

混合ガウス分布と多変量ガウス分布は違うものだよ。EMは実装が容易なので、手を動かすとすぐに理解できます。参考資料は自分で探すこと。

実装例

#! coding:utf-8
"""
gmm.py

混合ガウス分布のためのEMアルゴリズム
http://aidiary.hatenablog.com/entry/20100521/1274445325

"""

import numpy as np
from pylab import *

# 混合ガウス分布の数(固定)
K = 2


def scale(X):
    """
    データ行列Xを属性ごとに標準化したデータを返す
    :param X:
    :return:
    """
    # 属性の数(=列の数)
    col = X.shape[1]

    # 属性ごとに平均値と標準偏差を計算
    mu = np.mean(X, axis=0)
    sigma = np.std(X, axis=0)

    # 属性ごとにデータを標準化
    for i in range(col):
        X[:, i] = (X[:, i] - mu[i]) / sigma[i]  # TODO:Pysonicに改善あり

    return X


def gaussian(x, mean, cov):
    """
    多変量ガウス分布
    :param x: 観測データ
    :param mean: 平均
    :param cov: 分散
    :return:
    """
    A = 1. / ((2. * np.pi) ** (x.size / 2.0))
    B = 1. / (np.linalg.det(cov) ** 0.5)
    C = -0.5 * np.dot(np.dot(x - mean, np.linalg.inv(cov)), x - mean)
    return A * B * np.exp(C)


def likelihood(X, mean, cov, pi):
    """
    対数裕度関数
    :param X: 観測データ
    :param mean: 平均
    :param cov: 分散
    :param pi: 混合計数
    :return: 尤度
    """
    sum = 0.
    for n in range(len(X)):
        temp = 0.0
        for k in range(K):
            temp += pi[k] * gaussian(X[n], mean[k], cov[k])
        sum += np.log(temp)
    return sum


def main():
    # 訓練データをロード
    data = np.genfromtxt("faithful.txt")
    X = data[:, :]
    # データの標準化
    x = scale(X)
    # データ数
    N = len(x)

    # ---------------------------------------------
    # 訓練データから混合ガウス分布のパラメータを
    # EMアルゴリズムで推定する
    # ---------------------------------------------

    # 平均、分散、混合計数を初期化
    mean = np.random.rand(K, 2)
    cov = np.zeros((K, 2, 2))
    for k in range(K):
        cov[k] = [[1., 0.],
                  [0., 1.]]
    pi = np.random.rand(K)

    # 負担率の空配列を用意
    gamma = np.zeros((N, K))

    # 対数裕度の初期値を計算
    like = likelihood(x, mean, cov, pi)
    print 'like: ', like

    turn = 0

    MAX_ITERATE = 100
    THR_VALUE = 0.01
    for l in range(MAX_ITERATE):
        print 'like[%d] = %0.1f' % (l, like)

        # E-step : 現在のパラメータを使って、負担率を計算
        for n in range(N):
            # 分母はkによらないので最初に1回だけ計算
            denominator = 0.0
            for j in range(K):
                denominator += pi[j] * gaussian(x[n], mean[j], cov[j])
            # 各kについて負担率を計算
            for k in range(K):
                gamma[n, k] = pi[k] * gaussian(x[n], mean[k], cov[k]) / denominator

        # M-step : 現在の負担率を使って、パラメータを再計算
        for k in range(K):
            # Nkを計算
            Nk = 0.0
            for n in range(N):
                Nk += gamma[n, k]

            # 平均を再計算
            mean[k] = np.array([0., 0.])
            for n in range(N):
                mean[k] += gamma[n, k] * x[n]
            mean[k] /= Nk

            # 共分散を再計算
            cov[k] = np.array([[0., 0.], [0., 0.]])
            for n in range(N):
                # おそらくx[n,:] - mean[k]  (1xK)行列
                temp = x[n] - mean[k]
                # おそらく cov[k] += gamma[n,k] * np.dot(tmep.T, temp)
                cov[k] += gamma[n, k] * np.matrix(temp).reshape(2, 1) * np.matrix(temp).reshape(1, 2)
            cov[k] /= Nk

            # 混合計数を再計算
            pi[k] = Nk / N

        # 収束判定
        new_like = likelihood(x, mean, cov, pi)
        diff = new_like - like
        if diff < THR_VALUE:
            break
        like = new_like

    # # ガウス分布の平均を描画
    for k in range(K):
        scatter(mean[k, 0], mean[k, 1], c='r', marker='o')

    # 等高線を描画
    xlist = np.linspace(-2.5, 2.5, 50)
    ylist = np.linspace(-2.5, 2.5, 50)
    gx, gy = np.meshgrid(xlist, ylist)
    for k in range(K):
        gz = bivariate_normal(gx, gy, np.sqrt(cov[k, 0, 0]), np.sqrt(cov[k, 1, 1]), mean[k, 0], mean[k, 1], cov[k, 0, 1])
        cs = contour(gx, gy,gz, 3, colors='k', linewidths=1)

    # 訓練データを描画
    plot(x[:, 0], x[:, 1], 'gx')

    grid(True)
    savefig('gmm.png')
    show()
    pass


if __name__ == "__main__":
    main()

結果

gmm

使用しているデータセット

faithful.txt

その他データセットが欲しい方は

http://naoyat.hatenablog.jp/entry/2011/12/29/064419

俺の彼女と幼なじみが修羅場すぎる12 (GA文庫)

俺の彼女と幼なじみが修羅場すぎる12 (GA文庫)