[Python] J'ai essayé d'implémenter un échantillonnage de Gibbs marginalisé

Symboles utilisés dans cet article

Voici quelques-uns des symboles mathématiques utilisés dans cet article qui sont rarement trouvés ailleurs.

• $ x_ {1: N} $: $ x_1, \ cdots, x_N $
• $ x_ {1: N} ^ {\ backslash i} $: $ x_ {1: N} $ moins $ x_i $

En outre, le vecteur de cet article est un vecteur de colonne et n'apparaît pas en gras.

Examen de Gibbs Sampling

Avant de parler de marginalisation, revenons un peu plus en détail sur l'échantillonnage de Gibbs précédent.

Ce que nous avons fait la dernière fois, c'est de supposer que l'ensemble de données $ x_ {1: N} $ est représenté par un modèle gaussien mixte de clusters $ K $, et que la distribution gaussienne moyenne pour chaque cluster est $ \ mu_ {1: K. } $ Et le cluster $ z_ {1: N} $ de chaque donnée devait être estimé. La distribution simultanée de $ z_ {1: N} $ et $ \ mu_ {1: K} $ donne les probabilités conditionnelles suivantes, qui sont difficiles à gérer analytiquement.

python

P(\mu_{1: K}, z_{1:N} | x_{1: N}, \Sigma_{1: K}, \pi)

Dans l'équation ci-dessus, $ \ Sigma_ {1: K} $ est la matrice de covariance de la distribution gaussienne, et $ \ pi: = (\ pi_1, \ cdots, \ pi_K) ^ T $ est la distribution catégorique que chaque $ z_i $ suit. C'est un paramètre. [^ 2] Par conséquent, au lieu d'estimer $ \ mu_ {1: K}, z_ {1: N} $ en même temps, chaque variable $ \ mu_1, \ cdots, \ mu_K, z_1, \ cdots, z_N $ est estimée séquentiellement une par une. L'ensemble est estimé par échantillonnage et répétition. Cette technique est appelée ** échantillonnage de Gibbs **. Par conséquent, il sera estimé dans les deux étapes suivantes.

--Pour chaque $ k = 1, \ cdots, K $, $ P (\ mu_k | \ mu_ {1: K} ^ {\ backslash k}, z_ {1: N}, x_ {1: N}, \ Échantillon $ \ mu_k $ basé sur Sigma_ {1: K}, \ pi, \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri}) $ --Pour chaque $ i = 1, \ cdots, N $, $ P (z_i | z_ {1: N} ^ {\ backslash i}, \ mu_ {1: K}, x_ {1: N}, \ Sigma_ {1: K}, \ pi) Échantillon $ z_i $ basé sur $

$ \ Mu_ {\ rm pri}, \ Sigma_ {\ rm pri} $ sont des paramètres de la distribution antérieure de $ \ mu_k $. [^ 1]

Qu'est-ce que l'échantillonnage de Gibbs périphérique?

Dans l'échantillonnage de Gibbs ci-dessus, la partie qui estime la grappe de données est la deuxième étape, et la première étape est le processus auxiliaire requis pour calculer la grappe. Ici, en utilisant une technique appelée ** marginalisation **, du coup la deuxième étape "$ z_ {1: N} $" est effectuée sans la première étape "échantillonnage de $ \ mu_ {1: N} $". Pensez à faire un «échantillonnage». Cette technique est appelée ** échantillonnage de Gibbs marginalisé **. L'échantillonnage de Gibbs périphérique répète une étape:

--Pour chaque $ i = 1, \ cdots, N $, $ P (z_i | z_ {1: N} ^ {\ backslash i}, x_ {1: N}, \ Sigma_ {1: K}, \ pi , \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri}) Échantillon $ z_i $ basé sur $

Dans l'échantillonnage de Gibbs marginalisé, même si la variance de $ \ mu_ {1: K} $ est grande [^ 6], $ z_ {1 : N} $ peut être échantillonné. Vous vous demandez peut-être comment estimer sans $ \ mu_ {1: K} $, mais les informations pour $ \ mu_ {1: K} $ sont $ z_ {1: N} ^ {\ backslash i}, x_ {1: N}, \ Sigma_ {1: K}, \ pi, \ mu_ {\ rm pri}, \ Sigma_ {\ rm pri} $ les ont, alors utilisez-les Image pour estimer $ z_i $ (sans utiliser $ \ mu_ {1: K} $). Maintenant, sortons de l'histoire d'échantillonnage de Gibbs et discutons de la marginalisation.

Périphérialisation de la distribution simultanée

Supposons que $ x \ in X, y \ in Y $ est une variable stochastique continue et que la distribution simultanée $ P (x, y) $ est connue. À ce stade, l'élimination de la variable de probabilité à l'aide de l'intégration comme suit s'appelle ** marginalisation **.

python

P(x) = \int_{Y}P(x, y)dy

Transformer le côté droit en utilisant le théorème de Bayes et l'écrire avec une probabilité conditionnelle donne: Utilisez cette option si vous connaissez la probabilité conditionnelle plutôt que la distribution simultanée.

python

P(x) = \int_{Y}P(x|y)P(y)dy

En passant, lorsque la variable de probabilité $ y $ est discrète, la formulation de marginalisation est la suivante. Cela peut être plus facile à comprendre.

python

P(x) = \sum_{y \in Y}P(x, y) = \sum_{y \in Y}P(x|y)P(y)

