Модель многомерной смеси PyMC3: ограничение компонентов не пустыми

Question

Я внедряю персонализированную смесь многомерных гауссовых регрессий в pymc3 и столкнуться с проблемой с пустыми компонентами. Обратившись к связанному примеру модели смеси PyMC3, я попытался реализовать модель с использованием одномерных нормалей, но у меня также были некоторые проблемы.

Я испробовал несколько стратегий, чтобы каждый компонент не был пустым, но каждый потерпел неудачу. Это показано в коде ниже. Мой конкретный вопрос: каков наилучший способ заставить все компоненты быть непустыми в смеси многомерных гауссианов, используя pymc3 ?

Обратите внимание, что попытка №1 в приведенном ниже коде взята из Mixture Model в Примере PyMC3 и здесь не работает.

Вы можете скопировать синтетические данные, которые я использую, с помощью функции в этом разделе.

import pymc3 as pm
import numpy as np
import theano
import theano.tensor as T
from scipy import stats

# Extract problem dimensions.
N = X.shape[0]  # number of samples
F = X.shape[1]  # number of features
pids = I[:, 0].astype(np.int)  # primary entity ids
uniq_pids = np.unique(pids)  # array of unique primary entity ids
n_pe = len(uniq_pids)  # number of primary entities

with pm.Model() as gmreg:
    # Init hyperparameters.
    a0 = 1
    b0 = 1
    mu0 = pm.constant(np.zeros(F))
    alpha = pm.constant(np.ones(K))
    coeff_precisions = pm.constant(1 / X.var(0))

    # Init parameters.
    # Dirichlet shape parameter, prior on indicators.
    pi = pm.Dirichlet(
        'pi', a=alpha, shape=K)

    # ATTEMPT 1: Make probability of membership for each cluter >= 0.1
    # ================================================================
    pi_min_potential = pm.Potential(
        'pi_min_potential', T.switch(T.min(pi) < .1, -np.inf, 0))
    # ================================================================

    # The multinomial (and by extension, the Categorical), is a symmetric
    # distribution. Using this as a prior for the indicator variables Z
    # makes the likelihood invariant under the many possible permutations of
    # the indices. This invariance is inherited in posterior inference.
    # This invariance model implies unidentifiability and induces label
    # switching during inference.

    # Resolve by ordering the components to have increasing weights.
    # This does not deal with the parameter identifiability issue.
    order_pi_potential = pm.Potential(
        'order_pi_potential',
        T.sum([T.switch(pi[k] - pi[k-1] < 0, -np.inf, 0)
               for k in range(1, K)]))

    # Indicators, specifying which cluster each primary entity belongs to.
    # These are draws from Multinomial with 1 trial.
    init_pi = stats.dirichlet.rvs(alpha.eval())[0]
    test_Z = np.random.multinomial(n=1, pvals=init_pi, size=n_pe)
    as_cat = np.nonzero(test_Z)[1]
    Z = pm.Categorical(
        'Z', p=pi, shape=n_pe, testval=as_cat)

    # ATTEMPT 2: Give infinite negative likelihood to the case
    # where any of the clusters have no users assigned.
    # ================================================================
    # sizes = [T.eq(Z, k).nonzero()[0].shape[0] for k in range(K)]
    # nonempty_potential = pm.Potential(
    #     'comp_nonempty_potential',
    #     np.sum([T.switch(sizes[k] < 1, -np.inf, 0) for k in range(K)]))
    # ================================================================

    # ATTEMPT 3: Add same sample to each cluster, each has at least 1.
    # ================================================================
    # shared_X = X.mean(0)[None, :]
    # shared_y = y.mean().reshape(1)
    # X = T.concatenate((shared_X.repeat(K).reshape(K, F), X))
    # y = T.concatenate((shared_y.repeat(K), y))

    # Add range(K) on to the beginning to include shared instance.
    # Z_expanded = Z[pids]
    # Z_with_shared = T.concatenate((range(K), Z_expanded))
    # pid_idx = pm.Deterministic('pid_idx', Z_with_shared)
    # ================================================================

    # Expand user cluster indicators to each observation for each user.
    pid_idx = pm.Deterministic('pid_idx', Z[pids])

    # Construct masks for each component.
    masks = [T.eq(pid_idx, k).nonzero() for k in range(K)]
    comp_sizes = [masks[k][0].shape[0] for k in range(K)]

    # Component regression precision parameters.
    beta = pm.Gamma(
        'beta', alpha=a0, beta=b0, shape=(K,),
        testval=np.random.gamma(a0, b0, size=K))
    # Regression coefficient matrix, with coeffs for each component.
    W = pm.MvNormal(
        'W', mu=mu0, tau=T.diag(coeff_precisions), shape=(K, F),
        testval=np.random.randn(K, F) * std)

    # The mean of the observations is the result of a regression, with
    # coefficients determined by the cluster the sample belongs to.

    # Now we have K different multivariate normal distributions.
    X = T.cast(X, 'float64')
    y = T.cast(y, 'float64')
    comps = []
    for k in range(K):
        mask_k = masks[k]

        X_k = X[mask_k]
        y_k = y[mask_k]

        n_k = comp_sizes[k]
        precision_matrix = beta[k] * T.eye(n_k)

        comp_k = pm.MvNormal(
            'comp_%d' % k,
            mu=T.dot(X_k, W[k]), tau=precision_matrix,
            observed=y_k)
        comps.append(comp_k)

