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.
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
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)
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