Die gemischte Gaußsche Verteilung ist ein Modell der gewichteten Addition mehrerer Gaußscher Verteilungen und kann eine multimodale Verteilung darstellen. Sie können auch den EM-Algorithmus verwenden, um die Parameter dieses Modells zu optimieren. Der erwartete Wert der latenten Variablen wird durch den E-Schritt dieses EM-Algorithmus berechnet. Wenn die Berechnung jedoch zu diesem Zeitpunkt naiv durchgeführt wird, können Probleme mit Überlauf und Unterlauf auftreten. Zu diesem Zeitpunkt kann dieses Problem durch Verwendung einer bekannten numerischen Berechnungsmethode namens logsumexp vermieden werden.
In Bezug auf die Parameteroptimierung der gemischten Gaußschen Verteilung durch den EM-Algorithmus wird nur der Teil angezeigt, der sich auf logsumexp bezieht. Weitere Informationen finden Sie unter "Mustererkennung und maschinelles Lernen" 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? Dh = UTF8 & refRID = 1GS4SDYG8VEYR40YSADB). (Vorerst werde ich Ihnen die Materialien geben, als ich SlideShare PRML 9.0-9.2 einen Vortrag hielt.)
\begin{aligned}
p(x|\theta) &= \sum_k \pi_k N(x_n|\mu_k,\Sigma_k)\\
\pi_k&:Mischverhältnis
\end{aligned}
Die Verteilung der beobachteten Variablen $ x $ der gemischten Gaußschen Verteilung wird durch die obige Gleichung ausgedrückt. Um den EM-Algorithmus anzuwenden, führen wir die latente Variable $ z $ ein und drücken die gleichzeitige Verteilung von $ x und z $ wie folgt aus.
\begin{aligned}
p(x,z) &= \prod_k \pi_k^{z_k} N(x_n|\mu_k,\Sigma_k)^{z_k}\\
z&: 1ofK
\end{aligned}
Mit diesen kann der erwartete Wert der latenten Variablen z aus dem Bayes'schen Theorem wie folgt berechnet werden.
\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}
Im E-Schritt der Optimierung durch den gemischten Gaußschen EM-Algorithmus wird der erwartete Wert dieser latenten Variablen berechnet. Wenn wir uns den Nenner ansehen, können wir sehen, dass bei der Exponentialberechnung der Gaußschen Verteilung ein Unterlauf auftreten kann, selbst wenn die Operation auf der logarithmischen Skala ausgeführt wird, da es eine Summenoperation der Gaußschen Verteilung gibt. Daher ist es notwendig, etwas zu entwickeln, und eine Methode ist logsumexp.
logsumexp
\log(\sum^N_{i=1} \exp(x_i))
Selbst wenn das Ergebnis dieser Berechnung selbst nicht über- oder unterläuft, kann ein einzelnes $ \ exp (x_i) $ über- oder unterlaufen. Daher wird die Formel wie folgt transformiert.
\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}
Wenn Sie die Berechnung auf diese Weise durchführen, fällt das Argument der Exponentialoperation in den Bereich von $ [x_ {min} - x_ {max}, 0] $, sodass die Möglichkeit eines Überlaufs und Unterlaufs stark verringert wird.
Dieser logsumexp ist auch in Pythons Bibliothek für maschinelles Lernen scikitlearn implementiert und einfach zu verwenden (Scikit Learn Utilities for Developers. )). Die Implementierung von logsumexp in scicit learn lautet außerdem sklearn / utils / extmath.py. Es ist im Folgenden. Sie können sehen, dass es genau der oben gezeigten Formel entspricht, mit Ausnahme des Teils, der der Achse entspricht.
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
Basierend auf dem oben Gesagten habe ich das gemischte Gaußsche Modell implementiert, also werde ich es verlassen (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)
In logsumexp werden viele Operationen von log und exp angezeigt, sodass die Berechnungsgeschwindigkeit langsam zu sein scheint (logsumexp ist die schwarze Geschichte der Menschheit). Wenn Geschwindigkeit zu einem Engpass wird, ist es möglicherweise besser, eine andere Methode zu verwenden.
Recommended Posts