如何在PyMC3中从狄利克雷过程中提取无监督聚类？

Question

我刚刚完成了Python中的贝叶斯分析一书，作者是Osvaldo Martin（这是一本非常好的书，可以帮助您理解贝叶斯概念和一些花式numpy索引）。

我真的想扩展我的理解，以便于使用贝叶斯混合模型对样本进行无监督聚类。我所有的谷歌搜索都导向了Austin Rochford的教程，它真的很有启发性。我理解正在发生的事情，但我不清楚如何将其适应于聚类（特别是在使用多个属性进行聚类分配时，但这是一个不同的话题）。

我知道如何为Dirichlet分布分配先验概率，但我不知道如何在PyMC3中获取聚类。看起来大多数mus会收敛到质心（即我从中抽样的分布的均值），但它们仍然是独立的组件。我考虑过为权重（模型中的w）制定截止值，但这似乎并不像我想象的那样有效，因为多个具有略微不同平均参数mu的组件正在收敛。

我该如何从这个PyMC3模型中提取聚类（质心）？我将其设置为最多有15个组件，并希望收敛到3个。mus似乎位于正确位置，但权重混乱了，因为它们分布在其他聚类之间，所以我不能使用权重阈值（除非我合并它们，但我认为这不是正常做法）。

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 将我指向了Austin Rochford的上面的教程

- O.rka

---

Edward可能有一些使用变分推断进行狄利克雷过程混合的示例。http://edwardlib.org/ - rafaelvalle

我会去看一下，看看能否弄清楚如何移植它！谢谢。我从未听说过Edward，但到目前为止似乎很酷。 - O.rka

这是你要找的吗？ https://pymc-devs.github.io/pymc3/notebooks/dp_mix.html - rafaelvalle

@rafaelvalle 我在上面提供了Austin Rochford的教程链接。它解释了如何使用狄利克雷过程，但是并没有解释如何用它进行聚类。我尝试逐步按照教程操作，并在最后一步进行调整以获取聚类数量，但我无法使其正常运行。 - O.rka

Answer 1

使用pymc3的一些新近添加功能可以帮助澄清问题。我认为在这些功能被添加后，我更新了狄利克雷过程的示例，但在文档清理期间似乎已经恢复到旧版本了。我会尽快修复。

其中一个困难是，您生成的数据比组件均值的先验分布更分散。如果标准化数据，则样本应该更快地混合。

第二个困难是，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);

{{链接1：标准化数据}}

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

- Austin Rochford

---

哇，这太不可思议了。我还没有见过pm.NormalMixture，但我喜欢它！有趣的是，使用tau*lambda_比仅使用tau表现要好得多。我需要再恶补一下统计学知识。最后一个问题，如果你不知道有3个聚类，你会设置权重的截止值吗（例如，任何大于1e-3的都是一个聚类）？如果是这样，你推荐一个好的经验法则来确定截止值吗？再次感谢，这非常有用。 - O.rka

那可能是我会做的事情，不幸的是我没有一个好的经验法则。 - Austin Rochford

此外，pymc3文档已经根据这些更改进行了更新。 - Austin Rochford

是的，那是唯一真正的区别。 - Austin Rochford

我不完全确定你想做什么，但使用trace['w']应该可以找到权重最大的组件。 - Austin Rochford