LoginSignup
24
29

More than 3 years have passed since last update.

@isuya

【Python】混合ガウスモデルを使ったクラスタリングの実装

Last updated at Posted at 2020-01-16

機械学習プロフェッショナルシリーズ『ノンパラメトリックベイズ　点過程と統計的機械学習の数理』を読んで、本に書いてあるモデルをPythonで実装してみました。ノンパラメトリックベイズとは、混合ガウスモデルを拡張した手法です。本当はノンパラメトリックベイズを実装したいのですが、その準備として、混合ガウスモデルを実装します。

※Qiitaにて、同様のモデルや手法の説明、実装を説明されている記事がありましたので、紹介させていただきます。どちらも大変参考になりました。

※本記事に記載している数式には、書籍を参考に自分で計算したものも含まれています。

環境
Python: 3.7.5
numpy: 1.18.1
pandas: 0.25.3
matplotlib: 3.1.2
seaborn: 0.9.0

本記事で使われている記号

本編に入る前に、本記事で使用されている数学的記号のうち、他ではあまり見られないものを説明しておきます。

  • $x_{1:N}$ : $x_1, \cdots, x_N$の略
  • $x_{1:N}^{\backslash i}$ : $x_{1:N}$から$x_i$を除外したもの

また、本記事におけるベクトルは列ベクトルとし、太字で表記しません。

混合ガウスモデルとは

混合ガウスモデル(Gaussian Mixture Model, GMM)とは、複数のガウス分布から生成されたデータが混在しているデータセットのモデルです。アヤメのデータセットを例に説明します。これはアヤメの特徴と種類のデータセットであり、3種類のアヤメ「setosa」「versicolor」「virginica」が混在しています。機械学習の分類の練習問題によく使われるのでご存知の方も多いと思います。
このデータセットの説明変数は4つですが、簡略化のため「petal_length」「petal_width」の2つを取り出し、「色分けなしバージョン」と「アヤメの種類ごとに色分けしたバージョン」の散布図を作成します。

iris_valualization_1.png

ソースコード
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
df = sns.load_dataset('iris')
ax1.scatter(df['petal_length'], df['petal_width'], color='gray')
ax1.set_title('without label')
ax1.set_xlabel('petal_length')
ax1.set_ylabel('petal_width')
for sp in df['species'].unique():
    x = df[df['species'] == sp]['petal_length']
    y = df[df['species'] == sp]['petal_width']
    ax2.scatter(x, y, label=sp)
ax2.legend()
ax2.set_title('with label')
ax2.set_xlabel('petal_length')
ax2.set_ylabel('petal_width')
plt.show()

右図を見ると、アヤメの種類ごとに平均が異なる多変量ガウス分布1に従って、各データが生成されていると解釈できます。これが混合ガウスモデルです。大雑把な平均を下の表に示しました。

種類 petal_length petal_width
setosa 約1.5 約0.25
versicolor 約4.5 約1.3
virsinica 約6.0 約2.0

混合ガウスモデルを使うことで、ラベルのない左図のデータセットから複数のガウス分布の平均を推定し、データと各ガウス分布を紐付けてクラスタリングをすることができます。

数式的な説明

混合ガウスモデルは確率モデルとして記述されます。
データ$x_i$のクラスタを$z_i$、各クラスタ$1, \cdots, k$に対応する多変量ガウス分布を$N(\mu_k, \Sigma_k)$とします。つまり、$x_i$がクラスラ$k$に所属するとき、その生成確率は次のように書けます。$D$は$x_i$の次元で、$N$はデータの数です。。

\begin{align}
P(x_i | z_i=k, x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}) 
&= N(x | \mu_k, \Sigma_k)
\end{align}

一方、$z_i$も確率的に表現します。$z_i$を$\pi := (\pi_1, \cdots, \pi_K)^T$をパラメータとするカテゴリカル分布から生成されているとします。$K$はクラスタの数です。

P(z_i = k | x_k, \pi) = \pi_k

以上から、データ$x_i$のクラスタ$z_i$が$k$である確率は次のように書けます。2

P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi)
= \frac{\pi_k N(x_i| \mu_k, \Sigma_k)}{\sum_{l=1}^K \pi_l N(x_i| \mu_l, \Sigma_l)} \qquad\cdots (1)

