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

Posted 2023-02-19

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

【英文标题】：How to extract unsupervised clusters from a Dirichlet Process in PyMC3?

【发布时间】：2017-05-24 01:35:47

【问题描述】：

我刚刚完成了 Osvaldo Martin 的 Bayesian Analysis in Python 书（理解贝叶斯概念和一些花哨的 numpy 索引的好书）。

我真的很想将我的理解扩展到用于无监督样本聚类的贝叶斯混合模型。我所有的谷歌搜索都把我带到了Austin Rochford's tutorial，这真的很丰富。我了解正在发生的事情，但我不清楚这如何适应集群（尤其是使用多个属性进行集群分配，但这是一个不同的主题）。

我了解如何为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 用于回归而非聚类

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 的教程

【问题讨论】：

Edward 可能有一个使用变分推理来处理狄利克雷过程混合的示例。 edwardlib.org

我会检查一下，看看我是否能弄清楚如何移植它！谢谢。我从未听说过爱德华，但到目前为止看起来很酷。

这是您要找的吗？ pymc-devs.github.io/pymc3/notebooks/dp_mix.html

@rafaelvalle 我将上面的链接作为 Austin Rochford 的教程。它解释了如何使用狄利克雷过程，但没有解释如何使用它进行聚类。我尝试逐步完成教程并在最后一步对其进行调整以获得集群的数量，但我无法让它工作。

【参考方案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 ])

【讨论】：

哇，这太不可思议了。我还没有看到pm.NormalMixture，但我喜欢它！有趣的是，tau*lambda_ 的性能比 tau 好得多。我需要复习一下我的统计数据。最后一个问题，如果你不知道有 3 个集群，你会为权重设置一个截止值吗（例如，超过 1e-3 的任何东西都是一个集群）？如果是这样，您是否推荐一个好的经验法则来确定截止值？再次感谢您，这非常有用。

这可能是我会做的，不幸的是我没有一个好的经验法则。

另外，pymc3 documentation 已根据这些更改进行了更新。

是的，这是唯一真正的区别。

我不完全确定您要做什么，但是使用 trace['w'] 应该可以找到具有最大 weight 的组件。