[Python] Clustering avec un modèle gaussien infiniment mélangé

L'une des méthodes de regroupement consiste à estimer un modèle gaussien mixte. Il s'agit d'une méthode pour exprimer et regrouper un ensemble de données avec plusieurs distributions gaussiennes. Normalement, ces distributions gaussiennes sont supposées être finies, mais s'il y avait un nombre infini, quel genre de regroupement serait-il?

Dans cet article, en référence à Machine Learning Professional Series «Non-parametric Bayes Point Process and Statistical Machine Learning Mathematics» Présentation du contour et de la mise en œuvre du «modèle gaussien infiniment mélangé» et «non paramétrique». Je recommande ce livre car il est très facile à comprendre.

Je voudrais omettre la théorie détaillée et me concentrer sur le plan et l'histoire. En outre, l'implémentation doit être limitée au code source.

Article précédent [Python] J'ai implémenté un échantillonnage de Gibbs marginalisé Cependant, je suis conscient que le flux peut être véhiculé même si je ne lis que cet article.

Introduction du modèle gaussien mixte infini

Dans un modèle gaussien mixte normal, l'ensemble de données est représenté par un nombre fini de distributions gaussiennes, et chaque donnée est associée à la distribution gaussienne qu'il contient. Le nombre de claspers (= nombre de distributions gaussiennes) doit être déterminé à l'avance, mais le nombre approprié de clusters varie en fonction de l'ensemble de données, et il n'est pas facile de le saisir.

Si le nombre de clusters est automatiquement déterminé lors de la mise en cluster, vous n'aurez pas à vous soucier des problèmes ci-dessus. Le ** modèle gaussien infiniment mélangé ** rend cela possible.

Le modèle gaussien infiniment mélangé est basé sur l'idée qu'il existe un nombre infini de distributions gaussiennes et que les données de l'ensemble de données sont associées à un nombre fini d'entre elles. La figure suivante montre le résultat du clustering, qui sera montré plus tard. (Parce que c'est gênant, je vais utiliser à nouveau le chiffre. Veuillez comprendre.)

Cet ensemble de données est représenté par une distribution gaussienne de trois. Alors qu'un modèle gaussien mixte normal pense qu'il y a exactement trois distributions gaussiennes dans ce monde, un modèle gaussien mixte infini pense qu'il existe un nombre infini de distributions gaussiennes potentielles qui ne sont pas liées à un ensemble de données. En pensant de cette manière, en plus de la probabilité que les données soient connectées à trois clusters existants, la probabilité que les données soient connectées à un cluster inconnu où la moyenne est inconnue peut être calculée.

En d'autres termes, dans le clustering avec le modèle gaussien infiniment mélangé, les données sont affectées à un nouveau cluster, tandis que les clusters existants disparaissent sans être connectés à aucune donnée, augmentant ou diminuant le nombre de clusters, et finalement convergeant vers un nombre approprié. Vous pouvez.

Calcul de la distribution de probabilité

Je publierai la formule de calcul à partir d'ici. Dans une distribution gaussienne mixte normale, la *** distribution catégorielle *** est appliquée à la distribution de probabilité de $ z_i $ (lorsque $ x_ {1: N} $ [^ 4] n'est pas observé).

`python`


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

Où $ \ pi: = (\ pi_1, \ cdots, \ pi_K) $ est un point sur la dimension $ (K-1) $ seule. [^ 1] Comme méthode de détermination de la valeur de $ \ pi $, il existe une méthode pour supposer une distribution uniforme ou une méthode d'estimation bayésienne en supposant une distribution a priori [^ 2].

D'autre part, le modèle gaussien infiniment mélangé applique un processus stochastique appelé ** processus de cuisson chinois **. Il s'agit d'un mécanisme pour calculer la distribution de probabilité du cluster $ z_i $ des $ i $ th données du cluster $ z_ {1: N} ^ {\ backslash i} $ des données autres que les $ i $ th.

`python`


