I am trying to sample from a posterior distribution that, I believe, is multimodal. When I sample 5 randomly initialized chains, they converge into two different values. I think this happens because, in each sampled chain, the NUTS sampler gets stuck in a sharp density region of the posterior (one of the modes) and fails to explore the rest of the density areas. If this is the problem, an additional drawback can occur: during the warm-up phase, the mass matrix and the leapfrog step-size are only locally optimized and, consequently, might not be suitable for the other density curvatures.
Is there any efficient sampler, implemented in PyMC, that is suitable for sampling from a multimodal posterior distribution? Or is my best option to reparameterize the model?
Thank you in advance!