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 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 full distance-ladder dataset to correctly propagate uncertainties. We have developed a Bayesian hierarchical model that describes the full distance ladder, from nearby geometric anchors through Cepheids to Hubble-Flow SNe. This model does not rely on any distributions being Gaussian, allowing outliers to be modeled and obviating the need for arbitrary data cuts. Sampling from the $\sim$3000-parameter joint posterior using Hamiltonian Monte Carlo, we find $H_0$ = (72.72 $\pm$ 1.67) ${\rm km\,s^{-1}\,Mpc^{-1}}$ when applied to the outlier-cleaned Riess et al. (2016) data, and ($73.15 \pm 1.78$) ${\rm km\,s^{-1}\,Mpc^{-1}}$ with SN outliers reintroduced. Our high-fidelity sampling of the low-$H_0$ tail of the distance-ladder likelihood allows us to apply Bayesian model comparison to assess the evidence for deviation from $\Lambda$CDM. We set up this comparison to yield a lower limit on the odds of the underlying model being $\Lambda$CDM given the distance-ladder and Planck XIII (2016) CMB data. The odds against $\Lambda$CDM are at worst 10:1 or 7:1, depending on whether the SNe outliers are cut or modeled, or 60:1 if an approximation to the Planck Int. XLVI (2016) likelihood is used.