I am trying to define a model that looks like
nvars, ncats = 3, 5
with pm.Model():
α = pm.Exponential('α', lam=λ, shape=(nvars, ncats))
b = pm.Dirichlet('b', a=α)
trace = pm.sample()
This fails, because the Dirichlet distribution expects an array, not a matrix, in PyMC3 (and in numpy, and in scipy).
Has anyone had any success in modelling something like this? I guess I can use loops…
