Edited at

EMアルゴリズム実装してみた(?)


EMアルゴリズム実装してみた(?)

どうも19erのTeraponです。この記事はISer Advent Calender 2018の19日目として書かれたものです。昨日kammer0820さんの記事でした。面白そうな遊びがたくさん紹介されていたのでISのみんなで来年地下でやりたいですね。

さて、12日の記事に書いたEMアルゴリズムの実装をしようとしましたが間に合わない(ごめんなさい)のでクラスを作るところまでを載せていきたいと思います。とりあえずEMアルゴリズムのまとめを再掲したいと思います。


EMアルゴリズムのまとめ


  • ①パラメータを初期化する。

$$ \boldsymbol{\mu_{k}}=\boldsymbol{\mu_{k0}} $$

$$ \boldsymbol{\Sigma_{k}}=\boldsymbol{\Sigma_{k0}} $$

$$ \pi_{k}=\pi_{k0} $$


  • ②Eステップ:負担率を計算する。

$$ \gamma(z_{nk})=\frac{\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})}{\sum_{j=1}^{N}\pi_{j}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{j}},\boldsymbol{\Sigma_{j}})} $$


  • ③Mステップ:負担率を用いて3つのパラメータの値を上から順番に計算する。

$$ \boldsymbol{\mu_{k}}^{new}=\frac{1}{N_{k}}\sum_{n=1}^{N}\gamma(z_{nk})\boldsymbol{x_{n}} $$

$$ \boldsymbol{\Sigma_{k}}^{new}=\frac{1}{N_{k}}\sum_{n=1}^{N}\gamma(z_{nk})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}^{new}})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}^{new}})^{T} $$

$$ \pi_{k}^{new}=\frac{N_{k}}{N} $$

ただし、

$$ N_{k}=\sum_{n=1}^{N}\gamma(z_{nk}) $$

である。


  • ④収束の判定(初回はここには来ずに②に戻る)

$$ ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=\sum_{n=1}^{N}ln\,\Biggl(\sum_{k=1}^{K}\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})\Biggr) $$

を計算して、その変化量と閾値$\delta$との大小を判定。$\delta$より大きければ②に戻る。


 実装の流れ

まずはガウス分布を作ります。numpy.linalgは行列計算を行うpackageです。行列計算の実装は今回の肝ではないのでnumpy様のお力を借ります。

import numpy as np

from numpy import linalg as la

ガウス分布をいくつ混合するか、その数はself.n_gaussに格納します。

平均はself.muに格納します。

Xの形を(データの数, 変数の数)とすると、普通のガウス分布ではself.muは(1, 変数の数)という形をとりますが、ここでは3つ目の軸にself.n_gaussをとってテンソルにしたいと思います。つまり、(1, 変数の数, self.n_gauss)という形にします。

共分散行列はself.sigmaに格納します。普通のガウス分布ではself.sigmaは(変数の数、変数の数)という形をとりますが、同様にしてテンソル化し、(変数の数, 変数の数, n_gauss)とします。

指数部分の計算はテンソル積になります。データの方がわけわからなくなるのでテスト文も書いておきます。

テンソル積で掛け合わせるデータですが、初めの二つは(データの数, 変数の数, self.n_gauss)×(変数の数, 変数の数, self.n_gauss)=(データの数, 変数の数, n_gauss)としたいわけです。

これは、(i,j,k)×(j,l,k)=(i,l,k)型のテンソル積なので、公式ドキュメント曰くこのテンソル積の計算はnp.einsum('ijk,jlk->ilk', 要素1, 要素2)とかくようです。

そして次の演算は(データの数, 変数の数, self.n_gauss)とその転置から(データの数, self.n_gauss)を得たいのですが、この計算は、アダマール積を計算して変数全てについて足し合わせることでなされるので、アダマール積を変数軸(axis=1)について足すことで計算します。(テンソル難しい)

# ガウス関数の定義    

def gauss(self, X):

# データの数
data_num = len(X)

# 変数の数
variable_num = X.shape[1]

# 逆共分散行列=精度
sigma_inv = la.inv(self.sigma.T).T

# 平均からの変位ベクトル(Xのデータ型にn_gaussが含まれないのでX[:,:,None]として拡張した)
self.gap = X[:,:,None] - self.mu

# 共分散行列式
sigma_det = la.det(self.sigma.T).T

# テンソル積の部分
exp = np.sum(np.einsum('ijk,jlk->ilk', gap, sigma_inv)*gap, axis=1)

# 分母
denomin = np.sqrt(sigma_det*(2*np.pi)*self.variables)

#混合係数行列で重みづけして返して終了
return self.pi*(np.sqrt(-0.5*exp)/denomin)

さて、次は、①パラメータを初期化するを作っていきます。self.muについては、通常は想定される値の範囲のなかで一様分布からランダムに抽出します。今回は入ってくるデータに応じて値のスケールが変わってくるので、入力Xのmaxとminの間から(データの数, self.n_gauss)だけとってきます。

self.sigmaに関しては、ここでは全体の分散を中心にして0.5倍から1.5倍までの範囲で一様分布からランダムにとってきましょう。

self.pi(混合係数)に関しては普通に1/Nをとりましょう。そうすれば、自ずと正規化されます。

# 初期化関数の定義

def __init__(self, n_gauss, X):
# ガウス分布の個数
self.n_gauss = n_gauss

self.X = X

# データの数
self.data_num = len(self.X)

# 変数の数
self.variable_num = self.X.shape[1]

# 平均値ベクトルの初期化
self.mu = np.random.uniform(self.X.max(), self.X.min(), (self.data_num, self.n_gauss))