P(z_i = k | z_{1:N}^{\backslash i}, \alpha) =
\left\{
\begin{array}{ll}
\frac{n_k}{N - 1 + \alpha} & (k = 1, \cdots, K) \\
\frac{\alpha}{N - 1 + \alpha} & (k = K + 1)
\end{array}
\right.

$ k = 1, \ cdots, K $ correspond au cluster existant et $ k = K + 1 $ correspond au nouveau cluster. Ici, $ n_k $ est le nombre de $ j $ tel que $ i \ neq j, z_j = k $, et $ \ alpha $ est un hyper paramètre qui prend une valeur positive. Comme vous pouvez le voir dans la formule, plus $ \ alpha $ est grand, plus il est probable qu'il soit affecté à un nouveau cluster. À partir de maintenant, j'écrirai cette distribution de probabilité sous la forme $ {\ rm CRP} (k | z_ {1: N} ^ {\ backslash i}, \ alpha) $.

Ensuite, considérons la distribution de probabilité lors de l'observation de $ x_ {1: N} $. Dans un modèle gaussien mixte normal, il peut être calculé comme suit.

`python`


\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, \pi) 
&= \frac{P(z_i = k | \pi) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}P(z_i = l | \pi) N(x_i | \mu_l, \Sigma_l)} \\
&= \frac{\pi_k N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K}\pi_l N(x_i | \mu_l, \Sigma_l)}
\end{align}

$ N (x | \ mu, \ Sigma) $ est la fonction de densité de probabilité de la distribution gaussienne multivariée. Dans le modèle gaussien infiniment mélangé, la distribution de probabilité peut être obtenue en changeant $ P (z_i = k | \ pi) $ dans l'équation ci-dessus en l'équation pour le processus de cuisson chinois.

`python`


\begin{align}
P(z_i = k | x_{1:N}, \mu_{1:K}, \Sigma_{1:K}, z_{1:N}^{\backslash i}, \alpha) 
&= \frac{{\rm CRP}(k | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_k, \Sigma_k)}{\sum_{l=1}^{K+1}{\rm CRP}(l | z_{1:N}^{\backslash i}, \alpha) N(x_i | \mu_l, \Sigma_l)}
\end{align}

Il y a une mise en garde ici. Pour les clusters existants $ k = 1, \ cdots, K $, il est possible d'estimer la moyenne $ \ mu_k $ et la matrice distribuée co-distribuée $ \ Sigma_k $ à partir des données appartenant à chaque cluster. Cependant, pour le nouveau cluster $ k = K + 1 $, il n'est pas possible d'estimer directement la moyenne et la matrice distribuée co-distribuée car il n'y a pas de données appartenant au cluster. Par conséquent, cette fois, nous adopterons la méthode suivante.

-Supposant une distribution a priori pour $ \ mu_k $, marginalisant la distribution simultanée $ P (x_i, \ mu_k | \ cdots) $, et trouvant la distribution de probabilité sans estimer explicitement $ \ mu_k $ -Soit $ \ Sigma_k $ la matrice commune $ \ sigma ^ 2 I $ pour tous les clusters. ($ \ Sigma ^ 2 $ est un hyper paramètre)

Pour la moyenne, l'estimation bayésienne et la marginalisation calculent la distribution de probabilité uniquement à partir de la distribution antérieure. En ce qui concerne la matrice de co-distribution de distribution, la fixer comme ci-dessus est suffisante pour le clustering, nous ne l'évaluerons donc pas cette fois. [^ 5] Bien entendu, il est également possible d'estimer la variance et de la regrouper.

Comme il est absurde de séparer les cas, cette méthode est également appliquée lorsque $ k = 1, \ cdots, K + 1 $. Veuillez vous référer à Article précédent pour la marginalisation.

Clustering avec des baies non paramétriques

Afin d'estimer le cluster $ z_ {1: N} $ de $ N $ données, nous échantillonnerons de manière probabiliste un par un de $ z_1 $ à $ z_N $. Cette méthode est appelée ** échantillonnage de Gibbs **, et la méthode basée sur la distribution de probabilité dimensionnelle infinie comme cette fois est particulièrement appelée ** baies non paramétriques **.

