Inferring the eccentricity distribution

Standard maximum-likelihood estimators for binary-star and exoplanet eccentricities are biased high, in the sense that the estimated eccentricity tends to be larger than the true eccentricity. As with most non-trivial observables, a simple histogram of estimated eccentricities is not a good estimate of the true eccentricity distribution. Here we develop and test a hierarchical probabilistic method for performing the relevant meta-analysis, that is, inferring the true eccentricity distribution, taking as input the likelihood functions for the individual-star eccentricities, or samplings of the posterior probability distributions for the eccentricities (under a given, uninformative prior). The method is a simple implementation of a hierarchical Bayesian model; it can also be seen as a kind of heteroscedastic deconvolution. It can be applied to any quantity measured with finite precision–other orbital parameters, or indeed any astronomical measurements of any kind, including magnitudes, parallaxes, or photometric redshifts–so long as the measurements have been communicated as a likelihood function or a posterior sampling.

💡 Research Summary

The paper addresses a well‑known but often overlooked problem in the statistical analysis of binary‑star and exoplanet orbital eccentricities: standard maximum‑likelihood (ML) or least‑squares estimators are positively biased, especially when the true eccentricity is small and measurement uncertainties are comparable to the signal. Because eccentricity is bounded below by zero, noise pushes the estimated values upward, producing a distribution of estimated eccentricities that does not reflect the underlying population. A naïve histogram of these point estimates therefore misrepresents the true eccentricity distribution.

To correct this, the authors develop a hierarchical Bayesian framework that treats the eccentricities of individual systems as latent variables and the population distribution as a hyper‑parameterized prior. For each star i they define a likelihood function (L_i(e)=p(d_i|e)), where (d_i) denotes the observed data (radial‑velocity time series, transit timings, etc.) and (e) is the true eccentricity. The population distribution is denoted (\pi(e|\theta)), with (\theta) representing the hyper‑parameters (e.g., the shape parameters of a Beta distribution or the coefficients of a non‑parametric spline). The joint likelihood for the whole sample is

\