# 共分散行列の初期化
self.sigma = np.random.uniform(0.5*np.var(self.X), 1.5*np.var(self.X), (self.variable_num, self.variable_num, self.n_gauss))

# 混合係数の初期値(=1/n_gauss、各ガウス分布の重要度は全部同じにして、合計を1にする)
self.pi = np.ones(self.n_gauss)/self.n_gauss

②Eステップを次に実装していきます。

# Eステップ

def e_step(self, X):
self.gamma = self.gauss(X)/self.gauss(X).sum(axis=-1)

③Mステップも実装していきます。ここら辺は式通り書くだけなので楽です。

# Mステップ

def m_step(self, X):
self.N_k = np.sum(self.gamma)

# self.muの更新
self.mu = np.dot(X.T, self.gamma)/self.N_k

# self.sigmaの更新(アダマール積と内積の組み合わせによりテンソル積を計算)
self.sigma = np.einsum('ijk,ilk->jlk', self.gap, self.gap*np.expand_dims(self.gamma, axis=1))

#self.piの更新
self.pi = self.N_k/self.data_num

最後にEMアルゴリズムを回すfit関数を実装すれば終わりです。

# fit関数   

def fit(self, X):
while True:
# 現在の対数尤度関数の値
previous = np.sum(np.ln(np.sum(self.gauss(), axis=-1)))

# Eステップ->Mステップ
self.e_step(X)
self.m_step()

# 収束判定(さっきの対数尤度関数と今新しくできた対数尤度関数の差を見る)
now = np.sum(np.ln(np.sum(self.gauss(), axis=-1)))
if np.allclose(previous, now):
break

最後に全体を載せます。

import numpy as np

# numpyの線形代数を扱うpackage linear algebra(linalg)をimport
from numpy import linalg as la

# EMAlgorismを動かすclassをdefine
class EMAlgorism(object):

# 初期化関数の定義
def __init__(self, n_gauss, X):
# ガウス分布の個数
self.n_gauss = n_gauss

self.X = X

# データの数
self.data_num = len(self.X)

# 変数の数
self.variable_num = self.X.shape[1]

# 平均値ベクトルの初期化
self.mu = np.random.uniform(self.X.max(), self.X.min(), (self.data_num, self.n_gauss))

# 共分散行列の初期化
self.sigma = np.random.uniform(0.5*np.var(self.X), 1.5*np.var(self.X), (self.variable_num, self.variable_num, self.n_gauss))

# 混合係数の初期値(=1/n_gauss、各ガウス分布の重要度は全部同じにして、合計を1にする)
self.pi = np.ones(self.n_gauss)/self.n_gauss

# ガウス関数の定義
def gauss(self, X):

# データの数
data_num = len(X)

# 変数の数
variable_num = X.shape[1]

# 逆共分散行列=精度
sigma_inv = la.inv(self.sigma.T).T

# 平均からの変位ベクトル(Xのデータ型にn_gaussが含まれないのでX[:,:,None]として拡張した)
self.gap = X[:,:,None] - self.mu

# 共分散行列式
sigma_det = la.det(self.sigma.T).T

# テンソル積の部分
exp = np.sum(np.einsum('ijk,jlk->ilk', gap, sigma_inv)*gap, axis=1)

# 分母
denomin = np.sqrt(sigma_det*(2*np.pi)*self.variables)

assert sigma_inv.shape == (variable_num, variable_num, self.n_gauss)
assert self.gap.shape == (data_num, variable_num, self.n_gauss)
assert exp.shape == (data_num, self.n_gauss)

#混合係数行列で重みづけして返して終了
return self.pi*(np.sqrt(-0.5*exp)/denomin)

# Eステップ
def e_step(self, X):
self.gamma = self.gauss(X)/self.gauss(X).sum(axis=-1)

# Mステップ
def m_step(self, X):
self.N_k = np.sum(self.gamma)

# self.muの更新
self.mu = np.dot(X.T, self.gamma)/self.N_k

# self.sigmaの更新(アダマール積と内積の組み合わせによりテンソル積を計算)
self.sigma = np.einsum('ijk,ilk->jlk', self.gap, self.gap*np.expand_dims(self.gamma, axis=1))

#self.piの更新
self.pi = self.N_k/self.data_num

# fit関数
def fit(self, X):
while True:
# 現在の対数尤度関数の値
previous = np.sum(np.ln(np.sum(self.gauss(X), axis=-1)))

# Eステップ->Mステップ
self.e_step(X)
self.m_step(X)

# 収束判定
now = np.sum(np.ln(np.sum(self.gauss(X), axis=-1)))
if np.allclose(previous, now):
break


感想

せっかく作ったので動かしたいのですが、締め切り(12/19)が近づいてきました。(なめていた。テンソル難しい、時間取られた、着手が遅すぎた(12/17~)) 悲しいので暇なときにこれを動かす記事を書きたいと思います。中途半端になって申し訳ないです...m(_ _)m

明日はyoshiki_FFFさんの記事です。ありがとうございました。


参考にしたサイトなど

https://www.sejuku.net/blog/70006

https://deepage.net/features/numpy-dot.html

https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html#numpy.einsum

https://mathtrain.jp/tahenryogauss

http://save.sys.t.u-tokyo.ac.jp/~kawai/main_math/node226.html

https://teratail.com/questions/124492

http://www.procrasist.com/entry/einsum

http://python-remrin.hatenadiary.jp/entry/2017/04/26/233717

https://pythondatascience.plavox.info/numpy/%E7%B7%9A%E5%BD%A2%E4%BB%A3%E6%95%B0

https://docs.scipy.org/doc/numpy/reference/generated/numpy.expand_dims.html

http://oppython.hatenablog.com/entry/2015/05/29/020902