Как извлечь неконтролируемые кластеры из процесса Дирихле в PyMC3?

Question

Я только что закончил книгу Bayesian Analysis in Python от Osvaldo Martin (отличная книга, чтобы понять байесовские концепции и некоторые фантазии индексации).

Я действительно хочу расширить свое понимание до байесовских моделей смеси для неконтролируемой кластеризации образцов. Все мои поисковые запросы привели меня к Austin Rochford's tutorial, который действительно информативен. Я понимаю, что происходит, но Непонятно, как это можно адаптировать для кластеризации (особенно с использованием нескольких атрибутов для кластерных назначений, но это другая тема).

Я понимаю, как назначить приоритеты для Dirichlet distribution, но я не могу понять, как получить кластеры в PyMC3. Похоже, что большинство из mus сходятся к центроидам (т. Е. Средствам выборок I, отбираемых из), но они все еще являются отдельными components. Я думал о том, чтобы сделать отсечение для weights (w в модели), но это, похоже, не работает так, как я себе представлял, поскольку несколько components имеют несколько разные средние параметры mus, которые сходятся.

Как извлечь кластеры (центроиды) из этой модели PyMC3? Я дал ему максимум 15 компонентов, которые я хочу свести на 3. mus, кажется, находятся в правильном месте, но весы перепутаны b/c, они распределяются между другими кластерами, поэтому я не могу использовать порог веса (если я не объединю их, но я не думаю, что так оно и есть обычно делается).

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

Подобные вопросы ниже

https://stats.stackexchange.com/questions/120209/pymc3-dirichlet-distribution для регрессии и не кластеризацию

https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process теория о процессе DP

https://stats.stackexchange.com/questions/116311/draw-a-multinomial-distribution-from-a-dirichlet-distribution объясняет DP

Dirichlet process in PyMC 3 направляет меня учебник Остин Рокфорд в выше

Answer 1

Используя пару новых иш дополнений pymc3 поможет сделать это ясно. Я думаю, что я обновил пример Dirichlet Process после их добавления, но, похоже, он был возвращен к старой версии во время очистки документации; Я исправлю это в ближайшее время.

Одна из трудностей заключается в том, что данные, которые вы создали, намного более разбросаны, чем приспосабливаются к средним компонентам; если вы стандартизируете свои данные, образцы должны смешиваться намного быстрее.

Вторым является то, что pymc3 теперь поддерживает распределение смеси, где индикаторная переменная component была изолирована. Эти распределения предельных смесей помогут ускорить смешивание и позволяют использовать NUTS (инициализируется с помощью ADVI).

И наконец, с этими усеченными версиями бесконечных моделей при столкновении с вычислительными проблемами часто бывает полезно увеличить количество потенциальных компонентов. Я обнаружил, что K = 30 работает лучше для этой модели, чем K = 15.

Следующий код реализует эти изменения и показывает, как можно извлечь «активное» компонентное средство.

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

Мы видим, что три компонента по всей видимости, будет использоваться, и что их веса достаточно близки к истинным значениям.

Mixture weights

Наконец, мы видим, что задние ожидали средства этих трех компонентов соответствуют истинным (нормированный) означает достаточно хорошо.

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

массив ([- 0,73763891, -0,17284594, 2,10423978])

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

массив ([- 0,73017789, -0,16765707, 2,0824262])