Thanks @junpenglao!
Most of the times the problem is bad mixing (it has to be the tricky case, obviously). I’ve tried using DFM’s helper function and tune the off-diagonal entires, but I’m actually getting worse results.
My current understanding is that there are “multiple possible solutions”, and the bad mixing is an indication of the sampler finding different solutions from different starting points - is this a plausible explanation? Some of these solutions don’t make much sense; to constrain the search space some more (in addition to “tightening” the priors) I’m trying to add order constraints using pm.Potential('a',X-Y). Other ideas on how to improve mixing would be much appreciated. 