Einführung

Ich habe eine gemischte Normalverteilung in Python implementiert. Ich habe ["Erste Mustererkennung"](https://www.amazon.co.jp/First Pattern Recognition-Hirai-Yuzo / dp / 4627849710) als Lehrbuch verwendet.

Struktur dieses Artikels

Einführung
Gemischte Normalverteilung
Gemischtes Normalverteilungsmodell
Versteckte Variablen und hintere Wahrscheinlichkeiten
Log Likelihood und Q Funktion
Parameterschätzung durch EM-Algorithmus
Implementierung in Python
Ergebnis
abschließend

Gemischte Normalverteilung

Durch Anwenden eines Wahrscheinlichkeitsmodells auf die Datenverteilung ist es möglich, probabilistisch zu bestimmen, zu welchem Cluster die einzelnen Daten gehören. Da viele Wahrscheinlichkeitsmodelle nur monomodale Wahrscheinlichkeitsverteilungen darstellen können, ist es erforderlich, die gesamte Wahrscheinlichkeitsverteilung mit einer gewichteten linearen Summe mehrerer Wahrscheinlichkeitsmodelle zu modellieren. Die Anzahl der Cluster sei $ K $ und das Wahrscheinlichkeitsmodell des $ k $ -ten Clusters sei $ p_k (\ boldsymbol x) $, und die Gesamtwahrscheinlichkeitsverteilung wird wie folgt ausgedrückt.

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

$ \ Pi_k $ ist das Gewicht des $ k $ -ten Wahrscheinlichkeitsmodells. Ein Modell, das eine Normalverteilung als Wahrscheinlichkeitsmodell verwendet, wird als ** gemischtes Normalverteilungsmodell ** bezeichnet.

Gemischtes Normalverteilungsmodell

Die $ d $ dimensionale Normalverteilungsfunktion, die den $ k $ -ten Cluster darstellt, wird wie folgt ausgedrückt.

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)

Die Gesamtverteilung wird als lineare Summe davon ausgedrückt.

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 $ wird als Mischungsverhältnis bezeichnet, und die Parameter des gemischten Normalverteilungsmodells sind wie folgt.

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

Sie benötigen $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $, um die gesamten Daten gut darzustellen. Da jedoch nicht bekannt ist, zu welchem Cluster die einzelnen Daten gehören, ist es nicht möglich, jeden Parameter direkt abzurufen.

Versteckte Variablen und hintere Wahrscheinlichkeiten

Um im gemischten Normalverteilungsmodell $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $ aus den gesamten Daten zu schätzen, Führen Sie eine $ K $ dimensionale versteckte Variable $ \ boldsymbol z = (z_1, ..., z_K) ^ T $ ein, die darstellt, zu welchem Cluster die einzelnen Daten $ \ boldsymbol x $ gehören. $ z_k $ benötigt $ 1 $, wenn die Daten zum $ k $ -ten Cluster gehören, $ 0 $, wenn dies nicht der Fall ist.

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

Die gleichzeitige Verteilung des beobachteten $ \ boldsymbol x $ und der versteckten Variablen $ \ boldsymbol z $ $ p (\ boldsymbol x, \ boldsymbol z) $ ist wie folgt.

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

Die versteckte Variablenverteilung $ p (\ boldsymbol z) $ und die bedingte Verteilung $ p (\ boldsymbol x \ | \ \ boldsymbol z) $ aufgrund der versteckten Variablen der beobachteten Daten können wie folgt ausgedrückt werden.

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

Darüber hinaus kann $ p (\ boldsymbol x) $ erhalten werden, indem die gleichzeitige Verteilung $ p (\ boldsymbol x, \ boldsymbol z) $ für alle $ \ boldsymbol z $ addiert wird.

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

Unter Verwendung der obigen Verteilung wird die hintere Wahrscheinlichkeit versteckter Variablen $ \ gamma (z_k) $ wie folgt ausgedrückt.

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

