PRML第9章 混合ガウス分布 Python実装

PRML第９章ではEMアルゴリズムが紹介されています。EMアルゴリズム自体は様々なところで使える手法で、私自身はEMアルゴリズムとある分類器を組み合わせてよりノイズに頑健なクラス分類を行ったりしたことがあります。しかし、EMアルゴリズムを適用する応用例として一番有名な混合ガウス分布によるクラスタリングをしたことがなかったので、今回は**EMアルゴリズムによる混合ガウス分布の最尤推定を実装**しました。

# 混合ガウス分布の最尤推定
![blobs.png](https://qiita-image-store.s3.amazonaws.com/0/148017/88f68e74-dd32-3c98-6336-651bef00998a.png)
例えば、上の図の青点のようなデータ点を普通のガウス分布

```math
p({\bf x}) = \mathcal{N}({\bf x}|{\bf\mu},{\bf\Sigma})
```
でフィッティングすると、
![gauss_fitting.png](https://qiita-image-store.s3.amazonaws.com/0/148017/37430ba8-acd5-fb9b-28b5-3740174bc313.png)
このようになってしまい、このケースでは単峰のガウス分布は適切なモデルではありません。データ点が３つの群に分かれていることに着目すると、今回の場合は３つのガウス分布を用いた混合ガウス分布

```math
p({\bf x}) = \sum_{k=1}^3 \pi_k\mathcal{N}({\bf x}|{\bf\mu}_k,{\bf\Sigma}_k)
```
が適切なモデルであると考えられます。ただし、$\sum_k\pi_k=1$とします。それぞれ${\bf\mu}\_1$が上、${\bf\mu}\_2$が右下、${\bf\mu}\_3$が左下のデータ点の塊の平均のようにします。よって**塊ごとにガウス分布でフィッティングすればいいのですが、どのデータ点がどの塊に属するのかは私たち人間には一目瞭然ですが、機械にはそんなことは分かりません**。

N個のデータ点の座標$\\{{\bf x}\_n\\}\_{n=1}^N$を観測変数として、どのデータ点がK個の塊のうちどれに属するかを1-of-k符号化で表す潜在変数$\\{{\bf z}\_n\\}\_{n=1}^N$とパラメータ$\\{{\bf\pi}\_k,{\bf\mu}\_k,{\bf\Sigma}\_k\\}\_{k=1}^K$を同時に推定します。このように潜在変数がある場合はEMアルゴリズムを用いるのが一般的です。


# EMアルゴリズム
PRMLの第9.2.2項に混合ガウス分布をEMアルゴリズムで最尤推定する手順（PRML式(9.23)~(9.27)）がまとめられているので割愛します。

# コード

## ライブラリ
いつもどおりにNumpyとmatplotlibを使います。

```python
import matplotlib.pyplot as plt
import numpy as np
```

## 混合ガウス分布

```python
class GaussianMixture(object):

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

    # EMアルゴリズムを用いた最尤推定
    def fit(self, X, iter_max=10):
        # データの次元
        self.ndim = np.size(X, 1)
        # 混合係数の初期化
        self.weights = np.ones(self.n_component) / self.n_component
        # 平均の初期化
        self.means = np.random.uniform(X.min(), X.max(), (self.ndim, self.n_component))
        # 共分散行列の初期化
        self.covs = np.repeat(10 * np.eye(self.ndim), self.n_component).reshape(self.ndim, self.ndim, self.n_component)

        # EステップとMステップを繰り返す
        for i in xrange(iter_max):
            params = np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))
            # Eステップ、負担率を計算
            resps = self.expectation(X)
            # Mステップ、パラメータを更新
            self.maximization(X, resps)
            # パラメータが収束したかを確認
            if np.allclose(params, np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))):
                break
        else:
            print("parameters may not have converged")

    # ガウス関数
    def gauss(self, X):
        precisions = np.linalg.inv(self.covs.T).T
        diffs = X[:, :, None] - self.means
        assert diffs.shape == (len(X), self.ndim, self.n_component)
        exponents = np.sum(np.einsum('nik,ijk->njk', diffs, precisions) * diffs, axis=1)
        assert exponents.shape == (len(X), self.n_component)
        return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim)

    # Eステップ
    def expectation(self, X):
        # PRML式(9.23)
        resps = self.weights * self.gauss(X)
        resps /= resps.sum(axis=-1, keepdims=True)
        return resps

    # Mステップ
    def maximization(self, X, resps):
        # PRML式(9.27)
        Nk = np.sum(resps, axis=0)

        # PRML式(9.26)
        self.weights = Nk / len(X)

        # PRML式(9.24)
        self.means = X.T.dot(resps) / Nk

        diffs = X[:, :, None] - self.means
        # PRML式(9.25)
        self.covs = np.einsum('nik,njk->ijk', diffs, diffs * np.expand_dims(resps, 1)) / Nk

    # 確率分布p(x)を計算
    def predict_proba(self, X):
        # PRML式(9.7)
        gauss = self.weights * self.gauss(X)
        return np.sum(gauss, axis=-1)

    # クラスタリング
    def classify(self, X):
        joint_prob = self.weights * self.gauss(X)
        return np.argmax(joint_prob, axis=1)
```

## 全体のコード

```python:gaussian_mixture.py
import matplotlib.pyplot as plt
import numpy as np


class GaussianMixture(object):

    def __init__(self, n_component):
        self.n_component = n_component

    def fit(self, X, iter_max=10):
        self.ndim = np.size(X, 1)
        self.weights = np.ones(self.n_component) / self.n_component
        self.means = np.random.uniform(X.min(), X.max(), (self.ndim, self.n_component))
        self.covs = np.repeat(10 * np.eye(self.ndim), self.n_component).reshape(self.ndim, self.ndim, self.n_component)
        for i in xrange(iter_max):
            params = np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))
            resps = self.expectation(X)
            self.maximization(X, resps)
            if np.allclose(params, np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))):
                break
        else:
            print("parameters may not have converged")

    def gauss(self, X):
        precisions = np.linalg.inv(self.covs.T).T
        diffs = X[:, :, None] - self.means
        assert diffs.shape == (len(X), self.ndim, self.n_component)
        exponents = np.sum(np.einsum('nik,ijk->njk', diffs, precisions) * diffs, axis=1)
        assert exponents.shape == (len(X), self.n_component)
        return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim)

    def expectation(self, X):
        resps = self.weights * self.gauss(X)
        resps /= resps.sum(axis=-1, keepdims=True)
        return resps

    def maximization(self, X, resps):
        Nk = np.sum(resps, axis=0)
        self.weights = Nk / len(X)
        self.means = X.T.dot(resps) / Nk
        diffs = X[:, :, None] - self.means
        self.covs = np.einsum('nik,njk->ijk', diffs, diffs * np.expand_dims(resps, 1)) / Nk

    def predict_proba(self, X):
        gauss = self.weights * self.gauss(X)
        return np.sum(gauss, axis=-1)

    def classify(self, X):
        joint_prob = self.weights * self.gauss(X)
        return np.argmax(joint_prob, axis=1)


def create_toy_data():
    x1 = np.random.normal(size=(100, 2))
    x1 += np.array([-5, -5])
    x2 = np.random.normal(size=(100, 2))
    x2 += np.array([5, -5])
    x3 = np.random.normal(size=(100, 2))
    x3 += np.array([0, 5])
    return np.vstack((x1, x2, x3))


def main():
    X = create_toy_data()

    model = GaussianMixture(3)
    model.fit(X, iter_max=100)
    labels = model.classify(X)

    x_test, y_test = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
    X_test = np.array([x_test, y_test]).reshape(2, -1).transpose()
    probs = model.predict_proba(X_test)
    Probs = probs.reshape(100, 100)
    colors = ["red", "blue", "green"]
    plt.scatter(X[:, 0], X[:, 1], c=[colors[int(label)] for label in labels])
    plt.contour(x_test, y_test, Probs)
    plt.xlim(-10, 10)
    plt.ylim(-10, 10)
    plt.show()


if __name__ == '__main__':
    main()
```

# 結果
点を学習データとして混合ガウス分布のパラメータを最尤推定して、その確率分布を等高線で図示したものです。また、点の色はどのクラスタに属するかを表しています。
![result.png](https://qiita-image-store.s3.amazonaws.com/0/148017/8f9c552e-04fe-8bf5-497c-e7d3e52e5a2e.png)
これは成功したときの結果ですが、**たまに失敗します**。PRMLにも記載されていますが、今回の対数尤度関数の最大化は不良設定問題になっていて良い解にならないことがあります。それを回避するヒューリスティックスもありますが、今回はその回避方法を実装していないのでたまに失敗します。といってもあまり失敗しないはずです。

# 終わりに
混合ガウス分布によるフィッティングを行うことで、教師なしのクラスタリングを行いました。その際に用いるガウス分布の個数をこちらで指定しました。次の第１０章では適切な混合ガウス分布の要素数も自動的に推定する手法が紹介されているので、次回はそこで用いられている変分ベイズを実装します。