[Python] Modèle gaussien mixte avec Pyro

J'ai essayé d'estimer un modèle gaussien mixte avec Pyro. Basé sur Exemple officiel, il est exécuté avec un contenu supplémentaire, le cas échéant.

• Le code source de cet article s'exécute sur Jupyter Notebook. <détails>
  Environnement </ summary> Windows10 Python: 3.7.7 Jupyter Notebook: 1.0.0 PyTorch: 1.5.1 Pyro: 1.4.0 scipy: 1.5.2 numpy: 1.19.1 matplotlib: 3.3.0 seaborn: 0.10.1

• J'ai écrit d'autres articles sur Pyro, alors jetez un œil si vous le souhaitez. [Python] Estimation bayésienne avec Pyro

import os

import numpy as np
from scipy import stats

import torch
import pyro
import pyro.distributions as dist
from pyro import poutine
from pyro.infer.autoguide import AutoDelta
from pyro.optim import Adam
from pyro.infer import SVI, TraceEnum_ELBO, config_enumerate, infer_discrete

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

pyro.set_rng_seed(0)
pyro.enable_validation(True)

Préparation des données

Appelez le jeu de données iris à partir de seaborn et définissez la valeur de petal_length comme données cible.

df = sns.load_dataset('iris')
data = torch.tensor(df['petal_length'], dtype=torch.float32)
sns.swarmplot(data=df, x='petal_length')

En regardant l'intrigue, il semble bon de diviser le cluster en deux [^ 1]

paramètres du modèle

Dans Pyro, le modèle de distribution est décrit dans la méthode du modèle. Appliquez le modèle gaussien mixte avec chaque cluster de données $ x_1, \ cdots, x_n \ in \ mathbb {R} $ comme $ z_1, \ cdots, z_n \ in \ {1, \ cdots, K \} $ Je vais.

\begin{align}
p &\sim Dir(\tau_0/K, \cdots, \tau_0/K) \\
z_i &\sim Cat(p) \\

\mu_k &\sim N(\mu_0, \sigma_0^2) \\
\sigma_k &\sim InvGamma(\alpha_0, \beta_0) \\
x_i &\sim N(\mu_{z_i}, \sigma_{z_i}^2)
\end{align}

$ K $ est le nombre de clusters, $ \ tau_0, \ mu_0, \ sigma_0, \ alpha_0, \ beta_0 $ sont des paramètres de distribution antérieurs. [^ 2]
Estimations bayésiennes de $ \ mu_1, \ cdots, \ mu_K $ et $ \ sigma_1, \ cdots, \ sigma_K $, et crée un modèle qui calcule de manière probabiliste les clusters $ z_1, \ cdots, z_n $.

K = 2  # Fixed number of clusters
TAU_0 = 1.0
MU_0 = 0.0
SIGMA_0_SQUARE = 10.0
ALPHA_0 = 1.0
BETA_0 = 1.0

@config_enumerate
def model(data):
    alpha = torch.full((K,), fill_value=TAU_0)
    p = pyro.sample('p', dist.Dirichlet(alpha))
    with pyro.plate('cluster_param_plate', K):
        mu = pyro.sample('mu', dist.Normal(MU_0, SIGMA_0_SQUARE))
        sigma = pyro.sample('sigma', dist.InverseGamma(ALPHA_0, BETA_0))

    with pyro.plate('data_plate', len(data)):
        z = pyro.sample('z', dist.Categorical(p))
        pyro.sample('x', dist.Normal(locs[z], scales[z]), obs=data)

@ config_enumerate est un décorateur pour l'échantillonnage de variables discrètes` pyro.sample ('z', dist.Categorical (p)) ʻen parallèle.

Vérifiez la valeur échantillonnée

En utilisant poutine.trace, vous pouvez vérifier la valeur d'échantillonnage lorsque les données sont données à model.

trace_model = poutine.trace(model).get_trace(data)
tuple(trace_model.nodes.keys())

> ('_INPUT',
   'p',
   'cluster_param_plate',
   'mu',
   'sigma',
   'data_plate',
   'z',
   'x',
   '_RETURN')

Le type de trace_model.nodes est ʻOrderedDict, qui contient la clé ci-dessus. _INPUT fait référence aux données fournies à model, _RETURN fait référence à la valeur de retour de model(Aucune dans ce cas), et les autres se réfèrent aux paramètres définis dansmodel`.

En guise de test, vérifions la valeur de p. C'est un paramètre échantillonné à partir de $ Dir (\ tau_0 / K, ⋯, \ tau_0 / K) $. trace_model.nodes ['p'] est aussi dict, et vous pouvez voir la valeur avec valeur.

trace_model.nodes['p']['value']

> tensor([0.8638, 0.1362])

Ensuite, vérifions la valeur de z, qui représente le cluster de chaque donnée.

trace_model.nodes['z']['value']

> tensor([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1,
          0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
          0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
          0, 0, 0, 0, 1, 0])

On peut dire que 0 est facilement échantillonné à partir de la valeur de «p», mais le résultat est exactement cela. Notez qu'il s'agit d'un échantillonnage d'une distribution antérieure, nous ne pouvons donc toujours pas faire une estimation correcte.

paramètres du guide

Dans Pyro, définissez la distribution postérieure dans le guide. pyro.infer.autoguide.AutoDelta est une classe d'estimation MAP.

guide = AutoDelta(poutine.block(model, expose=['p', 'mu', 'sigma']))

poutine.block est une méthode qui sélectionne les paramètres à estimer. ʻAutoDelta ne semble pas être capable de gérer le paramètre discret z`, donc il n'est pas spécifié par expose. L'estimation de «z» est effectuée après l'ajustement de la distribution.

