[Here is a reply to my comments on THAMES sent by the first author of the paper, Martin Metodiev. The above replica of the cover of Rivers of London is obviously unrelated with the reply or the original blog, beyond presenting a fantasy map of the Thames!]
Thank you for your review of our article! Adapting your previous work in this field has been a pleasure. Before I respond to your comments, I would like to emphasize that the simplicity of our estimator lies in its simple analytic expression (a truncated harmonic mean of reciprocal unnormalized posterior density values). Indeed, our package “thamesmix” (recently submitted to CRAN!) has a function to compute the marginal likelihood of any mixture model. This function requires only two parameters: the unnormalized log-posterior function (the logarithm of the prior plus the log-likelihood) and the MCMC simulations from the posterior.
Regarding your main comments:
1. “the evacuation of earlier methods as not simple or not universal enough is rather disingenuous. For instance, software that do not return (latent) allocation vectors can easily be post-processed.”
I could not find an example of post-process simulations on top of MCMC outputs applied to compute these methods. It sounds really interesting, and I would be happy to cite it. Is there a reference that you can recommend?
In any case, the point still stands. Most estimators which we cite with regards to this point do not just need allocation samplers, but also the analytic expressions of the distribution of the allocation vectors or the distribution of the data conditional on these allocation vectors that come with them. I do not think that a closed form of this distribution is available in general.
2.“the handling of the label switching issue—the reason why Larry Wasserman saw mixtures at the same magnitude of evil as tequila!—is problematic for several reasons.”
The fact that our estimator is invariant to label-switching is indeed the core of our method. The simple Gibbs sampler gets stuck in one mode, and this is why the classical version of bridge sampling is biased by a factor of G! in the simulation setting. As you point out, this is successfully resolved when using fully symmetric bridge sampling in the experiment section. However, the computation cost of this fully symmetric estimator rises super-exponentially with G, so I do not see how it could be evaluated for G=15, where the number of symmetric modes is equal to 15! (over one trillion). One of the main points of our article is that the symmetric THAMES can be evaluated in a feasible amount of time, even in such a high-dimensional multivariate setting.
3. “the (legitimate) purpose of using marginal likelihoods for selecting the number G of components is weakened by the intrusion of alternate proposals to assess G from the data”
I would like to point out that these alternate proposals do not in any way impact the definition of the THAMES. It is the simple definition given in Equation (5). They are only used to speed up the computation.
4. “several mentions are made of the other estimators being biased, which is indeed the case for bridge sampling (if not necessarily for importance sampling), but not necessarily a central issue”
The problem that we see with the classical, non-symmetric bridge sampling method in the setting of mixture models is not simply that it is biased. The problem is that the bias is persistent and often roughly equal to the factor of G! when the MCMC sampler failed to switch between modes. We have not had this experience with the THAMES: it converged even when the MCMC was stuck.
Martin Metodiev and his coauthor(es)s have produced 
