[Python] Implémentation du clustering à l'aide d'un modèle gaussien mixte

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()

En regardant la figure de droite, on peut interpréter que chaque donnée est générée selon la distribution gaussienne multivariée [^ 1], qui a des moyennes différentes pour chaque type d'iris. C'est le modèle gaussien mixte. Une moyenne approximative est indiquée dans le tableau ci-dessous.

En utilisant le modèle gaussien mixte, il est possible d'estimer la moyenne de plusieurs distributions gaussiennes à partir de l'ensemble de données non étiqueté illustré à gauche, et de corréler les données avec chaque distribution gaussienne pour le regroupement.

Explication mathématique

Le modèle gaussien mixte est décrit comme un modèle probabiliste. Soit le cluster de données $ x_i $ $ z_i $ et la distribution gaussienne multivariée correspondant à chaque cluster $ 1, \ cdots, k $ soit $ N (\ mu_k, \ Sigma_k) $. En d'autres termes, lorsque $ x_i $ appartient au clasper $ k $, sa probabilité de génération peut s'écrire: $ D $ est la dimension de $ x_i $ et $ N $ est le nombre de données. ..

python

\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}

D'autre part, $ z_i $ est également exprimé de manière probabiliste. Supposons que $ z_i $ soit généré à partir d'une distribution catégorielle avec $ \ pi: = (\ pi_1, \ cdots, \ pi_K) ^ T $ comme paramètres. $ K $ est le nombre de clusters.

python

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

D'après ce qui précède, la probabilité que le cluster $ z_i $ de données $ x_i $ soit $ k $ peut être écrite comme suit. [^ 2]

python

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)

Le clustering est effectué en estimant chaque $ z_i $ de l'ensemble de données à l'aide de cette formule.

Comment estimer le cluster

Pour estimer le cluster $ z_ {1: N} $ de chaque donnée à partir de la formule $ (1) $, les paramètres de la distribution gaussienne multivariée correspondant à chaque cluster $ \ mu_ {1: K}, \ Sigma_ { Vous devez connaître 1: K} $ et le paramètre de distribution catégorielle $ \ pi $. Celles-ci ne sont pas données explicitement et doivent être estimées à partir de l'ensemble de données ou supposées être des valeurs appropriées. Cette fois, nous estimons la moyenne $ \ mu_ {1: K} $ de la distribution gaussienne multivariée à partir de l'ensemble de données et supposons les paramètres restants aux valeurs suivantes: $ I_D $ est la matrice unitaire de $ D \ fois D $ et $ \ sigma ^ 2 $ est l'hyper paramètre. [^ 6]

python

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

Ensuite, pour estimer $ \ mu_ {1: K} $ et $ z_ {1: N} $, nous adopterons la méthode d'échantillonnage probabiliste de ceux-ci un par un. Cette méthode est appelée ** échantillonnage de Gibbs **. [^ 9] [^ 10]

Puisque l'échantillonnage de $ z_i $ peut être fait selon la formule $ (1) $, considérons la distribution de probabilité de $ \ mu_k $ ci-dessous. Je veux estimer la distribution de probabilité, il semble donc que l'estimation bayésienne puisse être utilisée. La moyenne $ \ mu_ {\ rm pri} $ de la distribution a priori conjuguée de $ \ mu_k $ et la matrice de covariance $ \ Sigma_ {\ rm pri} $ sont définies comme suit. [^ 7] Ceux-ci sont communs à tous les $ k = 1, \ cdots, K $. $ \ Sigma_ {\ rm pri} ^ 2 $ est un hyper paramètre.

python

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

A ce moment, la distribution a posteriori de $ \ mu_k $ est la suivante. $ n_k $ est le nombre de données appartenant au cluster $ k $ dans $ x_ {1: N} $, et $ \ bar {x} _k $ est leur moyenne.

python

\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}

Cette fois, $ \ Sigma_k $ et $ \ Sigma_ {\ rm pri} ^ {-1} $ sont tous deux des multiples constants de $ I_D $, donc $ \ mu_ {k, {\ rm pos}} $ et $ \ Sigma_ { k, {\ rm pos}} $ peut être transformé comme suit.

python

\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}

D'après ce qui précède, on peut voir que l'échantillonnage de $ \ mu_k $ doit être effectué selon l'équation (2).

la mise en oeuvre

Voici quelques points importants de la mise en œuvre.

Implémentation d'un tableau capable de compter le nombre de données dans le cluster

Il est nécessaire de calculer le nombre de données de chaque cluster à partir de la formule $ (2) $. Par conséquent, implémentez la classe de tableau ClusterArray avec une méthode count qui renvoie le nombre. [^ 8] En plus de la méthode count, les méthodes susceptibles d'être utilisées sont déléguées de numpy.ndarray.

python

import numpy as np
from collections import Counter

class ClusterArray(object):
    def __init__(self, array):
        #array est une liste unidimensionnelle, array
        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):
        #Soi lors de l'exécution._le compteur est également mis à jour
        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]

Supprimer les calculs supplémentaires

Vous pouvez calculer la fonction de densité de probabilité directement dans la formule $ (1) $, mais vous pouvez réduire le montant du calcul en supprimant les facteurs non liés à $ k $.

python

\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}

Implémentez la méthode log_deformed_gaussian pour calculer le contenu de $ \ exp $ dans cette formule et utilisez-la pour le calcul.

python

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

À propos de logsumexp

Pour les calculs comme $ \ log (\ sum_ {i = 1} \ exp f_i (x)) $, logsumexp est une technique pour éviter les débordements et sous-débordements. Cependant, ce n'est pas très efficace car il utilise $ \ log $ et $ \ exp $ plusieurs fois. Par conséquent, nous n'utiliserons pas logsumexp cette fois. [^ 5]

Pour logsumexp, l'article suivant sera très utile. Distribution gaussienne mixte et logsumexp

Mise en œuvre globale

J'ai essayé de le mettre en œuvre en gardant à l'esprit ce qui précède. Avec scikit-learn à l'esprit, le clustering peut être exécuté avec la méthode fit. Le code source est long, donc je le plie.

<détails>

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 est une liste unidimensionnelle, array
        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):
        #Soi lors de l'exécution._le compteur est également mis à jour
        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  #Nombre de clusters
        self.D = D  #Dimensions variables explicatives(Défini au moment du constructeur pour faciliter la mise en œuvre)
        self.z = None

        #Paramétrage de la distribution de probabilité
        self.mu = np.zeros((self.K, self.D))
        self.var = var  #Fixe, commun à tous les clusters
        self.pi = np.full(self.K, 1 / self.K)  #Fixe, commun à tous les clusters

        #Définition de la distribution antérieure
        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

import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline


#Charger le jeu de données
df = sns.load_dataset('iris')
X = df[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].values

#Clustering avec modèle gaussien mixte
gmc = GMMClustering(K=3, D=4, var=0.1, seed=1)
gmc.fit(X, n_iter=10)
df['GMM_cluster'] = gmc.z.array

#Visualisation des résultats
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()

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