To add, I now have checked, scipy.stats.mode is making its calculations based on only 30 bins binning, which is super law value for a range of values my posterior is spanning over. Still it would be interesting to know how mode is calculated inside the PyMC3, how reliable it is and how to calculate it with the tools of PyMC3, rather than third-party libs, which results in such discrepancies.