Pour exécuter cette méthode, nous devons faire quelques hypothèses sur le clustering initial $ z_ {1: N} $. Cette fois, avec le même cluster $ k = 1 $ que le clustering initial, l'échantillonnage de $ z_1 $ à $ z_N $ sera répété plusieurs fois. Au départ, il n'y a qu'un seul cluster, mais c'est une image que des données éloignées de la moyenne sont combinées avec un nouveau cluster et un anneau de classe progresse tout en augmentant le nombre de clusters.

la mise en oeuvre

Implémentation du clustering avec un modèle gaussien infiniment mélangé. Il est basé sur l'implémentation de Dernière fois, mais il a été globalement repensé et la différence peut être importante.

<détails> <résumé> Code source </ résumé>

`python`


from collections import Counter
from functools import partial
from typing import Sequence, Tuple
import numpy as np
from scipy.special import logsumexp


def _log_gaussian(x: np.ndarray, mu: np.ndarray, var: float) -> np.ndarray:
    #2 dans l'article précédent*np.J'ai omis pi, mais par souci de clarté, je vais en tenir compte cette fois.
    d = x.shape[0]
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return - d * np.log(2 * np.pi * var) / 2 - norm_squared / (2 * var)


#Classe Array pour stocker le cluster
class InfiniteClusterArray(object):
    def __init__(self, array: Sequence[int]):
        self._array = np.array(array, dtype=int)
        self._counter = Counter(array)
        self._n_clusters = max(array) + 1

    def update(self, i: int, k: int):
        assert k <= self._n_clusters

        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

        if k == self._n_clusters:
            self._n_clusters += 1

        if not self._counter[pre_value]:
            self._array -= 1 * (self._array > pre_value)
            self._counter = Counter(self._array)
            self._n_clusters -= 1

    @property
    def array(self) -> np.ndarray:
        return self._array.copy()

    @property
    def n_clusters(self) -> int:
        return self._n_clusters

    def count(self, k: int, excluded=None) -> int:
        assert (excluded is None) or isinstance(excluded, int)
        if excluded is None:
            return self._counter[k]

        counted = self._counter[k] - bool(self._array[excluded] == k)
        return counted

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

            

#Classe Array pour stocker l'ensemble de données
class DataArray(object):
    def __init__(self, array: np.ndarray):
        assert array.ndim == 2
        self._array = array

    @property
    def n(self) -> int:
        return self._array.shape[0]

    @property
    def d(self) -> int:
        return self._array.shape[1]

    def mean(self, excluded=None) -> np.ndarray:
        assert (excluded is None) or isinstance(excluded, Sequence)

        if excluded is None:
            return self._array.mean(axis=0)
        else:
            return self._array[excluded].mean(axis=0)

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

    def __iter__(self):
        for i in range(self.n):
            yield self._array[i]

    def __str__(self) -> str:
        return f'DataArray({self._array})'


