我刚刚完成 Bayesian Analysis in Python通过 Osvaldo Martin 预订(了解贝叶斯概念和一些奇特的 numpy 索引的好书)。

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

我了解如何为 Dirichlet 分布 分配先验，但我不知道如何在 PyMC3 中获取簇。看起来大部分 mus 都收敛到质心(即我从中采样的分布的均值)，但它们仍然是独立的 组件。我考虑过为 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 上面的教程

最佳答案

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

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

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

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

关于python - 如何从 PyMC3 中的狄利克雷过程中提取无监督集群？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/41553988/