Первые два подхода не обеспечивают непустые кластеры; пытаясь сэмплировать результаты в LinAlgError:

with gmreg:
step1 = pm.Metropolis(vars=[pi, beta, W])
step2 = pm.ElemwiseCategoricalStep(vars=[Z], values=np.arange(K))
tr = pm.sample(100, step=[step1, step2])
...:     

Failed to compute determinant []
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-2-c7df53f4c6a5> in <module>()
      2     step1 = pm.Metropolis(vars=[pi, beta, W])
      3     step2 = pm.ElemwiseCategoricalStep(vars=[Z], values=np.arange(K))
----> 4     tr = pm.sample(100, step=[step1, step2])
      5 

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/sampling.pyc in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    155         sample_args = [draws, step, start, trace, chain,
    156                        tune, progressbar, model, random_seed]
--> 157     return sample_func(*sample_args)
    158 
    159 

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/sampling.pyc in _sample(draws, step, start, trace, chain, tune, progressbar, model, random_seed)
    164     progress = progress_bar(draws)
    165     try:
--> 166         for i, strace in enumerate(sampling):
    167             if progressbar:
    168                 progress.update(i)

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/sampling.pyc in _iter_sample(draws, step, start, trace, chain, tune, model, random_seed)
    246         if i == tune:
    247             step = stop_tuning(step)
--> 248         point = step.step(point)
    249         strace.record(point)
    250         yield strace

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/step_methods/compound.pyc in step(self, point)
     12     def step(self, point):
     13         for method in self.methods:
---> 14             point = method.step(point)
     15         return point

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/step_methods/arraystep.pyc in step(self, point)
     87             inputs += [point]
     88 
---> 89         apoint = self.astep(bij.map(point), *inputs)
     90         return bij.rmap(apoint)
     91 

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/step_methods/gibbs.pyc in astep(self, q, logp)
     38 
     39     def astep(self, q, logp):
---> 40         p = array([logp(v * self.sh) for v in self.values])
     41         return categorical(p, self.var.dshape)
     42 

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/blocking.pyc in __call__(self, x)
    117 
    118     def __call__(self, x):
--> 119         return self.fa(self.fb(x))

/home/mack/anaconda/lib/python2.7/site-packages/pymc3/model.pyc in __call__(self, *args, **kwargs)
    423     def __call__(self, *args, **kwargs):
    424         point = Point(model=self.model, *args, **kwargs)
--> 425         return self.f(**point)
    426 
    427 compilef = fastfn

