[Python] Gemischtes Gaußsches Modell mit Pyro

Ich habe versucht, ein gemischtes Gaußsches Modell mit Pyro zu schätzen. Basierend auf Offizielles Beispiel wird es mit entsprechendem Zusatzinhalt ausgeführt.

• Der Quellcode für diesen Artikel läuft auf Jupyter Notebook.

• Ich habe andere Artikel über Pyro geschrieben, schauen Sie also bitte nach, wenn Sie möchten. [Python] Bayesianische Schätzung mit 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)

Datenaufbereitung

Rufen Sie den Iris-Datensatz von seaborn auf und legen Sie den Wert von petal_length als Zieldaten fest.

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

Wenn man sich die Handlung ansieht, scheint es gut, den Cluster in zwei Teile zu teilen [^ 1]

Modelleinstellungen

In Pyro wird das Verteilungsmodell in der Modellmethode beschrieben. Wenden Sie das gemischte Gaußsche Modell mit jedem Datencluster $ x_1, \ cdots, x_n \ in \ mathbb {R} $ als $ z_1, \ cdots, z_n \ in \ {1, \ cdots, K \} $ an Ich werde.

\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 $ ist die Anzahl der Cluster, $ \ tau_0, \ mu_0, \ sigma_0, \ alpha_0, \ beta_0 $ sind vorherige Verteilungsparameter. [^ 2]
Bayesianische Schätzungen von $ \ mu_1, \ cdots, \ mu_K $ und $ \ sigma_1, \ cdots, \ sigma_K $ und Erstellung eines Modells, das die Cluster $ z_1, \ cdots, z_n $ wahrscheinlich berechnet.

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 ist ein Dekorator zum parallelen Abtasten diskreter Variablen pyro.sample ('z', dist.Categorical (p)).

Überprüfen Sie den Abtastwert

Mit poutine.trace können Sie den Abtastwert überprüfen, wenn Daten an model übergeben werden.

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

Der Typ von "trace_model.nodes" ist "OrderedDict", der den obigen Schlüssel enthält. "_INPUT" bezieht sich auf die Daten, die für "model" angegeben wurden, "_RETURN" bezieht sich auf den Rückgabewert von "model" (in diesem Fall keine), und die anderen beziehen sich auf die in "model" definierten Parameter.

Lassen Sie uns als Test den Wert von p überprüfen. Dies ist ein Parameter, der aus $ Dir (\ tau_0 / K, ⋯, \ tau_0 / K) $ abgetastet wurde. trace_model.nodes ['p'] ist auch dict und Sie können den Wert mit value sehen.

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

> tensor([0.8638, 0.1362])

Als nächstes überprüfen wir den Wert von "z", der den Cluster jeder Daten darstellt.

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

Man kann sagen, dass 0 leicht aus dem Wert von "p" abgetastet werden kann, aber das Ergebnis ist genau das. Beachten Sie, dass dies eine Stichprobe aus einer früheren Verteilung ist, sodass wir immer noch keine korrekte Schätzung vornehmen können.

Führungseinstellungen

Stellen Sie in Pyro die hintere Verteilung in der Führung ein. pyro.infer.autoguide.AutoDelta ist eine Klasse für die MAP-Schätzung.

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

poutine.block ist eine Methode, die die zu schätzenden Parameter auswählt. Es scheint, dass AutoDelta den diskreten Parameter z nicht verarbeiten kann, daher wird er nicht durch exponieren angegeben. Die Schätzung von "z" erfolgt nach dem Anpassen der Verteilung.

In diesem Handbuch werden die Parameter des geschätzten Werts für die Daten angegeben.

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

Im Moment wird nur der Anfangswert zurückgegeben, aber von nun an werden wir den MAP-Schätzwert durch Anpassen an SVI zurückgeben.

Verteilungsanpassung

In der Anleitung habe ich ein Modell erstellt, das andere Parameter schätzt, ohne "z" zu schätzen. Mit anderen Worten, "z" muss marginalisiert und berechnet werden. Setzen Sie dazu den Verlust der Wahrscheinlichkeitsvariationsschätzung auf "TraceEnum_ELBO ()".

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

Machen Sie eine Anpassung.

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}')

Die Darstellung des Verlusts bei jedem Schritt sieht folgendermaßen aus:

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

Es kann beurteilt werden, dass der Wert des Verlusts konvergiert hat und die Schätzung abgeschlossen ist.

Bestätigung der geschätzten Verteilung

Holen Sie sich Schätzungen von $ p, \ mu, \ sigma $ von 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>)

Zeichnen Sie die Verteilung. In der folgenden Abbildung bedeutet das Diagramm mit der x-Markierung den Wert der Daten.

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

Die Verteilung kann gut geschätzt werden.

Clusterschätzung

Stellen Sie zuerst die von guide geschätzten Parameter auf model ein. In Pyro werden die Parameter über die Ablaufverfolgung eingestellt.

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

Überprüfen Sie die in model eingestellten Parameter. Hier wird nur $ \ mu $ bestätigt.

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

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

Es entspricht dem Wert von $ \ mu $ in guide. Schätzen Sie dann den $ z $ -Wert für jedes Datenelement. Verwenden Sie zu diesem Zeitpunkt pyro.infer.infer_discrete.

model_map = infer_discrete(model_map, first_available_dim=-2)

first_available_dim = -2 ist eine Einstellung, um Konflikte mit der Dimension von data_plate zu vermeiden. Dies setzt die $ z $ -Schätzung auf "Modell", das aus der Ablaufverfolgung erhalten werden kann.

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

Zeichnen wir die Daten für jeden Wert von $ 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')

Sie können sehen, dass es gut geschätzt werden kann.

abschließend

Ich habe mit Pyro ein gemischtes Gauß gemacht und versucht, es anzupassen. Ich bin an objektorientiertes Denken gewöhnt, daher fand ich es etwas ärgerlich, "poutine.trace" zu verwenden, um den geschätzten Wert abzurufen. Wenn Sie es tatsächlich verwenden, scheint es besser, eine Klasse wie GaussianMixtureModel zu erstellen und den Prozess des internen Extrahierens des Werts zu beschreiben. Ich werde Pyro weiterhin berühren, um mein Verständnis zu vertiefen.

[^ 1]: Die Irisdaten sind ursprünglich ein 3-Klassen-Klassifizierungsdatensatz, aber wir werden die ursprüngliche Klasse hier nicht berücksichtigen. [^ 2]: In Pyros Beispiel wird LogNormal als Verteilung von $ \ sigma_k $ angewendet, aber dieses Mal wird InverseGamma angewendet, eine konjugierte vorherige Verteilung für Scala der Gaußschen Verteilung.