如何从 PyMC3 中的狄利克雷过程中提取无监督集群？

Question

我刚刚读完 Bayesian Analysis in Python book by Osvaldo Martin（了解贝叶斯概念和一些奇特的 numpy 索引的好书）。

我真的很想将我的理解扩展到用于无监督样本聚类的贝叶斯混合模型。我所有的 google 搜索都让我找到了 Austin Rochford's tutorial，这非常有用。我明白发生了什么，但 我不清楚如何将其适应聚类 （尤其是对聚类分配使用多个属性，但这是一个不同的主题）。

我了解如何为 Dirichlet distribution 分配先验，但我不知道如何在 PyMC3 中获取聚类。看起来 mus 的大部分都收敛到质心（即我从中采样的分布的均值），但它们仍然是分开的 components。我考虑过为 weights（模型中的 w）设置一个截止点，但这似乎并不像我想象的那样有效，因为多个 components 的平均参数 [=14] 略有不同=] 正在收敛。

如何从这个 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-processDP过程理论

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

Dirichlet process in PyMC 3 将我定向到上面 Austin Rochford 的教程

Answer 1

在 pymc3 的基础上添加几个 new-ish 将有助于阐明这一点。我想我在添加 Dirichlet Process 示例后更新了它们，但在文档清理期间它似乎已恢复为旧版本；我会尽快解决的。

困难之一是您生成的数据比组件均值上的先验可以容纳的数据分散得多；如果您对数据进行标准化，样本混合的速度应该会快得多。

第二个是 pymc3 现在支持指标变量 component 被边缘化的混合分布。这些边际混合分布将有助于加速混合并允许您使用 NUTS（使用 ADVI 初始化）。

最后，使用这些无限模型的截断版本，在遇到计算问题时，增加潜在组件的数量通常很有用。我发现 K = 30 比 K = 15.

更适合这个模型

以下代码实现了这些更改并展示了如何提取 "active" 分量均值。

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]

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

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

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

如何从 PyMC3 中的狄利克雷过程中提取无监督集群？

How to extract unsupervised clusters from a Dirichlet Process in PyMC3?

python

machine-learning

bayesian

unsupervised-learning

pymc3