Go to Qiita Advent Calendar 2024 Top
LoginSignup
5
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@tera_poon

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

Last updated at Posted at 2018-12-18

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ステップ：負担率を用いて３つのパラメータの値を上から順番に計算する。

$$ \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, 変数の数)という形をとりますが、ここでは３つ目の軸に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', 要素１, 要素２)**とかくようです。
そして次の演算は（データの数, 変数の数, 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に関しては、ここでは全体の分散を中心にして０.５倍から１.５倍までの範囲で一様分布からランダムにとってきましょう。
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、各ガウス分布の重要度は全部同じにして、合計を１にする)
        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、各ガウス分布の重要度は全部同じにして、合計を１にする)
        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

5
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?