Ce guide dicte les paramètres de la valeur estimée des données.

guide(data)

> {'p': tensor([0.5000, 0.5000], grad_fn=<ExpandBackward>),
   'mu': tensor([4.0607, 2.8959], grad_fn=<ExpandBackward>),
   'sigma': tensor([1.3613, 1.6182], grad_fn=<ExpandBackward>)}

Pour le moment, il ne renvoie que la valeur initiale, mais à partir de maintenant, nous retournerons la valeur estimée MAP en ajustant avec SVI.

Raccord de distribution

Dans le guide, j'ai construit un modèle qui estime d'autres paramètres sans estimer «z». En d'autres termes, «z» doit être marginalisé et calculé. Pour ce faire, définissez la perte de l'estimation de la variation probabiliste sur TraceEnum_ELBO ().

optim = pyro.optim.Adam({'lr': 1e-3})
svi = SVI(model, guide, optim, loss=TraceEnum_ELBO())

Faites un essayage.

NUM_STEPS = 3000
pyro.clear_param_store()

history = []
for step in range(1, NUM_STEPS + 1):
    loss = svi.step(data)
    history.append(loss)
    if step % 100 == 0:
        print(f'STEP: {step} LOSS: {loss}')

Le graphique de la perte à chaque étape ressemble à ceci:

plt.figure()
plt.plot(history)
plt.title('Loss')
plt.grid()
plt.xlim(0, 3000)
plt.show()

On peut juger que la valeur de la perte a convergé et l'estimation est terminée.

Confirmation de la distribution estimée

Obtenez des estimations de $ p, \ mu, \ sigma $ à partir de guide.

map_params = guide(data)
p = map_params['p']
mu = map_params['mu']
sigma = map_params['sigma']
print(p)
print(mu)
print(sigma)

> tensor([0.6668, 0.3332], grad_fn=<ExpandBackward>)
  tensor([4.9049, 1.4618], grad_fn=<ExpandBackward>)
  tensor([0.8197, 0.1783], grad_fn=<ExpandBackward>)

Tracez la distribution. Dans la figure ci-dessous, le graphique avec la marque x signifie la valeur des données.

x = np.arange(0, 10, 0.01)
y1 = p[0].item() * stats.norm.pdf((x - mu[0].item()) / sigma[0].item())
y2 = p[1].item() * stats.norm.pdf((x - mu[1].item()) / sigma[1].item())

plt.figure()
plt.plot(x, y1, color='red', label='z=0')
plt.plot(x, y2, color='blue', label='z=1')
plt.scatter(data.numpy(), np.zeros(len(data)), color='black', alpha=0.3, marker='x')
plt.legend()
plt.show()

La distribution peut être bien estimée.

Estimation de cluster

Tout d'abord, définissez les paramètres estimés par «guide» sur «modèle». Dans Pyro, les paramètres sont définis via trace.

trace_guide_map = poutine.trace(guide).get_trace(data)
model_map = poutine.replay(model, trace=trace_guide_map)

Vérifiez les paramètres définis dans model. Ici, seul $ \ mu $ est confirmé.

trace_model_map = poutine.trace(model_map).get_trace(data)
trace_guide_map.nodes['mu']['value']

>> tensor([4.9048, 1.4618], grad_fn=<ExpandBackward>)

Il correspond à la valeur de $ \ mu $ dans guide. Estimez ensuite la valeur de $ z $ pour chaque donnée. À ce stade, utilisez pyro.infer.infer_discrete.

model_map = infer_discrete(model_map, first_available_dim=-2)

first_available_dim = -2 est un paramètre pour éviter les conflits avec la dimension de data_plate. Cela définit l'estimation $ z $ sur model, qui peut être obtenue à partir de trace.

trace_model_map = poutine.trace(model_map).get_trace(data)
z_inferred = trace_model_map.nodes['z']['value']
z_inferred

> tensor([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
          1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0])

Tracons les données pour chaque valeur de $ z $.

df['z'] = trace_model_map.nodes['z']['value']
df['z'] = df['z'].apply(lambda z: f'z={z}')
sns.swarmplot(data=df, x='petal_length', y='z')

Vous pouvez voir que cela peut être bien estimé.

en conclusion

J'ai fait un gauss mixte avec Pyro et j'ai essayé de l'adapter. Je suis habitué à la pensée orientée objet, donc j'ai trouvé un peu fastidieux d'utiliser poutine.trace pour récupérer la valeur estimée. Lors de son utilisation, il semble préférable de créer une classe comme GaussianMixtureModel et de décrire le processus d'extraction de la valeur en interne. Je continuerai à toucher Pyro pour approfondir ma compréhension.

[^ 1]: Les données d'iris sont à l'origine un jeu de données de classification à 3 classes, mais nous ne considérerons pas la classe d'origine ici. [^ 2]: Dans l'exemple de Pyro, LogNormal est appliqué en tant que distribution de $ \ sigma_k $, mais cette fois, InverseGamma, qui est une distribution a priori conjuguée pour une échelle de distribution gaussienne, est appliqué.