#Classe d'échantillonnage Gaussian Gibbs infiniment mélangée
class IGMMGibbsSampling(object):
    def __init__(self, var: float, var_pri: float, alpha: float, seed=None):
        #La moyenne de la distribution précédente est définie sur la moyenne de l'ensemble de données.
        assert (var > 0) and (var_pri > 0)
        assert alpha > 0

        self.var = var
        self.var_pri = var_pri
        self.alpha = alpha
        self._random = np.random.RandomState(seed)

    def predict_clusters(self, X: np.ndarray, n_iter: int, init_z=None) -> np.ndarray:
        #Méthode d'estimation du cluster
        assert X.ndim == 2
        X = DataArray(X)
        mu_pri = X.mean()

        if init_z is None:
            init_z = np.zeros(X.n, dtype=int)
        z = InfiniteClusterArray(init_z)

        for _ in range(n_iter):
            compute_probs = partial(self._compute_pmf_zi, X, z, mu_pri)
            for i in range(X.n):
                probs = compute_probs(i)
                z_i = self._random.multinomial(1, probs)
                z_i = np.where(z_i)[0][0]  # One-hot to index
                z.update(i, z_i)

        return z.array

    def _compute_pmf_zi(self, X: DataArray, z: InfiniteClusterArray, mu_pri: np.ndarray, i: int) -> np.ndarray:
        mu, var = self._compute_marginal_distribution(X, z, mu_pri, i)
        pi = np.array([
            z.count(k, excluded=i) / (X.n - 1 + self.alpha)
            if k < z.n_clusters
            else self.alpha / (X.n - 1 + self.alpha)
            for k in range(z.n_clusters + 1)
        ])

        log_pdf_xi = _log_gaussian(X[i], mu, var)
        with np.errstate(divide='ignore'):
            pmf_zi = np.exp(np.log(pi) + log_pdf_xi - logsumexp(np.log(pi) + log_pdf_xi)[np.newaxis])

        return pmf_zi

    def _compute_marginal_distribution(self, X: DataArray, z: InfiniteClusterArray,
                                       mu_pri: np.ndarray, i: int) -> Tuple[np.ndarray, float]:

        n_vec = np.array([
            z.count(k, excluded=i)
            for k in range(z.n_clusters + 1)
        ], dtype=np.int)

        tmp_vec = 1 / (n_vec / self.var + 1 / self.var_pri)

        mu = np.tile(mu_pri / self.var_pri, (z.n_clusters + 1, 1))  # Shape (z.n_clusters + 1, X.d)
        for k in range(z.n_clusters):
            excluded_list = [
                j for j in range(X.n)
                if (j != i) and (z[j] == k)
            ]
            if excluded_list:
                mean = X.mean(excluded=excluded_list)
                mu[k] += n_vec[k] * mean / self.var
        mu *= tmp_vec[:, np.newaxis]
        var = tmp_vec + self.var

        return mu, var

Pratique de clustering

À titre de test, essayez le clustering en utilisant la classe auto-créée implémentée. Cette fois, appliquez à l'ensemble de données suivant.

<détails> <résumé> Code source </ résumé>

J'utilise un bloc-notes Jupyter.

`python`


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
%matplotlib inline

np.random.seed(0)
cmap = plt.rcParams['axes.prop_cycle'].by_key()['color']

pi = [0.6, 0.2, 0.2]
mu = np.array([[1, 2], [-2, 0], [3, -2]])
var = [0.8, 0.4, 0.2]
size = 100
ndim = 2

def make_data():
    z = np.random.multinomial(1, pvals=pi)
    z = np.where(z==1)[0][0]
    cov = var[z] * np.eye(2)
    
    return np.random.multivariate_normal(mu[z], cov)

X = np.array([make_data() for _ in range(100)])

fig = plt.figure()
ax = plt.axes()
ax.scatter(X[:, 0], X[:, 1], color='gray')
ax.set_aspect('equal')
plt.show()

Il s'agit de la version non clusterisée de l'ensemble de données que je vous ai montré plus tôt. Le nombre de données est de 100, et il est créé en sélectionnant au hasard parmi 3 distributions gaussiennes. Par conséquent, il convient de le diviser en trois groupes.

Maintenant, regroupons.

`python`


model = IGMMGibbsSampling(var=0.5, var_pri=4, alpha=0.1, seed=0)  #Passer les hyper paramètres au constructeur
clusters = model.predict_clusters(X, n_iter=10)  #Exécution du clustering

#Visualisation
fig = plt.figure()
ax = plt.axes()
for k in range(max(clusters) + 1):
    data = np.array([
        X[i, :] for i in range(X.shape[0])
        if clusters[i] == k
    ])
    ax.scatter(data[:, 0], data[:, 1], label=f'cluster {k + 1}', color=cmap[k])

ax.set_aspect('equal')
plt.legend()
plt.show()