データセットからこの式を使って各$z_i$を推定することで、クラスタリングを行います。

クラスタの推定方法

数式$(1)$より、各データのクラスタ$z_{1: N}$を推定するためには、各クラスタに対応する多変量ガウス分布のパラメータ$\mu_{1: K}, \Sigma_{1:K}$とカテゴリカル分布のパラメータ$\pi$を知っておく必要があります。これらは陽に与えられるものではないため、データセットから推定するか適切な値に仮定することになります。
今回、多変量ガウス分布の平均$\mu_{1: K}$をデータセットから推定し、残りのパラメータを次の値に仮定します。$I_D$は$D \times D$の単位行列で、$\sigma^2$はハイパパラメータです。3

\Sigma_1 = \cdots = \Sigma_K = \sigma^2 I_D \\
\pi = (1 / K, \cdots, 1 / K)^T

そして、$\mu_{1: K}$および$z_{1: N}$を推定するため、これらをひとつずつ確率的にサンプリングしていく方法を採用します。この方法をギブスサンプリングといいます。45

$z_i$のサンプリングは数式$(1)$に従えばいいので、以下、$\mu_k$の確率分布を考えます。
確率分布を推定したいので、ベイズ推定が使えそうですね。$\mu_k$の共役事前分布の平均$\mu_{\rm pri}$と共分散行列$\Sigma_{\rm pri}$を次の通りに定めます。6これらは、全ての$k=1, \cdots, K$で共通とします。$\sigma_{\rm pri}^2$はハイパパラメータです。

\begin{align}
\mu_{\rm pri} &= (0, \cdots, 0)^T \\
\Sigma_{\rm pri} &= \sigma_{\rm pri}^2I_D
\end{align}

このとき、$\mu_k$の事後分布は次のようになります。$n_k$は、$x_{1: N}$のうちクラスタ$k$に属するデータの数、$\bar{x}_k$はそれらの平均です。

\begin{align}
\mu_{k, {\rm pos}} &= \Sigma_{k, {\rm pos}}(n_k \Sigma_k^{-1} \overline{x}_k + \Sigma_{\rm pri}^{-1}\mu_{\rm pri})  \\
\Sigma_{k, {\rm pos}} &= (n_k \Sigma_{k}^{-1} + \Sigma_{\rm pri}^{-1})^{-1} \\
\mu_k &= N(\mu | \mu_{k, {\rm pos}}, \Sigma_{k, {\rm pos}}) \qquad\cdots (2)
\end{align}

今回は、$\Sigma_k$と$\Sigma_{\rm pri}^{-1}$がともに$I_D$の定数倍なので、$\mu_{k, {\rm pos}}$と$\Sigma_{k, {\rm pos}}$は次のように変形できます。

\begin{align}
\mu_{k, {\rm pos}} &= \Sigma_{k, {\rm pos}}\left(\frac{n_k}{\sigma^2} \overline{x}_k + \frac{1}{\sigma_{\rm pri}^2}\mu_{\rm pri}\right)  \\
\Sigma_{k, {\rm pos}} &= \left(\frac{n_k}{\sigma^2} + \frac{1}{\sigma_{\rm pri}^{2}}\right)^{-1}
\end{align}

以上より、$\mu_k$のサンプリングは数式(2)に従えばよいことが分かります。

実装

実装において重要なポイントをいくつかピックアップして説明します。

クラスタのデータ数をカウントできる配列の実装

数式$(2)$より、各クラスタが持っているデータの個数を計算する必要があります。そこで、その個数を返すcountメソッドが備わった配列クラスClusterArrayを実装します。7countメソッドの他、使いそうなメソッドをnumpy.ndarrayから委譲しています。

import numpy as np
from collections import Counter

class ClusterArray(object):
    def __init__(self, array):
        # arrayは1次元のリスト、配列
        self._array = np.array(array, dtype=np.int)
        self._counter = Counter(array)

    @property
    def array(self):
        return self._array.copy()

    def count(self, k):
        return self._counter[k]

    def __setitem__(self, i, k):
        # 実行されるとself._counterも更新される
        pre_value = self._array[i]
        if pre_value == k:
            return

        if self._counter[pre_value] > 0:
            self._counter[pre_value] -= 1
        self._array[i] = k
        self._counter[k] += 1

    def __getitem__(self, i):
        return self._array[i]

