Я пытаюсь использовать процессы Дирихле для идентификации кластеров в моих двоичных данных.Я использую учебник в качестве отправной точки, но учебник создан там, где результат представляет собой смесь одномерных нормальных или распределенных по Пуассону переменных.
У меня есть несколько бинарных переменных на наблюдение, 5 в приведенном ниже примере кода, и я не могу понять, как создать окончательный шаг смешивания.Из математического описания в этого отчета я вижу, что общая вероятность - это просто произведение вероятностей во всех назначенных кластерах.
Я не сформировал явно метки кластера (используя Categorical(w)), поскольку распределение pm.Mixture обрабатывает это, но не могу понять, как сформулировать вероятность как модель вероятности, которую понимает pymc3.
import numpy as np
import pandas as pd
import pymc3 as pm
from matplotlib import pyplot as plt
import seaborn as sns
from theano import tensor as tt
N = 100
P = 5
K_ACT = 3
# Simulate 5 variables with 100 observations of each that fit into 3 groups
mu_actual = np.array([[0.7, 0.8, 0.2, 0.1, 0.5],
[0.3, 0.4, 0.9, 0.8, 0.6],
[0.1, 0.2, 0.3, 0.2, 0.3]])
cluster_ratios = [0.6, 0.2, 0.2]
df = np.concatenate([np.random.binomial(1, mu_actual[0, :], size=(int(N*cluster_ratios[0]), P)),
np.random.binomial(1, mu_actual[1, :], size=(int(N*cluster_ratios[1]), P)),
np.random.binomial(1, mu_actual[2, :], size=(int(N*cluster_ratios[2]), P))])
# Deterministic function for stick breaking
def stick_breaking(beta):
portion_remaining = tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]])
return beta * portion_remaining
K_THRESH = 20
with pm.Model() as model:
# The DP priors to obtain w, the cluster weights
alpha = pm.Gamma('alpha', 1., 1.)
beta = pm.Beta('beta', 1, alpha, shape=K_THRESH)
w = pm.Deterministic('w', stick_breaking(beta))
# Each variable should have a probability parameter for each cluster
mu = pm.Beta('mu', 1, 1, shape=(K_THRESH, P))
obs = pm.Mixture('obs', w, pm.Bernoulli.dist(mu), observed=df)
with model:
step = pm.Metropolis()
trace = pm.sample(100, step=step, random_seed=17)
pm.traceplot(trace, varnames=['alpha', 'w'])
РЕДАКТИРОВАТЬ 28/01/2019
Я предоставил пользовательскую функцию правдоподобия, которая просто вычисляет вероятность смеси Бернулли после рисования меток компонентов из категориального распределения.Однако, хотя модель сейчас что-то делает, она не может идентифицировать 3 группы, а только находит 2. Я не могу сказать, нужна ли ей просто выборка / более эффективная параметризация, или определение модели некорректно.
import numpy as np
import pandas as pd
import pymc3 as pm
from matplotlib import pyplot as plt
import seaborn as sns
from theano import tensor as tt
N = 1000
P = 5
# Simulate 5 variables with 1000 observations of each that fit into 3 groups
mu_actual = np.array([[0.7, 0.8, 0.2, 0.1, 0.5],
[0.3, 0.4, 0.9, 0.8, 0.6],
[0.1, 0.2, 0.3, 0.2, 0.3]])
cluster_ratios = [0.4, 0.3, 0.3]
df = np.concatenate([np.random.binomial(1, mu_actual[0, :], size=(int(N*cluster_ratios[0]), P)),
np.random.binomial(1, mu_actual[1, :], size=(int(N*cluster_ratios[1]), P)),
np.random.binomial(1, mu_actual[2, :], size=(int(N*cluster_ratios[2]), P))])
# Deterministic function for stick breaking
def stick_breaking(beta):
portion_remaining = tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]])
return beta * portion_remaining
K_THRESH = 20
def bernoulli_mixture_loglh(comp, mus):
# K = maximum number clusters
# N = number observations (1000 here)
# P = number predictors (5 here)
# Shape of tensors:
# comp: K
# mus: K, P
# value (data): N, P
def loglh_(value):
mus_comp = mus[comp, :]
# These are (NxK) matrices giving likelihood contributions
# from each observation according to each component's probability
# parameter (mu)
pos = value * tt.log(mus_comp)
neg = (1-value) * tt.log((1-mus_comp))
comb = pos + neg
overall_sum = tt.sum(comb)
return overall_sum
return loglh_
with pm.Model() as model:
# The DP priors to obtain w, the cluster weights
alpha = pm.Gamma('alpha', 1., 1.)
beta = pm.Beta('beta', 1, alpha, shape=K_THRESH)
w = pm.Deterministic('w', stick_breaking(beta))
component = pm.Categorical('component', w, shape=N)
# Each variable should have a probability parameter for each cluster
mu = pm.Beta('mu', 1, 1, shape=(K_THRESH, P))
obs = pm.DensityDist('obs', bernoulli_mixture_loglh(component, mu), observed=df)
n_samples = 5000
burn = 500
thin = 10
with model:
step1 = pm.Metropolis(vars=[alpha, beta, w, mu])
step2 = pm.ElemwiseCategorical([component], np.arange(K_THRESH))
trace_ = pm.sample(n_samples, [step1, step2], sample=17)
trace = trace_[burn::thin]
pm.traceplot(trace)
plt.show()
Должен ли след быть стационарным?
# Plot weights per component
fig, ax = plt.subplots(figsize=(8, 6))
plot_w = np.arange(K_THRESH) + 1
ax.bar(plot_w - 0.5, trace['w'].mean(axis=0), width=1., lw=0);
ax.set_xlim(0.5, K_THRESH);
ax.set_xlabel('Component');
ax.set_ylabel('Posterior expected mixture weight');
Почему нет?т компоненты заказаны?В других примерах, которые я видел, они обычно располагаются в порядке убывания размера
Приведенный ниже код отображает параметр mu Bernoulliзначения для 2 основных компонентов, но они далеки от фактических значений
# Posterior values of mu for the non-zero components
mean_w = np.mean(trace['w'], axis=0)
nonzero_component = np.where(mean_w > 0.3)[0]
mean_mu = np.mean(trace['mu'], axis=0)
print(mean_mu[nonzero_component, :])
[[0.47587256 0.50065195 0.51081395 0.57693177 0.40762681]
[0.42596485 0.69626519 0.5629946 0.30185575 0.64322441]]
Фактические параметры, используемые для моделирования данных:
[[0.7, 0.8, 0.2, 0.1, 0.5],
[0.3, 0.4, 0.9, 0.8, 0.6],
[0.1, 0.2, 0.3, 0.2, 0.3]]