Problem with fitting Multivariate Mixture of Gaussians

Questions
2

Initially, I thought the problem is the function multivariatenormal
return pm.MvNormal.dist(mu, cov) if dist else pm.MvNormal('mvn' + suffix, mu, cov) where a shape should be explicitly specified. But running the code it turns out it is not the case. The problem is that the Mixture is trying to compute a mode by evaluating the logp, but did not catch the wrong shape in the mode. For now, wrap the mixture logp as a DensityDist or a Potential should allow you to continue with the sampling etc.

from pymc3.distributions.dist_math import bound
from pymc3.math import logsumexp

K = 3
with pm.Model() as model:
    #alpha = pm.Gamma('alpha', 1., 1.)
    #beta = pm.Beta('beta', 1, alpha, shape=K)
    #w = pm.Deterministic('w', stick_breaking(beta))
    w = pm.Dirichlet('w', np.ones(K))
    comp_dists = [multivariatenormal(init_mean, init_sigma, init_corr, k, dist=True) for k in range(K)]
    logpcomp = tt.stack([comp_dist.logp(theano.shared(data)) for comp_dist in comp_dists], axis=1)
    pm.Potential('logp', logsumexp(tt.log(w) + logpcomp, axis=-1).sum())

I will try to send a bugfix in the meantime.