Application à l'échantillonnage Gibbs

Revenons à l'histoire de l'échantillonnage Gibbs. Effectuez une marginalisation pour trouver la distribution de probabilité de la version $ z_i $ none de $ \ mu_ {1: K} $. $ D $ est la dimension de $ x_i $ et $ \ mu_k $. Puisqu'il existe de nombreuses variables conditionnelles de probabilité conditionnelle, elles sont omises.

python

P(z_i=k | \cdots)
= \int_{\mathbb{R}^D} P(z_i=k | \mu_{1:K}, \cdots)P(\mu_{1:K} | \cdots)d\mu_{1:K}

Comme précédemment, en supposant $ \ Sigma_k = \ sigma_k ^ 2 I_D, \ Sigma_ {\ rm pri}: = \ sigma_ {\ rm pri} ^ 2I_D $, ceci est calculé comme suit. [^ 3]

python

\begin{align}
P(z_i=k | \cdots)
&\propto \pi_k N\left(x_i | a_k\left(\frac{n_k^{\backslash i}}{\sigma^2}\overline{x_k}^{\backslash i} + \frac{1}{\sigma_{\rm pri}^2} \mu_{\rm pri}  \right),
 \left( \sigma^2 + a_k \right)I_D\right) \\
a_k &:= \left( \frac{n_k^{\backslash i}}{\sigma^2} + \frac{1}{\sigma_{\rm pri}^2} \right)^{-1}
\end{align}

$ n_k ^ {\ backslash i} $ est le nombre de données appartenant au cluster $ k $ dans l'ensemble de données obtenu en soustrayant $ x_i $ de $ x_ {1: N} $, $ \ override {x_k} ^ {\ la barre oblique inverse i} $ est leur moyenne.

la mise en oeuvre

Sur la base de l'implémentation précédente, nous implémenterons un échantillonnage de Gibbs marginalisé. Il y a deux changements majeurs par rapport à la dernière fois.

--Supprimer la variable de classe mu (inutile car elle n'est pas échantillonnée)

• Changement de la méthode log_deformed_gaussian

Le premier est l'idée centrale de l'échantillonnage de Gibbs périphérique et est le plus important. Laissez-moi vous expliquer le deuxième. La dernière fois, afin d'obtenir le rapport des valeurs de la fonction de densité de probabilité de la distribution normale, le calcul a été effectué en ignorant les facteurs qui ne dépendent pas de la grappe «k». Je republierai la formule précédente. [^ 4]

python

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

Ici, comme la formule- \frac{1}{2}\log|\Sigma_k|J'ai pu ignorer le facteur de\Sigma_kEst un clusterkC'est parce que c'était une valeur fixe qui ne dépendait pas de. Cependant, cette fois, $ \ Sigma_k = (\ sigma ^ 2 + a_k) I_D $, et $ a_k $ change avec $ k $, donc ce facteur ne peut pas être ignoré. Par conséquent, la formule pour ce temps est la suivante.

python

\begin{align}
P(z_i = k | x_{1:N}^{\backslash i}, z_{1:N}^{\backslash i},  \mu_{1:K}, \Sigma_{1:K}, \pi) 
&\propto \pi_k \exp\left\{ \frac{D}{2}\log \sigma_k^2 -\frac{1}{2 \sigma^2}\ \| x_i - \mu_k \|_2^2\right\}
\end{align}

L'implémentation de la méthode log_deformed_gaussian qui calcule le contenu de cette ʻexp` est la suivante.

python

def log_deformed_gaussian(x, mu, var):
    D = x.shape[0]
    norm_squared = ((x - mu) * (x - mu)).sum(axis=1)
    return -D * np.log(var) / 2 - norm_squared / (2 * var)

Sur la base de ce qui précède, j'ai essayé de le mettre en œuvre. C'est long donc c'est plié.

<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 -np.log(var) / 2 - 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]


#Classe de regroupement avec échantillonnage de Gibbs marginalisé
class GMMClustering(object):
    def __init__(self, X, K, var=1, var_pri=1, seed=None):
        self.X = X.copy()
        self.K = K
        self.N = X.shape[0]
        self.D = X.shape[1]
        self.z = None

        #Paramétrage de la distribution de probabilité
        # self.mu n'est pas utilisé, donc il n'est pas défini
        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, n_iter):
        init_z = self._random.randint(0, self.K, self.N)
        self.z = ClusterArray(init_z)

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

    def _sample_zi(self, i):
        n_vec = np.array([
            self.z.count(k) - bool(k == self.z.array[i])
            for k in range(self.K)
        ])
        x_bar_vec = np.array([self._mean_in_removed(k, i) for k in range(self.K)])

        tmp = 1/(n_vec/self.var + 1/self.var_pri)
        mu = np.tile(tmp, (self.D, 1)).T * (np.tile(n_vec, (self.D, 1)).T * x_bar_vec / self.var + self.mu_pri / self.var_pri)
        var = tmp + self.var

        log_probs_xi = log_deformed_gaussian(self.X[i], mu, 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 _mean_in_removed(self, k, i):
        # self.Calculer la moyenne des données appartenant au cluster k à partir de l'ensemble de données avec le i-ème soustrait de X
        mean = np.array([
            x for j, x in enumerate(self.X)
            if (j != i) and (self.z[j] == k)
        ]).mean(axis=0)
        return mean

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(X, K=3, var=0.05, seed=1)
gmc.fit(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'
)