I read this paper presenting a Bayesian strategy for inferring the form of the Milky Way’s distribution function given a mix of complete and incomplete velocity data for its satellites by Eadie et al. on astro ph this week. The methodology presented follows closely Gelman et al.’s Bayesian Data Analysis, including a rigorous checking of the Gibbs sampling MCMC convergence and an examination of the Bayesian coverage under the assumed model and a model mis-specification. My one ‘complaint’ would be that the presentation seems overly pedagogical given that Bayesian methods are already quite thoroughly adopted within this branch of astronomy (e.g. Magorrian’s Dirichlet Process mixture model version).

Some food for thought is included in the ‘future work’ section of the paper which ponders solutions to the problem that some popular theoretical models for the distribution function are unnormalizable (i.e., they effectively postulate an infinite halo with an infinite amount of mass in its tails). The natural solution would be to simply truncate these profiles at a sensible radius: maybe even allowing this truncation as a nuisance parameter of the model. But I wonder whether there’s any way to work directly with these unnormalised sampling distributions using something like the exchange algorithm. Probably not, but I’d be interested to hear suggestions.

Having read about a problem of inferring the distribution function of a physical system evolving under Hamiltonian dynamics gave a natural segue to the CSML reading group lecture this week by Michael Betancourt on Adiabatic Monte Carlo: a version of Hamiltonian Monte Carlo for exploring (e.g.) the thermodynamic path from prior to posterior. While the technique proposed seems to offer a clever & sophisticated solution to some of the difficulties of ordinary thermodynamic integration (e.g. the choice of temperature sequence and difficulties navigating a discrete set of transitions without losing equilibrium) I wonder how challenging it will prove to incorporate this within a general modeling tool like STAN, given the need for diagnostics and work around on the two types of metastability encountered.

The range of models to which this technique is applicable, while large, is obviously limited in some ways not clearly elaborated in the paper. Certainly the sample space must be a topological space, but then more specifically it must have locally Euclidean geometry almost everywhere, so we’re looking at R^N, and then I suspect the likelihood must also be continuous with connected support. Interestingly the last condition is precisely not true for general distribution models of gravitationally stable stellar systems considered for multi-component potentials in the field of the Eadie et al. paper—although the parameter space can be placed in R^N it contains islands of allowed distribution functions separated by tracts of strictly disallowed distribution functions.

One further thought I had was whether the adiabatic Monte Carlo dynamics can be run over the space of partial data posteriors for iid data problem a la Chopin 2001. Presumably so if the target family is set to be something along the lines of for ?