/home/mack/anaconda/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    604                         self.fn.nodes[self.fn.position_of_error],
    605                         self.fn.thunks[self.fn.position_of_error],
--> 606                         storage_map=self.fn.storage_map)
    607                 else:
    608                     # For the c linker We don't have access from

/home/mack/anaconda/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    593         t0_fn = time.time()
    594         try:
--> 595             outputs = self.fn()
    596         except Exception:
    597             if hasattr(self.fn, 'position_of_error'):

/home/mack/anaconda/lib/python2.7/site-packages/theano/gof/op.pyc in rval(p, i, o, n)
    766             # default arguments are stored in the closure of `rval`
    767             def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 768                 r = p(n, [x[0] for x in i], o)
    769                 for o in node.outputs:
    770                     compute_map[o][0] = True

/home/mack/anaconda/lib/python2.7/site-packages/theano/tensor/nlinalg.pyc in perform(self, node, (x,), (z,))
    267     def perform(self, node, (x,), (z, )):
    268         try:
--> 269             z[0] = numpy.asarray(numpy.linalg.det(x), dtype=x.dtype)
    270         except Exception:
    271             print 'Failed to compute determinant', x

/home/mack/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in det(a)
   1769     """
   1770     a = asarray(a)
-> 1771     _assertNoEmpty2d(a)
   1772     _assertRankAtLeast2(a)
   1773     _assertNdSquareness(a)

/home/mack/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in _assertNoEmpty2d(*arrays)
    220     for a in arrays:
    221         if a.size == 0 and product(a.shape[-2:]) == 0:
--> 222             raise LinAlgError("Arrays cannot be empty")
    223 
    224 

LinAlgError: Arrays cannot be empty
Apply node that caused the error: Det(Elemwise{Mul}[(0, 1)].0)
Inputs types: [TensorType(float64, matrix)]
Inputs shapes: [(0, 0)]
Inputs strides: [(8, 8)]
Inputs values: [array([], shape=(0, 0), dtype=float64)]

Backtrace when the node is created:
  File "/home/mack/anaconda/lib/python2.7/site-packages/pymc3/distributions/multivariate.py", line 66, in logp
    result = k * T.log(2 * np.pi) + T.log(1./det(tau))

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

... который указывает, что компонент пуст, так как матрица точности имеет форму (0, 0),

Третий метод фактически решает проблему пустых компонентов, но дает очень странное поведение вывода. Я выбрал выжигание на основе графиков и прореживал до каждого десятого образца. Образцы все еще сильно автокоррелированы, но гораздо лучше, чем без разбавления. В этот момент я суммировал значения Z по выборкам, и вот что я получаю:

In [3]: with gmreg:
    step1 = pm.Metropolis(vars=[pi, beta, W])
    step2 = pm.ElemwiseCategoricalStep(vars=[Z], values=np.arange(K))
    tr = pm.sample(1000, step=[step1, step2])
   ...:     
 [-----------------100%-----------------] 1000 of 1000 complete in 258.8 sec

...

In [24]: zvals = tr[300::10]['Z']

In [25]: np.array([np.bincount(zvals[:, n]) for n in range(nusers)])
Out[25]: 
array([[ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70],
       [ 0,  0, 70]])

Поэтому по какой-то причине все пользователи назначаются в последний кластер для каждого образца.

Answer 1

Я столкнулся с подобной проблемой. Нечто подобное сработало для смеси многомерных гауссовских моделей. Что касается того, является ли это лучшим, это, безусловно, лучшее решение, которое я нашел.

pm.Potential('pi_min_potential', T.switch(
                                           T.all(
                                                [pi[i, 0] < 0.1 for i in range(K)]), -np.inf, 0))

Ключевым моментом здесь является то, что вам необходимо учитывать каждый потенциал, который ниже вашего предела. Далее вам следует отрегулировать shape вашей pi Распределение, как указано в комментариях. Это повлияет на вашу индексацию в T.switch вызов (на pi[i,0]).