In order to estimate how well one estimated the mean of a distribution one could look at the Monte Carlo Standard Error (see Stan reference manual Section 16.4.3).

However in case of an asymmetric (or multimodal) distribution this seems to be unreasonable.

Is there an alternative way to estimate how well one has estimated the most probable value of the resulting parameter distribution.

Part of the issue is that the mean is not the most probable value. Indeed for continuous parameter spaces there isnâ€™t a well-defined notion of â€śmost probableâ€ť because all individual points have the same probability: zero!

characterizes the â€ścentralityâ€ť of the posterior distribution. The median provides a similar quantification. For a multimodal distribution that central value will likely be somewhere in between the modes, which is technically correct but also not a great summary of any individual modes.

To isolate a particular mode we need expectations values that vanish outside of the range of that mode. For example letâ€™s say that we have a one-dimensional parameter space \theta \in \mathbb{R} and we know that one model is centered at -5 and one mode is centered at +7 and there is negligible posterior probability in between. To isolate the negative mode we might consider a truncated mean,

\mathbb{E}_{\pi}[ \theta (1 - H (\theta)) ]
\approx
\frac{1}{N} \sum_{n = 1}^{N} \theta_{n} (1 - H (\theta_{n}))

where the step function H (\theta) is zero if \theta < 0 and one if \theta \ge 0. Likewise we can summarize the positive mode by

\mathbb{E}_{\pi}[ \theta H (\theta) ]
\approx
\frac{1}{N} \sum_{n = 1}^{N} \theta_{n} \cdot H (\theta_{n}).

To implement this we have to know where the modes are, but this will be true for any summary that can target individual modes.