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]]
ほぼ同じ値に落ち着いた。