La distribution commune var dans le cluster est un petit ensemble de données de 0,5 $, la distribution de la distribution antérieure moyenne var_pri est de 4 $ pour la faire ressembler à une distribution uniforme, et le paramètre du processus de cuisson chinois ʻalpha` est de 0,1 $. Il est défini sur $. À propos, la moyenne de la distribution antérieure de la moyenne est implémentée pour être la moyenne de l'ensemble de données. Nous transmettons l'ensemble de données et le nombre d'itérations «n_iter» à la méthode de clustering «predict_clusters». Ici, le nombre d'itérations est fixé à 10.

Le résultat est le suivant.

Même si je n'ai pas spécifié le nombre de clusters et que je suis parti du même état de cluster, ils ont été divisés en 3 clusters.

Comme c'est un gros problème, examinons la transition du clustering pour chaque nombre d'itérations. Les résultats de regroupement à «n_iter = 2, 4, 6, 8, 10» sont affichés.

n_iter = 2

n_iter = 4

n_iter = 6

n_iter = 8

n_iter = 10

Il semble qu'un grand nombre de clusters ont été créés à partir de l'état initial, puis le nombre de clusters a progressivement diminué à trois. Au fait, même si j'augmentais n_iter à partir d'ici, le nombre de clusters maintenait fondamentalement 3 états.

Parfois, un cluster avec une seule donnée a été créé et le nombre de clusters est devenu quatre.

Réglage des hyper paramètres

Comme mentionné ci-dessus, plus l'hyper paramètre $ \ alpha $ est grand dans le processus de cuisson chinois, plus il est facile de former de nouveaux clusters. À l'origine, il y avait une raison d'omettre l'estimation du nombre de clusters, et si vous avez du mal à ajuster ce $ \ alpha $, l'importance de son introduction sera diminuée. Quand j'ai essayé en changeant $ \ alpha $, j'ai eu l'impression qu'en gros, il n'y a pas de problème s'il est réduit. La figure suivante montre le résultat de $ \ alpha = 0,0001 $ et 10 itérations.

Il est fermement divisé en 3 clusters. Même si $ \ alpha $ est assez petit, il semble qu'il puisse être clusterisé correctement sans grand impact sur le nombre d'itérations.

À propos, lorsque $ \ alpha = 1 $ a été défini, 5 clusters ont été formés au moment de 10 itérations, et le nombre n'a pas diminué à 3 après cela. Il semble que le nombre de clusters ne puisse pas être bien déterminé si $ \ alpha $ est trop grand. Il a les mêmes propriétés que la largeur de pas de l'algorithme d'optimisation utilisé en apprentissage automatique.

en conclusion

Jusqu'à présent, j'ai essayé de regrouper avec un modèle gaussien infiniment mélangé. C'était intéressant car il était divisé en nombre de grappes comme prévu.

[^ 1]: Un ensemble de $ (\ pi_1, \ cdots, \ pi_K) $ qui satisfait $ \ pi_1, \ cdots, \ pi_K \ ge 0 $ et $ \ sum_ {k = 1} ^ K \ pi_k = 1 $ Est appelée l'unité de dimension $ (K-1) $. Initialement utilisée en géométrie, l'unité de dimension $ 0 $ représente un point, l'unité de dimension $ 1 $ représente un segment de ligne et l'unité de dimension $ 2 $ représente un triangle. Si vous lisez le livre d'estimation bayésienne, vous verrez cet ensemble utilisé dans le contexte de la probabilité. [^ 2]: Il existe aussi une méthode pour estimer $ \ pi $ le plus à partir de l'observation de $ z_ {1: N} ^ {\ backslash i} $, mais cette méthode toujours $ \ lorsque le cluster $ k $ disparaît. Notez que pi_k = 0 $ et le cluster ne sera pas réactivé. [^ 4]: $ N $ variables $ x_1, \ cdots, x_N $ seront abrégés en $ x_ {1: N} $. Aussi, j'écrirai $ x_ {1: N} ^ {\ backslash i} $ pour les variables $ N-1 $ à l'exclusion des $ i $ th. [^ 5]: Même si $ \ Sigma_k $ est une valeur fixe commune aux clusters, je pense que le clustering peut être effectué avec la même précision que k-means. C'est parce que la méthode des k-moyennes équivaut à l'estimation la plus probable de chaque $ \ mu_k $ et $ z_i $ dans le modèle gaussien mixte avec une variance fixe comme décrit ci-dessus.