余分な計算の削除

数式$(1)$において、確率密度関数を直接計算してもいいのですが、$k$と無関係な因子を削除することで計算量を減らすことができます。

\begin{align}
P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi)
&= \frac{\pi_k \exp\{- \frac{1}{2}\log|\Sigma_k| - \frac{1}{2} (x_i - \mu_k)^T\Sigma_k(x_i - \mu_k)\}}{\sum_{l=1}^K \pi_l \exp\{- \frac{1}{2}\log|\Sigma_l|- \frac{1}{2} (x_i - \mu_l)^T\Sigma_l(x_i - \mu_l)\}}  \\
&= \frac{\pi_k \exp\{- \frac{1}{2 \sigma^2}\ \| x_i - \mu_k \|_2^2\}}{\sum_{l=1}^{K}\pi_l \exp\{- \frac{1}{2 \sigma^2}\ \| x_i - \mu_l \|_2^2\}}
\end{align}

この式の$\exp$の中身を計算するメソッドlog_deformed_gaussianを実装し、計算に使います。

def log_deformed_gaussian(x, mu, var):
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -norm_squared / (2 * var)

logsumexpについて

$\log(\sum_{i=1} \exp f_i(x))$のような計算において、オーバーフローやアンダーフローを防ぐためのテクニックにlogsumexpがあります。しかし、これは$\log$や$\exp$を何度も使うため計算の効率は良くないと思います。そこで、今回はlogsumexpを使用しません。8

logsumexpについては、次の記事が大変参考になります。
混合ガウス分布とlogsumexp

全体の実装

上記に留意しつつ実装してみました。scikit-learnを意識して、fitメソッドでクラスタリングが実行できるようにしています。
ソースコードは、長いので折りたたみにしています。

ソースコード
import numpy as np
from collections import Counter


def log_deformed_gaussian(x, mu, var):
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -norm_squared / (2 * var)


class ClusterArray(object):
    def __init__(self, array):
        # arrayは1次元のリスト、配列
        self._array = np.array(array, dtype=np.int)
        self._counter = Counter(array)

    @property
    def array(self):
        return self._array.copy()

    def count(self, k):
        return self._counter[k]

    def __setitem__(self, i, k):
        # 実行されるとself._counterも更新される
        pre_value = self._array[i]
        if pre_value == k:
            return

        if self._counter[pre_value] > 0:
            self._counter[pre_value] -= 1
        self._array[i] = k
        self._counter[k] += 1

    def __getitem__(self, i):
        return self._array[i]


class GaussianMixtureClustering(object):
    def __init__(self, K, D, var=1, var_pri=1, seed=None):
        self.K = K  # クラスタ数
        self.D = D  # 説明変数の次元(実装しやすたのため、コンストラクタの時点で設定しておく)
        self.z = None

        # 確率分布のパラメータ設定
        self.mu = np.zeros((self.K, self.D))
        self.var = var  # 固定、すべてのクラスタで共通
        self.pi = np.full(self.K, 1 / self.K)  # 固定、すべてのクラスタで共通

        # 事前分布の設定
        self.mu_pri = np.zeros(self.D)
        self.var_pri = var_pri

        self._random = np.random.RandomState(seed)

    def fit(self, X, n_iter):
        init_z = self._random.randint(0, self.K, X.shape[0])
        self.z = ClusterArray(init_z)

        for _ in range(n_iter):
            for k in range(self.K):
                self.mu[k] = self._sample_mu_k(X, k)
            for i, x_i in enumerate(X):
                self.z[i] = self._sample_zi(x_i)

    def _sample_zi(self, x_i):
        log_probs_xi = log_deformed_gaussian(x_i, self.mu, self.var)

        probs_zi = np.exp(log_probs_xi) * self.pi
        probs_zi = probs_zi / probs_zi.sum()

        z_i = self._random.multinomial(1, probs_zi)
        z_i = np.where(z_i)[0][0]
        return z_i

    def _sample_mu_k(self, X, k):
        xk_bar = np.array([x for i, x in enumerate(X) if self.z[i] == k]).mean(axis=0)
        var_pos = 1 / (self.z.count(k) / self.var + 1 / self.var_pri)
        mu_pos = var_pos * (xk_bar * self.z.count(k) / self.var + self.mu_pri / self.var_pri)

        mu_k = self._random.multivariate_normal(mu_pos, var_pos * np.eye(self.D))
        return mu_k

