python - ¿Cómo extraer los clústeres no supervisados ​​de un proceso de Dirichlet en PyMC3?

machine-learning bayesian

Question

Acabo de terminar el análisis bayesiano en el libro de Python de Osvaldo Martin (gran libro para comprender los conceptos bayesianos y algunos indicios de indiferencia).

Realmente quiero ampliar mi comprensión a los modelos de mezclas bayesianas para la agrupación no supervisada de muestras. Todas mis búsquedas en google me llevaron al tutorial de Austin Rochford, que es realmente informativo. Entiendo lo que está sucediendo, pero no tengo claro cómo se puede adaptar a la agrupación (especialmente al utilizar atributos múltiples para las asignaciones de clúster, pero ese es un tema diferente).

Entiendo cómo asignar los priores para la Dirichlet distribution pero no puedo entender cómo obtener los clústeres en PyMC3 . Parece que la mayoría de los mus convergen a los centroides (es decir, los medios de las distribuciones de las que tomé muestras) pero siguen siendo components separados. Pensé en hacer un punto de corte para las weights ( w en el modelo), pero parece que no funciona de la manera que imaginaba, ya que los components múltiples tienen parámetros medios ligeramente diferentes que convergen.

¿Cómo puedo extraer los clusters (centroides) de este modelo PyMC3 ? Le di un máximo de 15 componentes con los que quiero convergir 3 . Parece que el mus se encuentra en el lugar correcto, pero los pesos están mal porque se están distribuyendo entre los otros grupos, así que no puedo usar un umbral de peso (a menos que los combine, pero no creo que así sea normalmente se hace).

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import multiprocessing
import seaborn as sns
import pandas as pd
import theano.tensor as tt
%matplotlib inline
# Clip at 15 components
K = 15
# Create mixture population
centroids = [0, 10, 50]
weights = [(2/5),(2/5),(1/5)]
mix_3 = np.concatenate([np.random.normal(loc=centroids[0], size=int(150*weights[0])), # 60 samples
np.random.normal(loc=centroids[1], size=int(150*weights[1])), # 60 samples
np.random.normal(loc=centroids[2], size=int(150*weights[2]))])# 30 samples
n = mix_3.size

# Create and fit model
with pm.Model() as Mod_dir:
alpha = pm.Gamma(''alpha'', 1., 1.)
beta = pm.Beta(''beta'', 1., alpha, shape=K)
w = pm.Deterministic(''w'', beta * tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]]))
component = pm.Categorical(''component'', w, shape=n)
tau = pm.Gamma("tau", 1.0, 1.0, shape=K)
mu = pm.Normal(''mu'', 0, tau=tau, shape=K)
obs = pm.Normal(''obs'',
mu[component],
tau=tau[component],
observed=mix_3)
step1 = pm.Metropolis(vars=[alpha, beta, w, tau, mu, obs])
# step2 = pm.CategoricalGibbsMetropolis(vars=[component])
step2 = pm.ElemwiseCategorical([component], np.arange(K)) # Much, much faster than the above
tr = pm.sample(1e4, [step1, step2], njobs=multiprocessing.cpu_count())
#burn-in = 1000, thin by grabbing every 5th idx
pm.traceplot(tr[1e3::5])

Preguntas similares a continuación

https://stats.stackexchange.com/questions/120209/pymc3-dirichlet-distribution para regresión y no agrupamiento

Teoría de https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process en el proceso de DP

https://stats.stackexchange.com/questions/116311/draw-a-multinomial-distribution-from-a-dirichlet-distribution explica DP

El proceso Dirichlet en PyMC 3 me dirige al tutorial de Austin Rochford arriba

Answer 1

Usar un par de adiciones nuevas a pymc3 ayudará a aclarar esto. Creo que actualicé el ejemplo del Proceso Dirichlet después de que se agregaron, pero parece haber sido revertido a la versión anterior durante la limpieza de la documentación; Lo arreglaré pronto.

Una de las dificultades es que los datos que ha generado están mucho más dispersos que los anteriores sobre los medios de los componentes; si estandariza sus datos, las muestras se mezclarán mucho más rápido.

El segundo es que pymc3 ahora admite distribuciones de mezcla donde el component variable del indicador ha sido marginado. Estas distribuciones de mezclas marginales ayudarán a acelerar la mezcla y le permitirán usar NUTS (inicializado con ADVI).

Finalmente, con estas versiones truncadas de modelos infinitos, cuando se encuentran con problemas de computación, a menudo es útil aumentar el número de componentes potenciales. He encontrado que K = 30 funciona mejor para este modelo que K = 15 .

El siguiente código implementa estos cambios y muestra cómo se puede extraer el componente "activo".

from matplotlib import pyplot as plt
import numpy as np
import pymc3 as pm
import seaborn as sns
from theano import tensor as T
blue = sns.color_palette()[0]
np.random.seed(462233) # from random.org
N = 150
CENTROIDS = np.array([0, 10, 50])
WEIGHTS = np.array([0.4, 0.4, 0.2])
x = np.random.normal(CENTROIDS[np.random.choice(3, size=N, p=WEIGHTS)], size=N)
x_std = (x - x.mean()) / x.std()
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist(x_std, bins=30);

Datos estandarizados

K = 30
with pm.Model() as model:
alpha = pm.Gamma(''alpha'', 1., 1.)
beta = pm.Beta(''beta'', 1., alpha, shape=K)
w = pm.Deterministic(''w'', beta * T.concatenate([[1], T.extra_ops.cumprod(1 - beta)[:-1]]))
tau = pm.Gamma(''tau'', 1., 1., shape=K)
lambda_ = pm.Uniform(''lambda'', 0, 5, shape=K)
mu = pm.Normal(''mu'', 0, tau=lambda_ * tau, shape=K)
obs = pm.NormalMixture(''obs'', w, mu, tau=lambda_ * tau,
observed=x_std)
with model:
trace = pm.sample(2000, n_init=100000)
fig, ax = plt.subplots(figsize=(8, 6))
ax.bar(np.arange(K) - 0.4, trace[''w''].mean(axis=0));

Vemos que se utilizan tres componentes y que sus ponderaciones están razonablemente cerca de los valores verdaderos.

Pesos de mezcla

Finalmente, vemos que los medios esperados posteriores de estos tres componentes coinciden bastante bien con los verdaderos (estandarizados).

trace[''mu''].mean(axis=0)[:3]

array ([- 0.73763891, -0.17284594, 2.10423978])

(CENTROIDS - x.mean()) / x.std()

array ([- 0.73017789, -0.16765707, 2.0824262])