So, it’s probably no secret to my close collaborators that this paper by Martin Weinberg (who is, indisputably, a highly talented astro-statistician & cosmologist) really Wangs my Chung. In my opinion, anyone searching for a generally-reliable marginal likelihood estimator using only samples from the MCMC posterior is seriously tilting at windmills. We know the harmonic mean estimator is a disaster and we know why very well (cf. Wolpert & Schmidler 2011). We also know of a few improvements (e.g. the Gelfand-Dey or one of the two Newton & Raftery stabilisations), but we do not have (and I expect never will have) any solution to its real convergence rate issues.
In effect, all the Weinberg paper does then is to propose yet another stablisation (based on volume-tesselation and truncation), which he freely admits, “does not, and indeed cannot, circumvent the limitations of the HMA“. Yet we are then later asked to accept the sweeping proposition that: “Together, the NLA and VTA provide robust alternatives to computing the marginal likelihood with errors of delta log Z < 0.5 for models with dimensionality k < 20.” For the target audience of astronomers and physicists skimming over the technical details of marginal likelihood computation and Lebesgue integration one suspects it will be the latter, not the former, ending up as the take-home message.
Then there is the fact that Weinberg’s method is embedded within his Bayesian Inference Engine software, which explores the posterior via simulated tempering between the prior and posterior. In any sensible Universe it is this complete collection of tempered draws that would be forming the basis of BIE’s marginal likelihood estimation routine (via e.g. reverse logistic regression or importance tempering). Instead, in this Bizarro Universe these valuable draws are simply discarded in favour of a HME proxy using only the posterior samples. Of course, the power of these alternative approaches is never explored in the BA paper; the presented comparisons being against an order one estimator, the Laplace approximation.
But what really grinds my gears is the circularity of the initial discussion. We’re given a an MCMC sample from the posterior and told (mistakenly attributed to Newton & Raftery 1994?) that such a sample is distributed as Z × π(θ|y) = π(y|θ) × π(θ), after which Weinberg asks under what additional conditions does the HME integral for Z exist (i.e., remain finite). But, in order for ergodicity such that the MCMC chain has a unique limiting distribution, its target, f(θ) [proportional to π(y|θ) × π(θ)], must already be proper in the first place (Tierney 1994), meaning that Z is very much finite with no additional conditions. Pointless, right?
UPDATE: After some discussion on ASAIP it seems that Martin & I are converging towards some points of agreement …