Protokollwahrscheinlichkeit und Q-Funktion

Die beobachteten Daten $ \ boldsymbol X $ und die versteckte Variable $ \ boldsymbol Z $ werden wie folgt dargestellt. $ C_k $ repräsentiert den 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}

Der Satz beobachteter Daten und versteckter Variablen wird als vollständige Daten bezeichnet.

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

Die Parameter der gemischten Normalverteilung $ \ boldsymbol \ pi, \ \ boldsymbol \ mu, \ \ boldsymbol \ Sigma $ sind die Parameter, die die Wahrscheinlichkeit vollständiger Daten maximieren. Die Wahrscheinlichkeit vollständiger Daten kann wie folgt ausgedrückt werden.

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

In logarithmische Wahrscheinlichkeit umwandeln, um die wahrscheinlichste Schätzung zu finden.

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

Nimmt den erwarteten Wert für die versteckte Variable der logarithmischen Wahrscheinlichkeit.

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

Der erwartete Wert für die versteckte Variable ergibt sich aus der Gleichung (*), die die hintere Wahrscheinlichkeit der versteckten Variablen darstellt.

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

Die Funktion, die den erwarteten Wert der versteckten Variablen von $ L $ durch die hintere Wahrscheinlichkeit ersetzt, heißt ** Q-Funktion **.

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)

Parameterschätzung durch EM-Algorithmus

Der ** EM-Algorithmus ** wird als Methode verwendet, um die wahrscheinlichste Schätzung der Parameter eines Wahrscheinlichkeitsmodells mit versteckten Variablen zu finden. Der EM-Algorithmus besteht aus zwei Schritten.

Erwartungsschritt: Eine Operation, die einen erwarteten Wert für eine versteckte Variable annimmt. Die Modellparameter werden festgelegt und die hintere Wahrscheinlichkeit versteckter Variablen berechnet.
Maximierungsschritt: Finden Sie den Modellparameter, der die Q-Funktion maximiert, unter Verwendung der posterioren Wahrscheinlichkeit der im E-Schritt erhaltenen versteckten Variablen.

Der EM-Algorithmus wiederholt die Schritte E und M, bis sich die logarithmische Wahrscheinlichkeit der vollständigen Daten nicht mehr ändert. Da der Konvergenzwert vom Anfangswert abhängt, kann nur die lokale optimale Lösung erhalten werden. Es ist notwendig, den Anfangswert zu ändern und den EM-Algorithmus mehrmals zu verwenden, um den wahrscheinlichsten geschätzten Wert zu erhalten und den besten aus ihnen auszuwählen. Der EM-Algorithmus kann wie folgt zusammengefasst werden.

Initialisierung von $ \ pi_k, \ \ boldsymbol \ mu_k, \ \ boldsymbol \ Sigma_k $
Schritt E: Schätzen Sie die hintere Wahrscheinlichkeit $ \ gamma (z_ {ik}) $ unter Verwendung der aktuellen Parameter $ \ 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)}$
M-Schritt: Modellparameter mit geschätzten $ \ gamma (z_ {ik}) $ neu schätzen $ \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} $
Beurteilung beenden: Wenn die logarithmische Wahrscheinlichkeit vollständiger Daten nicht konvergiert hat, berechnen Sie neu aus 2. Wenn es konvergiert hat, endet es.

Implementierung in Python

Ich habe es wie folgt implementiert.

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

Ergebnis

Die folgenden Ergebnisse wurden erhalten. Es scheint, dass eine vernünftige Schätzung vorgenommen wurde.

abschließend

Ich konnte eine gemischte Normalverteilung in Python implementieren. Da es einige Teile gibt, die ich nicht verstehe, wie die Ableitungsmethode jedes Parameters im M-Schritt und die Art des EM-Algorithmus. Ich werde es in Zukunft studieren und hinzufügen.

Implementierung einer gemischten Normalverteilung in Python