Estimates of the Hubble constant, H 0 , from the local distance ladder and the cosmic microwave background (CMB) are discrepant at the ∼3-σ level, indicating a potential issue with the standard ΛCDM cosmology. Interpreting this tension correctly requires a model comparison calculation which depends on not only the traditional 'n-σ' mismatch but also the tails of the likelihoods. Determining the form of the tails of the local H 0 likelihood is impossible with the standard Gaussian or least-squares approximation, as it requires using non-Gaussian distributions to faithfully represent anchor likelihoods and model outliers in the Cepheid and supernova (SN) populations, and simultaneous fitting of the complete distance-ladder dataset to ensure correct propagation of uncertainties. We have hence developed a Bayesian hierarchical model (BHM) that describes the full distance ladder, from nearby geometric-distance anchors through Cepheids to SNe in the Hubble flow. This model does not rely on any of the underlying distributions being Gaussian, allowing outliers to be modeled and obviating the need for any arbitrary data cuts. Sampling from the full ∼3000-parameter joint posterior distribution using Hamiltonian Monte Carlo and marginalizing over the nuisance parameters (i.e., everything bar H 0 ), we find H 0 = (72.72 ± 1.67) km s −1 Mpc −1 when applied to the outlier-cleaned Riess et al. (2016) data, and (73.15 ± 1.78) km s −1 Mpc −1 with SN outliers reintroduced (the pre-cut Cepheid dataset is not available). Our highfidelity sampling of the low-H 0 tail of the distance-ladder posterior allows us to apply Bayesian model comparison to assess the evidence for deviation from ΛCDM. We set up this comparison to yield a lower limit on the odds of the underlying model being ΛCDM given the distance-ladder and Planck Collaboration (2016b) CMB data. The odds against ΛCDM are at worst 10:1 when considering the outlier-free distanceladder data, or 7:1 when the SNe outliers are included and modeled, both considerably less dramatic than naïvely implied by the 2.8-σ discrepancy. These odds become ∼60:1 when an approximation to the more-discrepant Planck Collaboration (2016c) likelihood is included. The code used in this analysis is made publicly available at https://github.com/sfeeney/hh0.

It is well known in astronomy that propagating non-Gaussian prediction uncertainty in photometric redshift estimates is key to reducing bias in downstream cosmological analyses. Similarly, likelihoodfree inference approaches, which are beginning to emerge as a tool for cosmological analysis, require the full uncertainty landscape of the parameters of interest given observed data. However, most machine learning (ML) based methods with open-source software target point prediction or classification, and hence fall short in quantifying uncertainty in complex regression and parameter inference settings such as the applications mentioned above. As an alternative to methods that focus on predicting the response (or parameters) y from features x, we provide nonparametric conditional density estimation (CDE) tools for approximating and validating the entire probability density p(y | x) given training data for x and y. This density approach offers a more nuanced accounting of uncertainty in situations with, e.g., nonstandard error distributions and multimodal or heteroskedastic response variables that are often present in astronomical data sets. As there is no one-size-fits-all CDE method, and the ultimate choice of model depends on the application and the training sample size, the goal of this work is to provide a comprehensive range of statistical tools and open-source software for nonparametric CDE and method assessment which can accommodate different types of settings -involving, e.g., mixedtype input from multiple sources, functional data, and image covariates -and which in addition can easily be fit to the problem at hand. Specifically, we introduce CDE software packages in Python and R based on four ML prediction methods adapted and optimized for CDE: NNKCDE, RFCDE, FlexCode, and DeepCDE. Furthermore, we present the cdetools package, which includes functions for computing a CDE loss function for model selection and tuning of parameters, together with diagnostic functions for computing posterior quantiles and coverage probabilities. We provide sample code in Python and R as well as examples of applications to photometric redshift estimation and likelihood-free cosmology via CDE.

Many areas of science make extensive use of computer simulators that implicitly encode likelihood functions of complex systems. Classical statistical methods are poorly suited for these so-called likelihood-free inference (LFI) settings, outside the asymptotic and lowdimensional regimes. Although new machine learning methods, such as normalizing flows, have revolutionized the sample efficiency and capacity of LFI methods, it remains an open question whether they produce reliable measures of uncertainty.This paper presents a statistical framework for LFI that unifies classical statistics with modern machine learning to: (1) efficiently construct frequentist confidence sets and hypothesis tests with finite-sample guarantees of nominal coverage (type I error control) and power; (2) provide practical diagnostics for assessing empirical coverage over the entire parameter space. We refer to our framework as likelihood-free frequentist inference (LF2I). Any method that estimates a test statistic, like the likelihood ratio, can be plugged into our framework to create valid confidence sets and compute diagnostics, without costly Monte Carlo samples at fixed parameter settings. In this work, we specifically study the power of two test statistics (ACORE and BFF), which, respectively, maximize versus integrate an odds function over the parameter space. Our study offers multifaceted perspectives on the challenges in LF2I.

No abstract