クラスタリングを試してみる

実装した混合ガウスモデルを使って、冒頭のアヤメのデータセットをクラスタリングしてみます。ハイパパラメータは、$\sigma^2=0.1, \sigma_{\rm pri}^2=1$とし、ギブスサンプリングの反復回数は10回にしました。
実際のデータセットのラベルとクラスタリング結果を比較すると、次のようになりました。

clustering_valualization.png

ソースコード
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline


# データセットの読み込み
df = sns.load_dataset('iris')
X = df[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].values

# 混合ガウスモデルによるクラスタリング
gmc = GMMClustering(K=3, D=4, var=0.1, seed=1)
gmc.fit(X, n_iter=10)
df['GMM_cluster'] = gmc.z.array

# 結果の可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
for sp in df['species'].unique():
    x = df[df['species'] == sp]['petal_length']
    y = df[df['species'] == sp]['petal_width']
    ax1.scatter(x, y, label=sp)
ax1.legend()
ax1.set_title('species')
ax1.set_xlabel('petal_length')
ax1.set_ylabel('petal_width')
for k in range(gmc.K):
    x = df[df['GMM_cluster'] == k]['petal_length']
    y = df[df['GMM_cluster'] == k]['petal_width']
    ax2.scatter(x, y, label=k)
ax2.legend()
ax2.set_title('GMM cluster')
ax2.set_xlabel('petal_length')
ax2.set_ylabel('petal_width')
plt.show()

「versicolor」と「virginica」の境目が怪しいですが、ラベル通りのクラスタリングが概ねできています。
また、seabornのpairplotを使って、クラスタリング結果を可視化してみました。

clustering_valualization_2.png

ソースコード
sns.pairplot(
    df.drop(columns=['species']),
    vars=['sepal_length', 'sepal_width', 'petal_length', 'petal_width'],
    hue='GMM_cluster'
)

いい感じにクラスタリングできています。pairplotの対角線上に並んでいる図を見ると、データセットが混合ガウスモデルで表現されていることが分かると思います。

おわりに

機械学習プロフェッショナルシリーズ『ノンパラメトリックベイズ　点過程と統計的機械学習の数理』から混合ガウスモデルを実装し、簡単なデータセットでクラスタリングを試してみました。
今回は混合ガウスモデルを扱いましたが、「混合ベルヌーイモデル」や「混合ポアソンモデル」のようなモデルも考えることができ、クラスタリングに利用することができます。
次回は、周辺化ギブスサンプリングについて書く予定です。これは、ノンパラメトリックベイズを行う上で必要な技術です。

(2020/1/21追記)
次回記事できました。
【Python】周辺化ギブスサンプリングを実装してみた

  1. 分散も異なりますが、簡略化のため平均のみ考えています。 ↩

  2. ベイズの定理を使って導出することができます。 ↩

  3. $\sigma^2$は分散です。共分散はすべて$0$とします。雑な仮定に見えますが、$\sigma^2$を適切に決めることでこれでも十分クラスタリングできますし、簡略化したことで計算量を減らすことができます。当然、これらもデータセットから推定することで、より正確なクラスタリングが可能になるはずです。 ↩

  4. ギブスサンプリングは解析的な計算が困難な同時分布を近似するための手法です。今回は、$P(\mu_{1 : K}, z_{1: N} | \cdots)$を近似しています。 ↩

  5. 混合ガウス分布の推定には他の手法も使えますが、ノンパラメトリックベイズに拡張するためにはギブスサンプリングが必要だと認識しています。正直ちゃんと理解できていないので、分かる方いたら教えて下さい。 ↩

  6. $\mu_k$の共役事前分布は多変量ガウス分布です。 ↩

  7. わざわざこのようなクラスを作らずに、Counterインスタンスを外に保持してもいいと思います。私はこういうのをカプセル化してクラスを作るのが好きです。 ↩

  8. データセットに対して大きすぎず、小さすぎない$\sigma^2$を設定すればオーバーフローとアンダーフローを回避できるので、logsumexpを使わなくても問題ないと判断しました。 ↩

24
29
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
24
29