Distribution gaussienne mixte et logsumexp

Aperçu

La distribution gaussienne mixte est un modèle d'addition pondérée de multiples distributions gaussiennes et peut représenter une distribution multimodale. Vous pouvez également utiliser l'algorithme EM pour optimiser les paramètres de ce modèle. La valeur attendue de la variable latente est calculée par l'étape E de cet algorithme EM, mais si le calcul est effectué naïf à ce moment-là, des problèmes de débordement et de sous-débordement peuvent survenir. À ce moment-là, ce problème peut être évité en utilisant une célèbre méthode de calcul numérique appelée logsumexp.

Vue d'ensemble de la distribution gaussienne mixte et de logsumexp

Concernant l'optimisation des paramètres par l'algorithme EM de la distribution gaussienne mixte, seule la partie liée à logsumexp est brièvement représentée. Pour plus de détails, ["Reconnaissance de formes et apprentissage automatique"](http://www.amazon.co.jp/%E3%83%91%E3%82%BF%E3%83%BC%E3%83%B3% E8% AA% 8D% E8% AD% 98% E3% 81% A8% E6% A9% 9F% E6% A2% B0% E5% AD% A6% E7% BF% 92-% E4% B8% 8B-% E3% 83% 99% E3% 82% A4% E3% 82% BA% E7% 90% 86% E8% AB% 96% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E7% B5% B1% E8% A8% 88% E7% 9A% 84% E4% BA% 88% E6% B8% AC-CM-% E3% 83% 93% E3% 82% B7% E3% 83% A7% E3 % 83% 83% E3% 83% 97 / dp / 4621061240 / ref = pd_sim_14_1? Ie = UTF8 & refRID = 1GS4SDYG8VEYR40YSADB) Veuillez vous référer à. (Pour le moment, je vous donnerai le matériel lorsque je donnerai une conférence à SlideShare PRML 9.0-9.2)

\begin{aligned}
p(x|\theta) &= \sum_k \pi_k N(x_n|\mu_k,\Sigma_k)\\
\pi_k&:taux de mélange
\end{aligned}

La distribution de la variable observée $ x $ de la distribution gaussienne mixte est exprimée par l'équation ci-dessus. Afin d'appliquer l'algorithme EM, nous introduisons la variable latente $ z $ et exprimons la distribution simultanée de $ x et z $ comme suit.

\begin{aligned}
p(x,z) &= \prod_k \pi_k^{z_k} N(x_n|\mu_k,\Sigma_k)^{z_k}\\
z&: 1ofK
\end{aligned}

En utilisant ces derniers, la valeur attendue de la variable latente z peut être calculée à partir du théorème de Bayes comme suit.

\begin{aligned}
\gamma(z_{nk})
&= E[z_{nk}] \\ 
&=\frac{\pi_k N(x_n|\mu_k,\Sigma_k)}{\sum_{k^\prime}\pi_{k^\prime} N(x_n|\mu_{k^\prime},\Sigma_{k^\prime})}
\end{aligned}

En E-step de l'optimisation par l'algorithme EM gaussien mixte, la valeur attendue de cette variable latente est calculée. En regardant le dénominateur, nous pouvons voir que le sous-débit peut se produire dans le calcul exponentiel de la distribution gaussienne même lorsque l'opération est effectuée à l'échelle logarithmique car il y a une opération de somme de la distribution gaussienne. Par conséquent, il est nécessaire de concevoir quelque chose, et une méthode est logsumexp.

logsumexp

\log(\sum^N_{i=1} \exp(x_i))

Même si le résultat de ce calcul lui-même ne déborde pas ou ne sous-est pas, l'individu $ \ exp (x_i) $ peut déborder ou sous-déborder. Par conséquent, la formule est transformée comme suit.

\begin{aligned}
\log(\sum^N_{i=1} \exp(x_i)) 
&= \log\{\exp(x_{max})\sum^N_{i=1} \exp(x_i - x_{max})\} \\
& = \log\{\sum^N_{i=1} \exp(x_i - x_{max})\}  + x_{max}
\end{aligned}

En effectuant le calcul de cette manière, l'argument de l'opération exponentielle se situe dans la plage $ [x_ {min} --x_ {max}, 0] $, donc la possibilité de débordement et de sous-dépassement est considérablement réduite.

Ce logsumexp est également implémenté dans la bibliothèque d'apprentissage automatique de Python scikitlearn et est facile à utiliser (Scikit Learn Utilities for Developers. )). De plus, l'implémentation de logsumexp dans scikit learn est sklearn / utils / extmath.py C'est dans ce qui suit. Vous pouvez voir que c'est exactement la même chose que la formule ci-dessus, sauf pour la partie correspondant à l'axe.


def logsumexp(arr, axis=0):
    """Computes the sum of arr assuming arr is in the log domain.
    Returns log(sum(exp(arr))) while minimizing the possibility of
    over/underflow.
    Examples
    --------
    >>> import numpy as np
    >>> from sklearn.utils.extmath import logsumexp
    >>> a = np.arange(10)
    >>> np.log(np.sum(np.exp(a)))
    9.4586297444267107
    >>> logsumexp(a)
    9.4586297444267107
    """
    arr = np.rollaxis(arr, axis)
    # Use the max to normalize, as with the log this is what accumulates
    # the less errors
    vmax = arr.max(axis=0)
    out = np.log(np.sum(np.exp(arr - vmax), axis=0))
    out += vmax
    return out

Implémentation du modèle gaussien mixte

Sur la base de ce qui précède, j'ai implémenté le modèle gaussien mixte, donc je vais le laisser (https://github.com/seataK/machine-learning/blob/master/GaussianMixture/GaussianMixture.py).

import numpy as np
import random
import pylab as plt
from sklearn.utils import extmath
from sklearn.cluster import KMeans
import sys


na = np.newaxis


class DataFormatter:
    def __init__(self, X):
        self.mean = np.mean(X, axis=0)
        self.std = np.std(X, axis=0)

    def standarize(self, X):
        return (X - self.mean[na, :]) / self.std[na, :]


def log_gaussian(X, mu, cov):
    d = X.shape[1]
    det_sig = np.linalg.det(cov)
    A = 1.0 / (2*np.pi)**(d/2.0) * 1.0 / det_sig**(0.5)
    x_mu = X - mu[na, :]
    inv_cov = np.linalg.inv(cov)
    ex = - 0.5 * np.sum(x_mu[:, :, na] * inv_cov[na, :, :] *
                        x_mu[:, na, :], axis=(1, 2))
    return np.log(A) + ex


class GMM:
    def __init__(self,
                 K=2,
                 max_iter=300,
                 diag=False):
        self.K = K
        self.data_form = None
        self.pi = None
        self.mean = None
        self.cov = None
        self.max_iter = max_iter
        self.diag = diag

    def fit(self, _X):
        self.data_form = DataFormatter(_X)
        X = self.data_form.standarize(_X)
        N = X.shape[0]
        D = X.shape[1]
        K = self.K

        # init parameters using K-means
        kmeans = KMeans(n_clusters=self.K)

        kmeans.fit(X)

        self.mean = kmeans.cluster_centers_

        self.cov = np.array([[[1 if i == j else 0
                             for i in range(D)]
                             for j in range(D)]
                             for k in range(K)])

        self.pi = np.ones(K) / K

        # Optimization
        for _ in range(self.max_iter):
            # E-step

            gam_nk = self._gam(X)

            # M-step
            Nk = np.sum(gam_nk, axis=0)

            self.pi = Nk / N

            self.mean = np.sum(gam_nk[:, :, na] * X[:, na, :],
                               axis=0) / Nk[:, na]

            x_mu_nkd = X[:, na, :] - self.mean[na, :, :]

            self.cov = np.sum(gam_nk[:, :, na, na] *
                              x_mu_nkd[:, :, :, na] *
                              x_mu_nkd[:, :, na, :],
                              axis=0) / Nk[:, na, na]

            if(self.diag):
                for k in range(K):
                    var = np.diag(self.cov[k])
                    self.cov[k] = np.array([[var[i] if i == j else 0
                                             for i in range(D)]
                                            for j in range(D)])

    def _gam(self, X):
        log_gs_nk = np.array([log_gaussian(X, self.mean[i], self.cov[i])
                              for i in range(self.K)]).T

        log_pi_gs_nk = np.log(self.pi)[na, :] + log_gs_nk

        log_gam_nk = log_pi_gs_nk[:, :] - extmath.logsumexp(log_pi_gs_nk, axis=1)[:, na]

        return np.exp(log_gam_nk)

    def predict(self, _X):
        X = self.data_form.standarize(_X)

        gam_nk = self._gam(X)

        return np.argmax(gam_nk, axis=1)

Autre

Dans logsumexp, de nombreuses opérations de log et exp apparaissent, donc la vitesse de calcul semble lente (logsumexp est l'histoire noire de l'humanité). Si la vitesse devient un goulot d'étranglement, il peut être préférable d'utiliser une autre méthode.

Recommended Posts

Distribution gaussienne mixte et logsumexp
EM de distribution gaussienne mixte
Statistiques de base et distribution gaussienne
PRML Chapitre 9 Implémentation Python de distribution gaussienne mixte
Estimation de la distribution gaussienne mixte par la méthode de Bayes
Estimation de la distribution gaussienne mixte par algorithme EM
Distribution et test
Calcul d'algorithme EM pour un problème de distribution gaussienne mixte
Distribution OS et Linux
Implémentation Python Distribution mixte de Bernoulli
Test d'hypothèse et distribution de probabilité
Tri avec un mélange de chiffres et de lettres
Implémentation de distribution normale mixte en python
Un mélange d'arguments de position, de tapotement et de dictionnaire