如何从PyMC3中的Dirichlet过程中提取无监督的聚类？

我刚刚完成了Osvaldo Martin的Python书中的贝叶斯分析(理解贝叶斯概念和一些花哨的numpy索引的好书).

我真的想将我的理解扩展到贝叶斯混合模型,用于无监督的样本聚类.我所有的谷歌搜索都让我看到了Austin Rochford的教程,这本教程非常有用.我理解发生了什么,但我不清楚它如何适应群集(特别是使用群集分配的多个属性,但这是一个不同的主题).

我知道如何分配先验,Dirichlet distribution但我无法弄清楚如何获得集群PyMC3.看起来大多数mus会聚到质心(即我从中采样的分布方式),但它们仍然是分开的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 for regression and not clustering

关于DP流程的https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process理论

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

PyMC 3中的Dirichlet过程指导我上面的Austin Rochford教程

使用一些新增功能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);

标准化数据

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

我们看到似乎使用了三个组件,并且它们的权重合理地接近真实值.

混合物重量

最后,我们看到这三个组成部分的后验预期方法与真实(标准化)方法相当匹配.

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

数组([ - 0.73763891,-0.17284594,2.10423978])

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

数组([ - 0.73017789,-0.16765707,2.0824262])