introduction

J'ai implémenté une distribution normale mixte en python. J'ai utilisé ["Première reconnaissance de formes"](https://www.amazon.co.jp/First reconnaissance de formes-Hirai-Yuzo / dp / 4627849710) comme manuel.

Structure de cet article

Introduction
Distribution normale mixte
Modèle de distribution normale mixte
Variables cachées et probabilités postérieures
Log de vraisemblance et fonction Q
Estimation des paramètres par algorithme EM
Implémentation en python
Résultat
en conclusion

Distribution normale mixte

En appliquant un modèle probabiliste à la distribution des données, il est possible de déterminer de manière probabiliste à quel cluster chaque donnée appartient. Étant donné que de nombreux modèles probabilistes ne peuvent représenter que des distributions probabilistes monomodales, il est nécessaire de modéliser la distribution probabiliste globale avec une somme linéaire pondérée de plusieurs modèles probabilistes. Soit le nombre de clusters $ K $ et le modèle de probabilité du $ k $ th cluster $ p_k (\ boldsymbol x) $, et la distribution de probabilité globale est exprimée comme suit.

p(\boldsymbol x) = \sum_{k = 1}^K \pi_kp_k(\boldsymbol x)

$ \ Pi_k $ est le poids du modèle probabiliste $ k $ th. Un modèle qui utilise une distribution normale comme modèle probabiliste est appelé un ** modèle de distribution normale mixte **.

Modèle de distribution normale mixte

La fonction de distribution normale dimensionnelle $ d $ qui représente le cluster $ k $ ème est exprimée comme suit.

N(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k) = \frac{1}{(2\pi)^{1/d} \ |\boldsymbol \Sigma_k|^{1/2}}\exp\bigl(-\frac{1}{2}(\boldsymbol x - \boldsymbol \mu_k)^T \boldsymbol \Sigma_k^{-1}(\boldsymbol x - \boldsymbol \mu_k)\bigr)

La distribution globale est exprimée comme la somme linéaire de ceux-ci.

p(\boldsymbol x) = \sum_k^K \pi_k N(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k), \ 0 \leq \pi_k \leq 1, \ \sum_{k = 1}^K \pi_k = 1

$ \ Pi_k $ est appelé le rapport de mélange, et les paramètres du modèle de distribution normale mixte sont les suivants.

\boldsymbol \pi = (\pi_1, ..., \pi_K), \ \boldsymbol \mu = (\boldsymbol \mu_1, ..., \boldsymbol \mu_K), \ \boldsymbol \Sigma = (\boldsymbol \Sigma_1, ..., \boldsymbol \Sigma_K)

Vous aurez besoin de $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $ pour bien représenter l'ensemble des données. Cependant, comme on ne sait pas à quel cluster chaque donnée appartient, il n'est pas possible d'obtenir chaque paramètre directement.

Variables cachées et probabilités postérieures

Dans le modèle de distribution normale mixte, afin d'estimer $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $ à partir de l'ensemble des données, Introduisez une variable cachée $ \ boldsymbol z = (z_1, ..., z_K) ^ T $ de la dimension $ K $ qui représente à quel cluster chaque donnée $ \ boldsymbol x $ appartient. $ z_k $ prend $ 1 $ si les données appartiennent au $ k $ ème cluster, $ 0 $ si ce n'est pas le cas.

\sum_{k = 1}^K z_k = 1, \ \boldsymbol z = (0, ..., 0, 1, 0, ..., 0)^T

La distribution simultanée du $ \ boldsymbol x $ observé et de la variable cachée $ \ boldsymbol z $ $ p (\ boldsymbol x, \ boldsymbol z) $ est la suivante.

p(\boldsymbol x, \boldsymbol z) = p(\boldsymbol z)p(\boldsymbol x \ | \ \boldsymbol z)

La distribution des variables cachées $ p (\ boldsymbol z) $ et la distribution conditionnelle $ p (\ boldsymbol x \ | \ \ boldsymbol z) $ dues aux variables cachées des données observées peuvent être exprimées comme suit.

\begin{align}
& p(\boldsymbol z) = \prod_{k = 1}^K \pi_k^{z_k} \\
& p(\boldsymbol x \ | \ \boldsymbol z) = \prod_{k = 1}^K \big[N(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)\bigr]^{z_k}
\end{align}

De plus, $ p (\ boldsymbol x) $ peut être obtenu en ajoutant la distribution simultanée $ p (\ boldsymbol x, \ boldsymbol z) $ pour tout $ \ boldsymbol z $.

\begin{align}
p(\boldsymbol x) &= \sum_{\substack{k = 1 \\ (z_k = 1, z_{j \neq k} = 0)}}^K p(\boldsymbol x, \boldsymbol z) \\
&= \sum_{\substack{k = 1 \\ (z_k = 1, z_{j \neq k} = 0)}}^K p(\boldsymbol z)p(\boldsymbol x \ | \ \boldsymbol z) \\
&= \sum_{\substack{k = 1 \\ (z_k = 1, z_{j \neq k} = 0)}}^K \prod_{k = 1}^K \pi_k^{z_k} \prod_{k = 1}^K \big[N(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)\bigr]^{z_k} \\
&= \sum_{k = 1}^K \pi_kN(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)
\end{align}

En utilisant la distribution ci-dessus, la probabilité postérieure des variables cachées $ \ gamma (z_k) $ est exprimée comme suit.

\begin{align}
\gamma(z_k) &\equiv p(z_k = 1 \ | \ \boldsymbol x) \\
&\\
&= \frac{p(z_k = 1)p(\boldsymbol x \ | \ z_k = 1)}{\sum_{j = 1}^K p(z_j = 1)p(\boldsymbol x \ | \ z_j = 1)} \\
&\\
&= \frac{\pi_kN(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)}{\sum_{j = 1}^K \pi_jN(\boldsymbol x \ | \ \boldsymbol \mu_j, \boldsymbol \Sigma_j)} \tag{*}
\end{align}

Log de vraisemblance et fonction Q

Les données observées $ \ boldsymbol X $ et la variable cachée $ \ boldsymbol Z $ sont représentées comme suit. $ C_k $ représente le cluster $ k $.

\begin{align}
& \boldsymbol X = (\boldsymbol x_1, ..., \boldsymbol x_N), \ \boldsymbol x_i = (x_{i1}, ..., x_{id})^T \\
& \boldsymbol Z = (\boldsymbol z_1, ..., \boldsymbol z_N), \ \boldsymbol z_i = (z_{i1}, ..., z_{iK})^T \\
& z_{ik} = \begin{cases}
1 \ (x_i \in C_k) \\
0 \ (x_i \notin C_k)
\end{cases}
\end{align}

L'ensemble des données observées et des variables cachées est appelé données complètes.

\boldsymbol Y = (\boldsymbol x_1, ..., \boldsymbol x_N, \boldsymbol z_1, ..., \boldsymbol z_N) = (\boldsymbol X, \boldsymbol Z)

Les paramètres de la distribution normale mixte, $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $, sont les paramètres qui maximisent la probabilité de données complètes. La probabilité d'obtenir des données complètes peut être exprimée comme suit.

\begin{align}
p(\boldsymbol Y \ | \ \boldsymbol \pi, \ \boldsymbol \mu, \ \boldsymbol \Sigma) &= p(\boldsymbol Z \ | \ \boldsymbol \pi, \ \boldsymbol \mu, \ \boldsymbol \Sigma) p(\boldsymbol X \ | \ \boldsymbol Z, \ \boldsymbol \pi, \ \boldsymbol \mu, \ \boldsymbol \Sigma) \\
&\\
&= \prod_{i = 1}^N \prod_{k = 1}^K \bigl[\pi_kN(\boldsymbol x_i \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k) \bigr]^{z_{ik}}
\end{align}

Transformez-vous en probabilité logarithmique pour trouver l'estimation la plus probable.

\begin{align}
\tilde{L} &= \ln p(\boldsymbol Y \ | \ \boldsymbol \pi, \ \boldsymbol \mu, \ \boldsymbol \Sigma) \\
&\\
&= \sum_{i = 1}^N \sum_{k = 1}^K z_{ik} \ln \pi_k + \sum_{i = 1}^N \sum_{k = 1}^K z_{ik} \ln N(\boldsymbol x_i \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k) \\
&\\
&= \sum_{i = 1}^N \sum_{k = 1}^K z_{ik} \ln \pi_k + \sum_{i = 1}^N \sum_{k = 1}^K z_{ik} \bigl(-\frac{d}{2}\ln (2\pi) + \frac{1}{2} \ln \bigl|\boldsymbol \Sigma_k \bigr|^{-1} - \frac{1}{2}(\boldsymbol x_i - \boldsymbol \mu_k)^T \boldsymbol \Sigma_k^{-1}(\boldsymbol x_i - \boldsymbol \mu_k) \bigr)
\end{align}

Prend la valeur attendue de la variable cachée de vraisemblance logarithmique.

\begin{align}
L &= E_z \bigl\{\tilde{L} \bigr\} \\
&\\
&= \sum_{i = 1}^N \sum_{k = 1}^K E_{z_{ik}} \bigl\{z_{ik} \bigr\} \ln \pi_k + \sum_{i = 1}^N \sum_{k = 1}^K E_{z_{ik}} \bigl\{z_{ik} \bigr\} \bigl(-\frac{d}{2}\ln (2\pi) + \frac{1}{2} \ln \bigl|\boldsymbol \Sigma_k \bigr|^{-1} - \frac{1}{2}(\boldsymbol x_i - \boldsymbol \mu_k)^T \boldsymbol \Sigma_k^{-1}(\boldsymbol x_i - \boldsymbol \mu_k) \bigr)
\end{align}

La valeur attendue de la variable masquée est la suivante à partir de l'équation (*), qui est la probabilité postérieure de la variable masquée.

\begin{align}
E_{z_{ik}} \bigl\{z_{ik} \bigr\} &= \sum_{z_{ik} = \{0, 1\}} z_{ik}p(z_{ik} \ | \ \boldsymbol x_i, \pi_k, \boldsymbol \mu_k, \boldsymbol \Sigma_k) \\
&\\
&= 1 \times p(z_{ik} = 1 \ | \ \boldsymbol x_i) \\
&\\
& = \frac{p(z_{ik} = 1)p(\boldsymbol x_i \ | \ z_{ik} = 1)}{\sum_{j = 1}^K p(z_{ij} = 1)p(\boldsymbol x \ | \ z_{ij} = 1)} \\
&\\
&= \frac{\pi_kN(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)}{\sum_{j = 1}^K \pi_jN(\boldsymbol x \ | \ \boldsymbol \mu_j, \boldsymbol \Sigma_j)} \\
&\\
&= \gamma(z_{ik})
\end{align}

La fonction qui remplace la valeur attendue de la variable cachée de $ L $ par la probabilité postérieure est appelée ** fonction Q **.

Q = \sum_{i = 1}^N \sum_{k = 1}^K \gamma(z_{ik}) \ln \pi_k + \sum_{i = 1}^N \sum_{k = 1}^K \gamma(z_{ik}) \ln \pi_k \bigl(-\frac{d}{2}\ln (2\pi) + \frac{1}{2} \ln \bigl|\boldsymbol \Sigma_k \bigr|^{-1} - \frac{1}{2}(\boldsymbol x_i - \boldsymbol \mu_k)^T \boldsymbol \Sigma_k^{-1}(\boldsymbol x_i - \boldsymbol \mu_k) \bigr)

Estimation des paramètres par algorithme EM

L '** algorithme EM ** est utilisé comme méthode pour trouver l'estimation la plus probable des paramètres d'un modèle de probabilité avec des variables cachées. L'algorithme EM comprend deux étapes.

Étape d'attente: opération qui prend une valeur attendue pour une variable cachée. Les paramètres du modèle sont fixes et la probabilité postérieure des variables cachées est calculée.
Étape de maximisation: trouvez le paramètre de modèle qui maximise la fonction Q en utilisant la probabilité postérieure de la variable cachée obtenue à l'étape E.

L'algorithme EM répète les étapes E et M jusqu'à ce qu'il n'y ait aucun changement dans la probabilité logarithmique des données complètes. Puisque la valeur de convergence dépend de la valeur initiale, seule la solution optimale locale peut être obtenue. Il est nécessaire de changer la valeur initiale et d'utiliser l'algorithme EM plusieurs fois pour obtenir la valeur estimée la plus probable et de sélectionner la meilleure parmi celles-ci. L'algorithme EM peut être résumé comme suit.

Initialisation de $ \ pi_k, \ \ boldsymbol \ mu_k, \ \ boldsymbol \ Sigma_k $
Étape E: Estimer la probabilité postérieure $ \ gamma (z_ {ik}) $ en utilisant les paramètres courants $ \ pi_k, \ \ boldsymbol \ mu_k, \ \ boldsymbol \ Sigma_k $ $\gamma(z_{ik}) = \frac{\pi_kN(\boldsymbol x \ | \ \boldsymbol \mu_k, \boldsymbol \Sigma_k)}{\sum_{j = 1}^K \pi_jN(\boldsymbol x \ | \ \boldsymbol \mu_j, \boldsymbol \Sigma_j)}$
Étape M: ré-estimer les paramètres du modèle en utilisant l'estimation $ \ gamma (z_ {ik}) $ $ \begin{align} & N_k = \sum_{i = 1}^N \gamma(z_{ik}) \\\ & \\\ & \boldsymbol \mu_k = \frac{1}{N_k} \sum_{i = 1}^N \gamma(z_{ik}) \boldsymbol x_i \\\ & \\\ & \boldsymbol \Sigma_k = \frac{1}{N_k} \sum_{i = 1}^N \gamma(z_{ik}) (\boldsymbol x_i - \boldsymbol \mu_k)(\boldsymbol x_i - \boldsymbol \mu_k)^T \\\ & \\\ & \pi_k = \frac{N_k}{N} \end{align} $
Jugement final: lorsque la probabilité logarithmique de données complètes n'a pas convergé, recalculez à partir de 2. S'il a convergé, il se termine.

Implémentation en python

Je l'ai implémenté comme suit.

`mixturesofgaussian.py`


import numpy
from matplotlib import pyplot
import sys

def create_data(N, K):
    X, mu_star, sigma_star = [], [], []
    for i in range(K):
        loc = (numpy.random.rand() - 0.5) * 10.0 # range: -5.0 - 5.0
        scale = numpy.random.rand() * 3.0 # range: 0.0 - 3.0
        X = numpy.append(X, numpy.random.normal(loc = loc, scale = scale, size = N / K))
        mu_star = numpy.append(mu_star, loc)
        sigma_star = numpy.append(sigma_star, scale)
    return (X, mu_star, sigma_star)

def gaussian(mu, sigma):
    def f(x):
        return numpy.exp(-0.5 * (x - mu) ** 2 / sigma) / numpy.sqrt(2 * numpy.pi * sigma)
    return f

def estimate_posterior_likelihood(X, pi, gf):
    l = numpy.zeros((X.size, pi.size))
    for (i, x) in enumerate(X):
        l[i, :] = gf(x)
    return pi * l * numpy.vectorize(lambda y: 1 / y)(l.sum(axis = 1).reshape(-1, 1))

def estimate_gmm_parameter(X, gamma):
    N = gamma.sum(axis = 0)
    mu = (gamma * X.reshape((-1, 1))).sum(axis = 0) / N
    sigma = (gamma * (X.reshape(-1, 1) - mu) ** 2).sum(axis = 0) / N
    pi = N / X.size
    return (mu, sigma, pi)

def calc_Q(X, mu, sigma, pi, gamma):
    Q = (gamma * (numpy.log(pi * (2 * numpy.pi * sigma) ** (-0.5)))).sum()
    for (i, x) in enumerate(X):
        Q += (gamma[i, :] * (-0.5 * (x - mu) ** 2 / sigma)).sum()
    return Q

if __name__ == '__main__':

    # data
    K = 2
    N = 1000 * K
    X, mu_star, sigma_star = create_data(N, K)

    # termination condition
    epsilon = 0.000001

    # initialize gmm parameter
    pi = numpy.random.rand(K)
    mu = numpy.random.randn(K)
    sigma = numpy.abs(numpy.random.randn(K))
    Q = -sys.float_info.max
    delta = None

    # EM algorithm
    while delta == None or delta >= epsilon:
        gf = gaussian(mu, sigma)

        # E step: estimate posterior probability of hidden variable gamma
        gamma = estimate_posterior_likelihood(X, pi, gf)

        # M step: miximize Q function by estimating mu, sigma and pi
        mu, sigma, pi = estimate_gmm_parameter(X, gamma)

        # calculate Q function
        Q_new = calc_Q(X, mu, sigma, pi, gamma)
        delta = Q_new - Q
        Q = Q_new

    # result
    print u'mu*: %s, simga*: %s' % (str(numpy.sort(numpy.around(mu_star, 3))), str(numpy.sort(numpy.around(sigma_star, 3))))
    print u'mu : %s, sgima : %s' % (str(numpy.sort(numpy.around(mu, 3))), str(numpy.sort(numpy.around(sigma, 3))))

    # plot
    n, bins, _ = pyplot.hist(X, 50, normed = 1, alpha = 0.3)
    seq = numpy.arange(-15, 15, 0.02)
    for i in range(K):
        pyplot.plot(seq, gaussian(mu[i], sigma[i])(seq), linewidth = 2.0)
    pyplot.savefig('gmm_graph.png')
    pyplot.show()

résultat

Les résultats suivants ont été obtenus. Il semble qu'une estimation raisonnable ait été faite.

en conclusion

J'ai pu implémenter une distribution normale mixte en python. Puisqu'il y a certaines parties que je ne comprends pas, comme la méthode de dérivation de chaque paramètre dans l'étape M et la nature de l'algorithme EM. Je vais l'étudier et l'ajouter à l'avenir.

Implémentation de distribution normale mixte en python