Hi! I am using a Hierachical Mixture model (trying the linear local Dirichlet model described here, see page 10). I have xdata in the form of an array with entries (N, Z), and about 2000 datapoints. The model is defined as follows:
mixture_model = pm.Model()
with mixture_model:
sigma = pm.HalfNormal("sigma", sigma=10, shape=n_models) # Tried gamma, gamma slow. HalfNormal with 10 -- representative of approx. rms
beta = pm.Normal("beta", mu=0, sigma=1, shape=(n_models, 2))
x = pm.Data("x", xdata)
logalpha = pm.math.matmul(beta, x.T)
weights = pm.Dirichlet("w", pm.math.exp(logalpha.T)) #, a=np.ones(2))#
y_obs = pm.NormalMixture("y", w=weights, mu=mu, sigma=sigma, observed=BEA, shape=1)
The problem I am having is an extremely slow sampling time, without very good results either. Running ~1000 draws and 1000 tuning for four chains takes about 1-2 hours, and running more, say 25 000 draws and 25 000 tuning, takes about 24 hours, but returns and ESS of ~4, and rhat >~ 3. I should also point out that I have managed to remove any divergences by adjusting the prior distributions of sigma to the reasonable sigma=10, but nonetheless the sampling is very slow.
I realize that the 4000 weight parameters here are a lot, but since the paper I am following must have more, as they are using multiple models, I am wondering whether I might have made some mistake.
I am thankful for any help that you can give. Please let me know if any additional information is needed.