Problems with Hierarchical Normal Mixture Model

[iavicenna]

What is the mu that goes into NormalMixture? It does not seem to be defined in the model. Also how many clusters are we talking about? Is it really 4000 or did you broadcast your weights into some shape compatible with the data and that is why it is so big? Cant tell what is going on in your dirichlet definition cause I dont know the shape of logalpha but for instance if you have 2000 datapoints and you want to define two clusters per point with possibly weights depending on data points that should not be a major cause of slow down. That would mean your weights will have shape 2000,2 (I think) which is what a weight of shape (2,) would have been broadcasted to too. Perhaps you can write the shape of your observed, mu and weights for us to get a better idea (or better yet post a complete code).

In my experience for instance when running about 200 data points, 9 clusters on a six dimensional problem with ~3000 draws and tune can take like up to an hour on a pretty decent computer. Is your problem 1D? Never worked with 1D problems much but maybe it is more tame there.

Also normal mixtures are very prone to mixing so you need to use tricks like ordering (and for higher dimensional problems, on a possibly choice of coordinates more suited for your problem), see for instance:

Why does MAP vastly outperform sample in bayesian clustering?

Typically, if you are getting better results with a “simpler” method (MAP, metropolis, etc), your model has some kind of mis-specification. One of the best and most frustrating things about NUTS is that it fails very loudly. I’m a beginner when it comes to mixture models, but I spent some time playing with your model. As written, you end up with really nonsense results in the non-sorted dimension. After some consultation with @ricardoV94 , I think it’s really important that you use a mixture of…

Also to possibly understand whether or not it is a problem with your data not being compatible with the model or not, you can maybe generate data using sample_prior_predictive and use the data generated to train the model and see if it runs faster. See:

Bayesian Non-negative Matrix Factorization using a Gamma-Poisson model

One general point on workflow: I strongly suggest you ditch generating test data with numpy, and instead use your PyMC model directly to do it. It’s cleaner and helps avoid errors that can crop up whenever one has to do the same thing twice (to say nothing of subtle differences in parameterization/broadcasting between pymc and numpy). Here’s how things could look in your case: import pymc as pm import arviz as az import numpy as np seed = sum(map(ord, 'Gamma Poisson Matrix')) rng = np.random.…