After examining my trace plots, I can clearly see that the posterior distribution of a certain parameter is bimodal. However, R-hat statistics are 1. How is it even possible if the chains are obviously not mixing?
You may have indeed hit the stationary distribution and you’re sampling from it, but you’re just not sampling it efficiently. It looks like your chains are able to move amongst the modes but only after several iterations. I’m guessing your ESS is probably low compared to the number of iterations. I bet if you ran for only 100 samples the Rhats would be way worse.
ESS is indeed lower than for other parameters, but according to bin/diagnose of cmdstan, all measures are satisfactory.
Rhat looks at the magnitude of variablility between the chains as compared to the magnitude of variability within the chains. So when all chains explore both modes, as is the case here, Rhat of 1 is expected.
This makes sense. But then, taking the average of all four chains seems fatally wrong.