EMアルゴリズムによるGMMパラメータの推定

トピックモデルでは EM アルゴリズムを用いるのが一般的なのに,EM の理屈を理解できていない。
そこで,GMM (Gaussian Mixture Model) のパラメータ推定を対象とし,
Simon J.D. Prince, Computer vision: models, learning and inference (2012)
で勉強しつつの,
http://nbviewer.ipython.org/github/tritemio/notebooks/blob/master/Mixture_Model_Fitting.ipynb
に掲載されていたコードを打ち込んでみた。

#coding: utf-8

'''
下記サイトに掲載されていたコード
http://nbviewer.ipython.org/github/tritemio/notebooks/blob/master/Mixture_Model_Fitting.ipynb
'''

import numpy as np
from numpy import random
from scipy.optimize import minimize, show_options
import matplotlib.mlab as mlab

np.random.seed(1)

# 2つの正規分布を重み (0.3, 0.7) で混合する
N=1000
a=0.3
s1 = random.normal(0, 0.08, size=N*a)
s2 = random.normal(0.6, 0.12, size=N*(1-a))
s = np.concatenate([s1, s2])

'''
import matplotlib.pyplot as plt
plt.hist(s, bins=20)
plt.show()
'''

# matplotlib.mlab.normpdf を用いて,混合確率密度を得る
def pdf_model(x, p):
    mu1, sig1, mu2, sig2, pi_1 = p
    return pi_1 * mlab.normpdf(x, mu1, sig1) + (1-pi_1) * mlab.normpdf(x, mu2, sig2)

# 繰返し回数を固定(本来は値が収束するまで反復)
max_iter = 100

# パラメータの初期設定(テキトー?)
p0 = np.array([-0.2, 0.2, 0.8, 0.2, 0.5])
mu1, sig1, mu2, sig2, pi_1 = p0
mu = np.array([mu1, mu2])
sig = np.array([sig1, sig2])
pi_ = np.array([pi_1, 1-pi_1])

gamma = np.zeros((2, s.size))  # 2 * サンプル数のゼロ行列
N_ = np.zeros(2)   # サイズ 2 のゼロベクトル
p_new = p0

# ここから EM ループ...
counter = 0
converged = False
while not converged:
    # Compute the responsibility function and new parameters
    for k in [0, 1]:
        # E-step
        gamma[k,:] = pi_[k] * mlab.normpdf(s, mu[k], sig[k]) / pdf_model(s, p_new)

        # M-step
        N_[k] = 1. * gamma[k].sum()
        mu[k] = sum(gamma[k] * s) / N_[k]
        sig[k] = np.sqrt(sum(gamma[k] * (s - mu[k])**2) / N_[k])
        pi_[k] = N_[k] / s.size
    p_new = [mu[0], sig[0], mu[1], sig[1], pi_[0]]
    # assert abs(N_.sum() - N) / float(N) < 1e-6
    # assert abs(pi_.sum() - 1) < 1e-6
    counter += 1
    converged = counter >= max_iter

print "parameters : ", p_new

実行結果を以下に示す。

parameters :  [0.0063749972625032486, 0.076249720004177929, 0.60298724972290296, 0.11940589580298418, 0.30040038558122012]

サンプル生成の際に仮定した密度関数のパラメータを推定できているようだ。

次に scikit learn を用いた実装を試した。

import numpy as np
from numpy import random
from sklearn import mixture

np.random.seed(1)

# 2つの正規分布を重み (0.3, 0.7) で混合する
N=1000
a=0.3
s1 = random.normal(0, 0.08, size=N*a)
s2 = random.normal(0.6, 0.12, size=N*(1-a))
s = np.concatenate([s1, s2])

#s = np.vstack(s)

clf = mixture.GMM(n_components=2, covariance_type='diag')
clf.fit(s)

print "weight : ", clf.weights_
print "mean : ", clf.means_
print "sigma : ", np.sqrt(clf.covars_)

sklearn を用いた場合の実行結果を以下に示す。

weight :  [ 0.69922406  0.30077594]
mean :  [[ 0.60316138] [ 0.00671512]]
sigma :  [[ 0.12332669] [ 0.08305453]]

ほぼ同じ値に落ち着いた。