How Do You Compare Models When Their Best Answers Come in Bunches?

Imagine two rival theories of a distant galaxy, each able to explain the same telescope data. To choose between them, astronomers often turn to a single number: the Bayesian evidence, also called the marginal likelihood. It measures how probable the observed data are under a whole model, after averaging over the model's parameters weighted by their prior. The ratio of two evidences, the Bayes factor, gives a quantitative score for which model the data favor.

Useful, yes. Easy, no. The evidence is a multidimensional integral over the entire parameter space. For simple models, numerical integration can work. For the complex models now common in astrophysics and cosmology, it becomes expensive or intractable. Nested sampling can compute evidence directly, but it restricts sampling to a nested strategy. Meanwhile, powerful samplers such as Hamiltonian Monte Carlo, NUTS, and newer machine-learning-assisted methods can explore posterior distributions efficiently but do not themselves provide the marginal likelihood.

A posterior distribution is the probability of a model's parameters after seeing the data. In a simple case, it is a single bump. In hard cases, it can be an archipelago: several separated islands of high probability, called modes. Such multimodal posteriors are where previous density estimators have often struggled.

The Harmonic Mean, Relearned

In earlier work, researchers introduced the learned harmonic mean estimator. It estimates the evidence using only posterior samples and their probability densities. That is important because it lets scientists use whatever sampler is best for the problem. The method builds on the classical harmonic mean estimator, which is notoriously unstable because its variance can be infinite. The learned version overcomes that by using machine learning to find an internal target density that is concentrated inside the posterior and has thinner tails.

The catch was the density estimator. Previous versions used discrete normalizing flows, neural networks that transform a simple distribution into a complex one. They can struggle with multimodal posteriors. When trained on samples from separated peaks, they tend to place probability mass in the bridges between peaks, violating the requirement that the target stays within the posterior. A previous highly multimodal gravitational-wave application required a hand-picked Gaussian mixture for the base distribution to cope.

A Flow That Learns the Whole Landscape

In a preprint posted to arXiv on 26 January 2026, Alicja Polanska and Jason D. McEwen introduce a new approach. They use flow matching-based continuous normalizing flows as the internal density estimator. A continuous normalizing flow describes a transformation between distributions using an ordinary differential equation, a path from a simple base distribution to the target. Flow matching trains the flow by teaching a neural network a time-dependent velocity field: the direction and speed that move samples along that path. Unlike earlier training schemes, flow matching is simulation-free. It does not require costly simulation of the differential equation during training, making it efficient and scalable.

The authors use a multilayer perceptron for the velocity field and train it on posterior samples. They reserve half of the posterior samples for training and half for inference. They then evaluate the learned density by solving the probability flow ODE backward. To ensure the target density is concentrated within the posterior, they scale the base distribution's variance with a temperature parameter T between 0 and 1, giving the target thinner tails. They report that this setup handles challenging multimodal posteriors without fine-tuning or heuristic modifications to the base distribution.

Two Benchmarks, and What They Do Not Yet Show

The paper presents numerical experiments, not new observations of the sky. The first is the Rastrigin function in two dimensions, a standard benchmark whose landscape has many local peaks. The authors draw posterior samples using the emcee MCMC sampler and repeat the experiment 100 times with different seeds. They find that the flow captures all represented modes, including modes with very low weighting, and that the estimated evidence is accurate with reliable error estimates. The flow was used at temperature T=0.98.

The second example is a mixture of five Gaussian components in 20 dimensions. The modes are non-overlapping, and the covariance matrices include correlations between adjacent dimensions, making the geometry deliberately difficult. Again, the authors repeat the experiment 100 times, comparing their estimates with an analytic ground truth. They report accurate evidence estimates and reasonable error estimates. The flow was used at temperature T=0.95.

The motivation is the data-rich era in astronomy and cosmology. As next-generation experiments deliver larger datasets, researchers can build more comprehensive models, but comparing those models becomes computationally harder. The learned harmonic mean is designed to work with accelerated sampling methods, including gradient-based samplers and simulation-based inference. Previous work applied related methods to cosmological parameter estimation in up to 159 parameter dimensions and to gravitational-wave model selection. The new architecture aims to make such evidence calculations more robust when posteriors are complex.

The tests are numerical and moderate in dimension. The paper notes that evidence estimation can be challenging even in moderate parameter spaces, and the authors plan to refine the architecture for very high-dimensional inference tasks such as imaging or field-level cosmological inference. They also plan to apply the method to real-world problems, including revisiting gravitational-wave astrophysics. The estimator has finite variance when the target density is contained within the posterior, and the temperature parameter enforces that containment.

The broader promise is simple. If evidence can be estimated from the same posterior samples that modern samplers already produce, then astronomers can compare models without forcing every problem into one sampling framework. When a model's best answers come in bunches, the scorekeeper needs to see all the bunches, and not mistake the empty space between them for part of the answer.