Jak wyodrębnić nienadzorowane klastry z procesu Dirichlet w PyMC3?

Question

Właśnie ukończyłem książkę Bayesian Analysis in Python autorstwa Osvaldo Martin (świetna książka do zrozumienia pojęć bayesowskich i fantazyjnych indeksów numpy).

Naprawdę chcę rozszerzyć moje zrozumienie na bayesowskie modele miksowania dla nienadzorowanego grupowania próbek. Wszystkie moje wyszukiwania google doprowadziły mnie do Austin Rochford's tutorial, która jest naprawdę pouczająca. Rozumiem, co się dzieje, ale Nie jestem jasne, w jaki sposób można go dostosować do klastrowania (zwłaszcza przy użyciu wielu atrybutów dla przypisania klastra, ale to jest inny temat).

Rozumiem, jak przypisywać priors dla Dirichlet distribution, ale nie mogę dowiedzieć się, jak uzyskać klastry w PyMC3. Wygląda na to, że większość z mus zbiegają się do centroidów (tj. Środków z dystrybucji, z których próbkowałem), ale nadal są one oddzielne components. Pomyślałem o zrobieniu cutoff dla weights (w w modelu), ale nie wydaje się działać tak, jak sobie wyobrażałem, ponieważ wiele components mają nieco różniące się średnie parametry mus, które są zbieżne.

Jak mogę wyodrębnić klastry (centroidy) z tego modelu PyMC3? Dałem mu maksymalnie 15 komponentów, które chcę zjednoczyć do 3. mus wydaje się być we właściwym miejscu, ale wagi są pomieszane b/c są one dystrybuowane między innymi klastrami, więc nie mogę użyć progu wagi (chyba że je scalę, ale nie sądzę, że tak to jest zwykle jest zrobione).

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

Podobne pytania poniżej

https://stats.stackexchange.com/questions/120209/pymc3-dirichlet-distribution dla regresji i nie grupowanie

https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process teoria na temat procesu DP

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

Dirichlet process in PyMC 3 kieruje mnie do samouczka Austin Rochford w powyższym

Answer 1

Korzystanie kilka nowych-owski dodatków do pymc3 pomogą uczynić to jasne. Myślę, że zaktualizowałem przykład procesu Dirichlet po ich dodaniu, ale wydaje się, że powrócił do starej wersji podczas czyszczenia dokumentacji; Naprawię to wkrótce.

Jedną z trudności jest to, że dane, które wygenerowałeś, są znacznie bardziej rozproszone niż te, które są dostępne na elementach składowych; jeśli ustandaryzujesz swoje dane, próbki powinny mieszać się znacznie szybciej.

Po drugie, pymc3 obsługuje teraz dystrybucje mieszanek, w których zmienna wskaźnika component została zmarginalizowana. Te marginalne rozkłady mieszanin pomogą przyspieszyć mieszanie i pozwolą na użycie NUTS (zainicjowanego przy pomocy ADVI).

Wreszcie, dzięki tym okrojonym wersjom nieskończonych modeli, przy napotykaniu problemów obliczeniowych często przydatne jest zwiększenie liczby potencjalnych komponentów. Zauważyłem, że K = 30 działa lepiej dla tego modelu niż K = 15.

Poniższy kod implementuje te zmiany i pokazuje, jak można wyodrębnić "aktywny" składnik.

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

Standardized data

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

Widzimy, że trzy elementy wydają się być używane, a ich waga są dość zbliżone do prawdziwych wartości.

Mixture weights

Wreszcie widzimy, że tylne oczekiwać pomocy tych trzech elementów pokrywa się z prawdą (znormalizowany) oznacza dość dobrze.

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

tablica ([- 0,73763891, -0,17284594, 2,10423978])

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

tablica ([- 0,73017789, -0,16765707, 2,0824262])