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

# EMアルゴリズム実装してみた(?)
どうも19erのTeraponです。この記事は[**ISer Advent Calender 2018**](https://adventar.org/calendars/3598)の19日目として書かれたものです。昨日**kammer0820**さんの記事でした。面白そうな遊びがたくさん紹介されていたのでISのみんなで来年地下でやりたいですね。

さて、[**12日の記事**](https://qiita.com/tera_poon/items/03ca723c029c9a456e0e)に書いた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様のお力を借ります。

```python
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)型のテンソル積なので、[**公式ドキュメント**](https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html#numpy.einsum)曰くこのテンソル積の計算は**np.einsum('ijk,jlk->ilk', 要素１, 要素２)**とかくようです。
そして次の演算は（データの数, 変数の数, **self.n_gauss**)とその転置から（データの数, **self.n_gauss**)を得たいのですが、この計算は、アダマール積を計算して変数全てについて足し合わせることでなされるので、アダマール積を変数軸(axis=1)について足すことで計算します。（テンソル難しい)

```python
# ガウス関数の定義    
    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をとりましょう。そうすれば、自ずと正規化されます。

```python
# 初期化関数の定義
    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ステップ**を次に実装していきます。

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

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

```python
# 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関数を実装すれば終わりです。

```python
# 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
```
最後に全体を載せます。

```python
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