PRML Chapitre 10 Implémentation Python de distribution gaussienne mixte
Cette fois, j'ai implémenté la distribution gaussienne mixte variable mentionnée à la fin de l'article précédent. La dernière fois, nous avons utilisé l'algorithme EM pour estimer le plus les paramètres de la distribution gaussienne mixte, mais cette fois nous utiliserons la méthode de Bayes variante pour estimer les paramètres de la distribution gaussienne mixte. Dans l'implémentation précédente, cela définissait le nombre d'éléments mixtes dans la distribution gaussienne mixte, mais en utilisant la méthode de Bayes variante, le nombre d'éléments mixtes est automatiquement déterminé. Il y a d'autres avantages à le traiter comme un bayésien.
Variante de distribution gaussienne mixte
Dans l'estimation la plus probable d'une distribution gaussienne mixte normale, le coefficient de mélange $ \ pi_k $, la moyenne $ \ mu_k $ et la covariance $ \ Sigma_k $ (ou la précision $ \ Lambda_k $) ont été très probablement estimés, mais cette fois tous sont Estimer la distribution de probabilité comme une variable stochastique.
La figure ci-dessous (extraite de PRML Figure 10.5) est un modèle graphique du modèle utilisé cette fois.
La distribution simultanée de toutes ces variables stochastiques
p({\bf X}, {\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda}) = p({\bf X}|{\bf Z},{\bf\mu},{\bf\Lambda})p({\bf Z}|{\bf\pi})p({\bf\pi})p({\bf\mu}|{\bf\Lambda})p({\bf\Lambda})
Ce sera.
Le point de cette variante de la méthode Bayes est la distribution postérieure $ p ({\ bf Z}, {\ bf \ pi}, {\ bf \ mu}, {\ après avoir observé ** $ {\ bf X} $. \ bf \ Lambda} | {\ bf X}) $ est une approximation différentielle de la distribution de probabilité décomposée en la variable latente $ {\ bf Z} $ et les paramètres comme suit **.
\begin{align}
p({\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda}|{\bf X}) &\approx q({\bf Z}, {\bf \pi}, {\bf \mu}, {\bf \Lambda})\\
&= q({\bf Z})q({\bf\pi},{\bf\mu},{\bf\Lambda})
\end{align}
Comme il est difficile de calculer exactement la distribution postérieure d'origine, nous utilisons une approximation de variante.
Le flux général de cette variante de la méthode de Bayes est
- Initialisez les distributions de probabilité $ q ({\ bf Z}) $ et $ q ({\ bf \ pi}, {\ bf \ mu}, {\ bf \ Lambda}) $
- Mettre à jour $ q ({\ bf Z}) $ (correspond à l'étape E de l'algorithme EM)
- Mettre à jour $ q ({\ bf \ pi}, {\ bf \ mu}, {\ bf \ Lambda}) $ (correspond à l'étape M de l'algorithme EM)
- Répétez les étapes 2 et 3 jusqu'à convergence Ça ressemble à ça. Si vous notez les formules des étapes 2 et 3, ce sera un montant considérable, je vais donc l'omettre.
code
Bibliothèque
En plus des habituels matplotlib et numpy, j'utilise également scipy. Puisque nous avons besoin d'une fonction digamma et d'une fonction gamma, nous utiliserons cette fois une bibliothèque autre que Numpy pour la partie algorithme.
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import digamma, gamma
Variante mixte Gauss
class VariationalGaussianMixture(object):
def __init__(self, n_component=10, alpha0=1.):
#Nombre d'éléments mixtes dans une distribution gaussienne mixte
self.n_component = n_component
#Paramètres de la distribution antérieure du coefficient de mélange pi
self.alpha0 = alpha0
def init_params(self, X):
self.sample_size, self.ndim = X.shape
#Définir les paramètres de distribution précédents
self.alpha0 = np.ones(self.n_component) * self.alpha0
self.m0 = np.zeros(self.ndim)
self.W0 = np.eye(self.ndim)
self.nu0 = self.ndim
self.beta0 = 1.
self.component_size = self.sample_size / self.n_component + np.zeros(self.n_component)
#Initialiser les paramètres de la distribution de probabilité
self.alpha = self.alpha0 + self.component_size
self.beta = self.beta0 + self.component_size
indices = np.random.choice(self.sample_size, self.n_component, replace=False)
self.m = X[indices].T
self.W = np.tile(self.W0, (self.n_component, 1, 1)).T
self.nu = self.nu0 + self.component_size
#Renvoie les paramètres de la distribution de probabilité
def get_params(self):
return self.alpha, self.beta, self.m, self.W, self.nu
#Méthode variante Bayes
def fit(self, X, iter_max=100):
self.init_params(X)
#Mettre à jour jusqu'à ce que la distribution de probabilité converge
for i in xrange(iter_max):
params = np.hstack([array.flatten() for array in self.get_params()])
#Distribution de probabilité q(z)Mise à jour
r = self.e_like_step(X)
#Distribution de probabilité q(pi, mu, Lambda)Mise à jour
self.m_like_step(X, r)
#Si ça converge, ça finit
if np.allclose(params, np.hstack([array.ravel() for array in self.get_params()])):
break
else:
print("parameters may not have converged")
#Distribution de probabilité q(z)Mise à jour
def e_like_step(self, X):
d = X[:, :, None] - self.m
gauss = np.exp(
-0.5 * self.ndim / self.beta
- 0.5 * self.nu * np.sum(
np.einsum('ijk,njk->nik', self.W, d) * d,
axis=1)
)
pi = np.exp(digamma(self.alpha) - digamma(self.alpha.sum()))
Lambda = np.exp(digamma(self.nu - np.arange(self.ndim)[:, None]).sum(axis=0) + self.ndim * np.log(2) + np.linalg.slogdet(self.W.T)[1])
#Formule PRML(10.67)Calcul du taux de charge
r = pi * np.sqrt(Lambda) * gauss
#Formule PRML(10.49)Normaliser le taux de charge
r /= np.sum(r, axis=-1, keepdims=True)
#Au cas où 0% se produirait
r[np.isnan(r)] = 1. / self.n_component
return r
#Distribution de probabilité q(pi, mu, Lambda)Mise à jour
def m_like_step(self, X, r):
#Formule PRML(10.51)
self.component_size = r.sum(axis=0)
#Formule PRML(10.52)
Xm = X.T.dot(r) / self.component_size
d = X[:, :, None] - Xm
#Formule PRML(10.53)
S = np.einsum('nik,njk->ijk', d, r[:, None, :] * d) / self.component_size
#Formule PRML(10.58) q(pi)Mettre à jour les paramètres pour
self.alpha = self.alpha0 + self.component_size
#Formule PRML(10.60)
self.beta = self.beta0 + self.component_size
#Formule PRML(10.61)
self.m = (self.beta0 * self.m0[:, None] + self.component_size * Xm) / self.beta
d = Xm - self.m0[:, None]
#Formule PRML(10.62)
self.W = np.linalg.inv(
np.linalg.inv(self.W0)
+ (self.component_size * S).T
+ (self.beta0 * self.component_size * np.einsum('ik,jk->ijk', d, d) / (self.beta0 + self.component_size)).T).T
#Formule PRML(10.63)
self.nu = self.nu0 + self.component_size
# p(données de test|Données d'entraînement)Pi,mu,Approximatif avec l'estimation de Lambda
def predict_proba(self, X):
#Formule PRML(B.80)Valeur attendue de la matrice de précision Lambda
covs = self.nu * self.W
precisions = np.linalg.inv(covs.T).T
d = X[:, :, None] - self.m
exponents = np.sum(np.einsum('nik,ijk->njk', d, precisions) * d, axis=1)
#Formule PRML(10.38)
gausses = np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(covs.T).T * (2 * np.pi) ** self.ndim)
#Formule PRML(10.69)Valeur attendue du coefficient de mélange pi
gausses *= (self.alpha0 + self.component_size) / (self.n_component * self.alpha0 + self.sample_size)
#Périphérialisation de la variable latente z
return np.sum(gausses, axis=-1)
#Clustering
def classify(self, X):
#Distribution de probabilité q(Z)Argmax
return np.argmax(self.e_like_step(X), 1)
#Calculer la distribution t des élèves
def student_t(self, X):
nu = self.nu + 1 - self.ndim
#Formule PRML(10.82)
L = nu * self.beta * self.W / (1 + self.beta)
d = X[:, :, None] - self.m
#Formule PRML(B.72)
maha_sq = np.sum(np.einsum('nik,ijk->njk', d, L) * d, axis=1)
#Formule PRML(B.68)
return (
gamma(0.5 * (nu + self.ndim))
* np.sqrt(np.linalg.det(L.T))
* (1 + maha_sq / nu) ** (-0.5 * (nu + self.ndim))
/ (gamma(0.5 * nu) * (nu * np.pi) ** (0.5 * self.ndim)))
#Distribution prévue p(données de test|Données d'entraînement)Calculer
def predict_dist(self, X):
#Formule PRML(10.81)
return (self.alpha * self.student_t(X)).sum(axis=-1) / self.alpha.sum()
Code entier
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import digamma, gamma
class VariationalGaussianMixture(object):
def __init__(self, n_component=10, alpha0=1.):
self.n_component = n_component
self.alpha0 = alpha0
def init_params(self, X):
self.sample_size, self.ndim = X.shape
self.alpha0 = np.ones(self.n_component) * self.alpha0
self.m0 = np.zeros(self.ndim)
self.W0 = np.eye(self.ndim)
self.nu0 = self.ndim
self.beta0 = 1.
self.component_size = self.sample_size / self.n_component + np.zeros(self.n_component)
self.alpha = self.alpha0 + self.component_size
self.beta = self.beta0 + self.component_size
indices = np.random.choice(self.sample_size, self.n_component, replace=False)
self.m = X[indices].T
self.W = np.tile(self.W0, (self.n_component, 1, 1)).T
self.nu = self.nu0 + self.component_size
def get_params(self):
return self.alpha, self.beta, self.m, self.W, self.nu
def fit(self, X, iter_max=100):
self.init_params(X)
for i in xrange(iter_max):
params = np.hstack([array.flatten() for array in self.get_params()])
r = self.e_like_step(X)
self.m_like_step(X, r)
if np.allclose(params, np.hstack([array.ravel() for array in self.get_params()])):
break
else:
print("parameters may not have converged")
def e_like_step(self, X):
d = X[:, :, None] - self.m
gauss = np.exp(
-0.5 * self.ndim / self.beta
- 0.5 * self.nu * np.sum(
np.einsum('ijk,njk->nik', self.W, d) * d,
axis=1)
)
pi = np.exp(digamma(self.alpha) - digamma(self.alpha.sum()))
Lambda = np.exp(digamma(self.nu - np.arange(self.ndim)[:, None]).sum(axis=0) + self.ndim * np.log(2) + np.linalg.slogdet(self.W.T)[1])
r = pi * np.sqrt(Lambda) * gauss
r /= np.sum(r, axis=-1, keepdims=True)
r[np.isnan(r)] = 1. / self.n_component
return r
def m_like_step(self, X, r):
self.component_size = r.sum(axis=0)
Xm = X.T.dot(r) / self.component_size
d = X[:, :, None] - Xm
S = np.einsum('nik,njk->ijk', d, r[:, None, :] * d) / self.component_size
self.alpha = self.alpha0 + self.component_size
self.beta = self.beta0 + self.component_size
self.m = (self.beta0 * self.m0[:, None] + self.component_size * Xm) / self.beta
d = Xm - self.m0[:, None]
self.W = np.linalg.inv(
np.linalg.inv(self.W0)
+ (self.component_size * S).T
+ (self.beta0 * self.component_size * np.einsum('ik,jk->ijk', d, d) / (self.beta0 + self.component_size)).T).T
self.nu = self.nu0 + self.component_size
def predict_proba(self, X):
covs = self.nu * self.W
precisions = np.linalg.inv(covs.T).T
d = X[:, :, None] - self.m
exponents = np.sum(np.einsum('nik,ijk->njk', d, precisions) * d, axis=1)
gausses = np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(covs.T).T * (2 * np.pi) ** self.ndim)
gausses *= (self.alpha0 + self.component_size) / (self.n_component * self.alpha0 + self.sample_size)
return np.sum(gausses, axis=-1)
def classify(self, X):
return np.argmax(self.e_like_step(X), 1)
def student_t(self, X):
nu = self.nu + 1 - self.ndim
L = nu * self.beta * self.W / (1 + self.beta)
d = X[:, :, None] - self.m
maha_sq = np.sum(np.einsum('nik,ijk->njk', d, L) * d, axis=1)
return (
gamma(0.5 * (nu + self.ndim))
* np.sqrt(np.linalg.det(L.T))
* (1 + maha_sq / nu) ** (-0.5 * (nu + self.ndim))
/ (gamma(0.5 * nu) * (nu * np.pi) ** (0.5 * self.ndim)))
def predict_dist(self, X):
return (self.alpha * self.student_t(X)).sum(axis=-1) / self.alpha.sum()
def create_toy_data():
x1 = np.random.normal(size=(100, 2))
x1 += np.array([-5, -5])
x2 = np.random.normal(size=(100, 2))
x2 += np.array([5, -5])
x3 = np.random.normal(size=(100, 2))
x3 += np.array([0, 5])
return np.vstack((x1, x2, x3))
def main():
X = create_toy_data()
model = VariationalGaussianMixture(n_component=10, alpha0=0.01)
model.fit(X, iter_max=100)
labels = model.classify(X)
x_test, y_test = np.meshgrid(
np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
X_test = np.array([x_test, y_test]).reshape(2, -1).transpose()
probs = model.predict_dist(X_test)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap=cm.get_cmap())
plt.contour(x_test, y_test, probs.reshape(100, 100))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
if __name__ == '__main__':
main()
résultat
La vidéo ci-dessous montre la mise en cluster à l'aide de la méthode de Bayes. La couleur des points indique à quel cluster il appartient, et au fur et à mesure que l'apprentissage progresse, le nombre d'éléments mixtes valides diminue, et finalement il devient trois. (Le code ci-dessus ne montre pas une telle progression et ne montre que le résultat final)
À la fin
Le chapitre 10 de PRML présente des exemples d'utilisation de la méthode de Bayes variante pour la régression linéaire et la régression logistique ainsi que le gauss mixte, nous les mettrons donc en œuvre si nous en avons l'occasion.
Recommended Posts