CosmoFit.stats

Statistical tools: priors, posteriors, the high-level Fitter, its MCMC sampling backend and consolidated result objects, model comparison (AIC/BIC/LRT), saved (HDF5) chains, CPL-specific posterior diagnostics (w(z) bands, w(z)=-1 crossing, Mahalanobis distance from LCDM), and expansion-history derived posteriors (transition redshift z_t, q0) that apply to every model.

Priors and the posterior

Prior distributions.

CosmoFit currently implements independent uniform (“top-hat”) priors, which is what essentially every MCMC-based cosmology analysis uses for the background parameters (H0, Omega_m, w0, wa, rd, …). More informative priors (Gaussian, etc.) can be added later following the same interface.

class CosmoFit.stats.priors.UniformPrior(names, bounds)

Bases: object

Independent uniform priors over a set of parameters.

Parameters:

• names (list[str]) – Parameter names, in the order matching the theta vector this prior will be evaluated on.

• bounds (dict[str, tuple[float, float]]) – Mapping of parameter name -> (lower, upper). Every name in names must have an entry.

log_prior(theta)

Return 0.0 if every entry of theta lies within its bounds, otherwise -inf.

Return type:

float

sample(n_samples, rng=None)

Draw n_samples independent samples from the prior.

Returns:

ndarray, shape (n_samples, ndim)

Parameters:

n_samples (int)

Return type:

ndarray

Posterior probability.

Glues a JointLikelihood, a prior, and a mutable set of cosmological parameters together into a single callable log_probability(theta) function of the kind samplers like emcee expect.

class CosmoFit.stats.posterior.LogPosterior(joint, cosmology, prior)

Bases: object

Callable log-posterior for a subset of “free” cosmological parameters, with the rest held fixed.

Parameters:

• joint (likelihoods.joint.JointLikelihood) – Joint likelihood (already built from one or more BaseLikelihood instances sharing a single cosmology object).

• cosmology (cosmology.core.base.Cosmology) – The (mutable) cosmology instance the likelihoods above were built with. Its .params will be updated in place at every evaluation.

• prior (statistics.priors.UniformPrior) – Prior over the free parameters. prior.names defines which entries of theta map to which parameter.

chi2(theta)

Convenience wrapper returning chi2(theta) directly – useful for scipy.optimize.minimize best-fit searches.

Some models have prior-bounded regions where the background is unphysical (e.g. E(z)^2 < 0 for a modified-gravity model at extreme coupling values) – log_likelihood already treats that as -inf log-probability during MCMC; mirror that here as +inf chi2 (worst possible fit) rather than letting scipy.optimize.minimize crash on a NaN/exception when its search steps into that region.

Two guards, and the first one exists because the second was not enough. A theta that is itself not finite is rejected before it reaches the cosmology at all: L-BFGS-B started at a point where chi2 is already inf computes a finite-difference gradient of inf - inf = nan, takes a nan search direction, and evaluates the objective at [nan, nan, nan]. Writing that into the parameters builds an interpolation table full of nan, and the interpolator raises rather than returning anything – from inside refresh(), which used to sit outside the try below.

So the try now covers _apply as well. A cosmology that cannot even be constructed is as excluded as one that fits badly, and neither should be able to crash a sampler that merely proposed it.

Return type:

float

Model comparison

Model comparison statistics.

Implements the standard chi2-based model selection tools used throughout cosmology: AIC, BIC, and the likelihood-ratio test for nested models (e.g. CPL vs its LCDM limit w0=-1, wa=0).

CosmoFit.stats.model_comparison.aic(chi2, k)

Akaike Information Criterion: chi2 + 2k.

Parameters:

• chi2 (float)

• k (int)

Return type:

float

CosmoFit.stats.model_comparison.bic(chi2, k, n_data)

Bayesian Information Criterion: chi2 + k * ln(n_data).

Parameters:

• chi2 (float)

• k (int)

• n_data (int)

Return type:

float

CosmoFit.stats.model_comparison.NESTED_CHI2_TOLERANCE = 0.001

How negative a nested delta_chi2 may be before it is treated as a real failure rather than convergence noise.

Not zero, and the reason is measured. When the general model’s best fit is the nested limit – LsCDM running off to z_dagger ~ 98, CPL sitting at w0 = -1, wa = 0 – the two optimizers stop at the same minimum by different routes and disagree at their own convergence tolerance. Observed across a 32-fit scan: -8e-06 to -5e-04, every one of them at the limit. Warning about those would train the reader to ignore the warning.

CosmoFit.stats.model_comparison.likelihood_ratio_test(chi2_null, k_null, chi2_alt, k_alt, tolerance=0.001)

Likelihood-ratio test between two nested models.

“null” is the simpler model (e.g. LCDM), “alt” is the more general model that reduces to it for a special parameter choice (e.g. CPL with w0=-1, wa=0).

Returns:

dict with keys – delta_chi2, delta_k, p_value, sigma

Parameters:

• chi2_null (float)

• k_null (int)

• chi2_alt (float)

• k_alt (int)

• tolerance (float)

Return type:

dict

Notes

A negative delta_chi2 is reported and warned about rather than passed through. Between nested models it cannot happen: the general model contains the simple one, so its minimum is at worst equal. Seeing one means an optimizer stopped somewhere that is not the minimum – typically in a second basin, which is a failure best_fit’s stall detection cannot see, since the run converged and reported success.

The formula would otherwise absorb it silently: chi2.sf(negative) == 1.0 and norm.isf(1.0) == -inf, so the result reads as “no evidence” instead of “this number is impossible”. sigma is therefore clamped to zero – a nested model cannot be evidence against the general one – while delta_chi2 is returned as measured, so the caller can see how bad it was. Refit with best_fit(restarts=...).

Only a shortfall larger than tolerance is warned about. A general model whose best fit is the nested limit reaches the same minimum by a different route and disagrees at the optimizer’s own convergence level – see NESTED_CHI2_TOLERANCE. sigma is clamped for any negative value regardless, since none of them is evidence.

CosmoFit.stats.model_comparison.compare_models(*, name_null, chi2_null, k_null, name_alt, chi2_alt, k_alt, n_data)

Full AIC/BIC/LRT comparison between two nested models, matching the “CPL vs LCDM” analysis block of the CPL_MCMC notebook.

Parameters:

• name_null (str)

• chi2_null (float)

• k_null (int)

• name_alt (str)

• chi2_alt (float)

• k_alt (int)

• n_data (int)

Return type:

dict

Bayesian evidence

Bayesian evidence by nested sampling.

Every model comparison in this library so far has been a delta chi2 at the best fit, wrapped in AIC, BIC or a likelihood-ratio test. Those are approximations, and the notebooks say so repeatedly – “these are delta chi2 at the best fit, not evidence ratios”. This computes the thing itself:

Z = Int L(theta) pi(theta) d theta

and ln Z is what a Bayes factor is a difference of.

Two reasons it is worth having rather than another information criterion.

It does not care about boundaries. Wilks’ theorem needs the null to be interior to the parameter space, and three of this library’s most interesting comparisons violate that: LsCDM reduces to LCDM only as z_dagger -> infinity, GEDE at an edge of its parameter, DGP not at all. The notebooks correctly refuse to quote a sigma in those cases. An evidence ratio is defined regardless.

It integrates rather than maximizes, so a parameter that improves the fit only in a sliver of its prior is penalized for it – the Occam factor that AIC and BIC only approximate, and that counts the volume rather than the parameter.

Which brings the warning that has to come with it.

Prior sensitivity

Z depends on the prior, and this library’s priors are uniform over DEFAULT_BOUNDS. Widening the range of a parameter the data does not constrain divides its evidence by roughly the widening factor, with the fit unchanged. That is not a defect of the method – it is the Occam factor doing its job – but it does mean a Bayes factor here is a statement about a model plus its priors, and comparing two models means having defended both. bayes_factor() reports the prior volumes alongside the ratio so that the dependence is visible rather than implicit.

How much it matters, measured on this library’s own case. LsCDM against LCDM on CC + DESI DR2 + low-z BAO + DES-SN5YR + Planck priors + BBN, changing nothing but the upper edge of the z_dagger prior:

z_dagger prior	ln B	label
(0.5, 10)	+1.158	positive
(0.5, 100)	+0.494	inconclusive

The same data, the same model, the same fit – and a different verdict, because above the sign-switch cliff the parameter is unconstrained and the model is charged for whatever room it was given. (The shift is 0.66 rather than the ln 10 = 2.30 a completely flat direction would give, because the likelihood does decline slowly towards the z_dagger -> infinity limit.)

For contrast, the frequentist number on the same fit is delta chi2 = 5.40, which Wilks would read as ~2 sigma and which does not apply here at all, LCDM being on the boundary. Neither statistic is wrong; they answer different questions, and only one of them has to be told what the prior is.

Optional dependency: pip install "cosmofit[evidence]".

CosmoFit.stats.evidence.interpret(ln_bayes_factor)

Kass & Raftery’s label for a log Bayes factor, on its absolute value – the sign says which model, the size says how much.

Parameters:

ln_bayes_factor (float)

Return type:

str

CosmoFit.stats.evidence.bayes_factor(alt, null)

Compare two completed nested-sampling runs.

Parameters:

• alt (NestedResult) – The more general model and the simpler one. Order only fixes the sign.

• null (NestedResult) – The more general model and the simpler one. Order only fixes the sign.

Returns:

dict – ln_B (positive favours alt), its uncertainty propagated from both runs, the Kass & Raftery label, and the two prior volumes – because a Bayes factor computed against uniform priors is a statement about those priors as much as about the models.

Return type:

dict

The nested-sampling run itself, and its result object.

Kept apart from stats.evidence, which only compares finished runs and needs no optional dependency to do it.

class CosmoFit.stats.nested.NestedResult(log_evidence, log_evidence_error, samples, free_params, prior_volume, n_live, n_evaluations, information=nan)

Bases: object

What a nested-sampling run produced.

Parameters:

• log_evidence (float)

• log_evidence_error (float)

• samples (ndarray)

• free_params (list)

• prior_volume (float)

• n_live (int)

• n_evaluations (int)

• information (float)

log_evidence, log_evidence_error

ln Z and its uncertainty, which for nested sampling is an honest sampling error rather than a scatter estimate.

Type:

float

samples

Posterior samples, (n, ndim), already resampled to equal weight – so they can be summarized exactly like MCMC output.

Type:

ndarray

free_params

Type:

list of str

prior_volume

Volume of the uniform prior box the evidence was computed against. Reported because ln Z moves with it, and a Bayes factor that does not say which box it used is not a reproducible number.

Type:

float

n_live

Type:

int

n_evaluations

Type:

int

summary(percentiles=(2.5, 16, 50, 84, 97.5))

Per-parameter posterior percentiles, in the same shape stats.fitter.Fitter.summary() returns.

Return type:

dict

CosmoFit.stats.nested.run_nested(logpost, prior, free_params, n_live=500, dlogz=0.05, seed=42, progress=True, **dynesty_kwargs)

Integrate the posterior with dynesty.

Parameters:

• logpost (LogPosterior) – Only its log_likelihood is used. The prior enters through the unit-cube transform instead, which is what makes the integral an evidence rather than an unnormalized posterior.

• prior (UniformPrior)

• free_params (list of str)

• n_live (int, optional) – Live points. The evidence error scales roughly as sqrt(information / n_live), so this is the accuracy knob.

• dlogz (float, optional) – Stopping criterion: the estimated remaining evidence.

• seed (int, optional)

• progress (bool, optional)

Returns:

NestedResult

Return type:

NestedResult

Tension

How much two measurements disagree.

This library quotes tensions constantly – 4.1 sigma in the Hubble constant, 2.9 sigma in S8 – and until now computed every one of them by hand as

sigma = |a - b| / sqrt(sigma_a^2 + sigma_b^2)

That formula is right, and it assumes three things that are worth being asked about each time rather than once: that both posteriors are Gaussian, that they are one-dimensional, and that they are independent. The functions here are the same arithmetic where those hold, and the alternatives where they do not.

Four of them, in increasing order of what they need:

gaussian_tension()

Two numbers with error bars. The formula above, named, with its assumptions written down.

sample_tension()

Two sets of posterior samples of the same parameter. Makes no Gaussian assumption – it builds the distribution of the difference and asks how much of it lies further from zero than zero itself. The right tool for the skewed posteriors this library keeps producing.

gaussian_tension_nd()

Two multi-dimensional posteriors, summarized by mean and covariance. A tension in a plane is not the larger of its two projections, and can be much bigger than either.

suspiciousness()

Two evidences plus their joint. This is the one that does not care about the prior – which matters, because stats.evidence has to warn that a Bayes factor does.

References

Handley & Lemos (2019), Phys. Rev. D 100, 043504, arXiv:1902.04029 (suspiciousness).

Raveri & Hu (2019), Phys. Rev. D 99, 043506, arXiv:1806.04649 (the parameter-difference family this borrows from).

CosmoFit.stats.tension.gaussian_tension(a, sigma_a, b, sigma_b)

Tension between two independent Gaussian measurements of one parameter.

Parameters:

• a (float)

• sigma_a (float)

• b (float)

• sigma_b (float)

Returns:

dict – difference, combined_sigma, n_sigma, p_value.

Return type:

dict

Notes

Assumes both posteriors are Gaussian, one-dimensional and independent. The first of those is the one that usually fails here – see sample_tension(), which does not need it.

CosmoFit.stats.tension.sample_tension(samples_a, samples_b, n_pairs=200000, seed=0, bins=200)

Tension between two posteriors of the same parameter, from samples, with no Gaussian assumption.

The parameter-difference construction: if the two are independent, the distribution of a - b is what you get by pairing their samples at random. Perfect agreement puts that distribution’s peak at zero; a tension pushes zero into its tail. The quoted probability is the fraction of the difference distribution at a higher density than zero – so it works for a skewed or double-peaked difference, where “how many standard deviations from zero” would not.

Parameters:

• samples_a (array_like) – Posterior samples, equally weighted. Need not be the same length.

• samples_b (array_like) – Posterior samples, equally weighted. Need not be the same length.

• n_pairs (int, optional) – Random pairs drawn to build the difference distribution.

• seed (int, optional)

• bins (int, optional) – Histogram resolution for the density comparison.

Returns:

dict – n_sigma, p_value, median_difference, and the difference samples themselves as difference.

Return type:

dict

CosmoFit.stats.tension.gaussian_tension_nd(mean_a, cov_a, mean_b, cov_b)

Tension between two multi-dimensional Gaussian posteriors.

chi2 = (mu_a - mu_b)^T (C_a + C_b)^-1 (mu_a - mu_b)

read against a chi-square with n degrees of freedom.

Worth having separately because a tension in a plane is not the larger of its two one-dimensional projections: two posteriors can overlap in every parameter separately and still be far apart jointly, if their degeneracy directions differ.

Returns:

dict – chi2, dof, p_value, n_sigma.

Return type:

dict

CosmoFit.stats.tension.suspiciousness(joint, first, second)

Tension from evidences, with the prior dependence divided out.

The evidence ratio

ln R = ln Z_AB - ln Z_A - ln Z_B

measures whether two datasets prefer to be described together, but it moves with the prior volume – widen a prior and R rises, with nothing about the data changed. stats.evidence documents that at length, and measures it.

Handley & Lemos’ suspiciousness removes it by subtracting the information the data gained,

ln S = ln R - ln I, ln I = D_A + D_B - D_AB

where D is each run’s Kullback-Leibler divergence from prior to posterior – which carries the same prior dependence and cancels it. What is left is a statement about the data.

Parameters:

• joint (stats.nested.NestedResult) – Nested-sampling runs of the two datasets together and separately, over the same parameters and priors.

• first (stats.nested.NestedResult) – Nested-sampling runs of the two datasets together and separately, over the same parameters and priors.

• second (stats.nested.NestedResult) – Nested-sampling runs of the two datasets together and separately, over the same parameters and priors.

Returns:

dict – ln_S, ln_R, ln_I, d (the effective number of constrained parameters), chi2, p_value, n_sigma.

Return type:

dict

Notes

The degrees of freedom use d = d_A + d_B - d_AB with d = 2 * (<ln L> - ln L_max)-style Bayesian model dimensionality; here it is approximated by the number of parameters, which is exact when every one of them is constrained by both datasets and an overestimate otherwise. Pass runs over the same parameter set and this is the standard result.

Derived quantities

Posteriors for quantities derived from the expansion history.

stats.cpl_diagnostics covers quantities that are closed-form functions of CPL’s w0/wa (the w(z)=-1 crossing, the Mahalanobis distance from LCDM). The quantities here instead need the model’s actual E(z)/dE/dz, so they are computed by pushing posterior samples back through the cosmology – which means they work for every model in the library, not just CPL.

z_t   the acceleration transition redshift, where the
      deceleration parameter q(z) changes sign
q0    the present-day deceleration parameter, q(z=0)
r_d   the BAO sound horizon at the drag epoch, when it is
      computed from the physical densities rather than fitted

The first two are standard reported numbers in dark-energy papers, and both are what the fitter.plots.deceleration() figure shows graphically. The third exists because compute_rd=True (cosmology.calculators.sound_horizon) turns r_d from a sampled parameter into a derived one – so it disappears from fitter.summary(), and the only honest way to quote it with an error bar is to push every posterior sample back through the sound-horizon integral, exactly as z_t and q0 are handled.

Why not just evaluate at the best fit

Because q(z) is nonlinear in the parameters, z_t evaluated at the posterior-median parameters is not the median of the z_t posterior, and it carries no uncertainty. Quoting a derived quantity with an error bar means mapping every posterior sample through the transformation and taking percentiles of the result, which is what these functions do.

Example

>>> from CosmoFit.stats import derived
>>> fit.run_mcmc(...)
>>> z_t = derived.transition_redshift(fit)
>>> derived.summarize(z_t)
{'median': 0.73..., 'plus': 0.04..., 'minus': 0.04..., ...}

CosmoFit.stats.derived.q_of_z(fit, z, burnin=None, max_samples=5000)

Deceleration parameter q(z) for every (thinned) posterior sample.

Parameters:

• fit (stats.fitter.Fitter) – A fitter with a completed run_mcmc().

• z (array_like) – Redshifts to evaluate q at.

• burnin (int, optional) – Steps to discard. Defaults to the fitter’s own burnin.

• max_samples (int or None, optional) – Cap on the number of posterior samples used (see _MAX_SAMPLES). None uses every sample.

Returns:

ndarray, shape (n_samples, len(z))

CosmoFit.stats.derived.transition_redshift(fit, burnin=None, max_samples=5000, z_max=3.0, n_grid=601)

Posterior of the acceleration transition redshift z_t, where the deceleration parameter changes sign (q(z_t) = 0) – the epoch the universe switched from decelerating to accelerating.

The root is bracketed on a uniform grid over [0, z_max] and then refined by linear interpolation across the bracketing interval. With the default 601-point grid this agrees with the closed-form flat-LCDM result z_t = (2 Omega_Lambda / Omega_m)^(1/3) - 1 to 4e-6 in z_t across a real posterior – four orders of magnitude below the ~0.03 width of the posterior itself.

Parameters:

• fit – See q_of_z().

• burnin – See q_of_z().

• max_samples – See q_of_z().

• z_max (float, optional) – Upper end of the search range.

• n_grid (int, optional) – Grid resolution used to bracket the sign change.

Returns:

ndarray – One z_t per posterior sample. Samples with no sign change in [0, z_max] – an expansion history that never transitions in range – are nan; summarize() drops them, and their count is reported as n_undefined.

CosmoFit.stats.derived.deceleration_today(fit, burnin=None, max_samples=5000)

Posterior of the present-day deceleration parameter, q(z=0).

Returns:

ndarray, one value per posterior sample.

CosmoFit.stats.derived.sound_horizon(fit, burnin=None, max_samples=5000)

Posterior of the BAO sound horizon at the drag epoch, r_d [Mpc], computed from each sample’s physical densities.

Useful in two different situations, and it is worth being clear about which one you are in:

• With compute_rd=True, r_d is what the fit actually used, and it is derived rather than sampled – so this is the only place it appears with an error bar.

• Without it, r_d was a free parameter and this returns what the early-universe physics would have predicted for the same densities. Comparing the two is a real consistency test: a fitted r_d that disagrees with the computed one is the standard signature of new physics before recombination, and it is one of the main ways the Hubble tension is diagnosed.

Parameters:

• fit (stats.fitter.Fitter) – A fitter with a completed run_mcmc().

• burnin (int, optional) – Steps to discard. Defaults to the fitter’s own burnin.

• max_samples (int or None, optional) – Cap on the number of posterior samples used.

Returns:

ndarray, one value per posterior sample.

Notes

r_d depends only on Omega_m, Omega_b, H0, N_eff and m_nu, so if none of those is free this returns a constant array – correctly, since nothing in the posterior can move it.

CosmoFit.stats.derived.summarize(values, percentiles=(2.5, 16, 50, 84, 97.5))

Median and credible intervals of a derived-quantity posterior, in the same median/plus/minus shape summary() uses.

Non-finite entries (e.g. samples with no transition redshift, see transition_redshift()) are dropped and counted separately rather than silently poisoning the percentiles.

Returns:

• dict with keys median, plus, minus, mean, std,

• lower95, upper95, n, n_undefined.

Return type:

dict

CPL diagnostics

CPL-specific posterior diagnostics.

Utility functions that turn MCMC posterior samples of the CPL equation-of-state parameters (w0, wa) into the derived quantities used to characterize deviations from LCDM:

• w(z) posterior bands

• the “w(z) = -1 crossing” redshift distribution

• the crossing direction (quintessence -> phantom or the reverse)

• which dark-energy region of the (w0, wa) plane a point falls in – phantom, quintessence, quintom-A or quintom-B – and how the posterior splits between them

• the Mahalanobis distance of the LCDM point (w0, wa) = (-1, 0) from the 2D posterior

These reproduce the final analysis cells of the CPL_MCMC notebook.

CosmoFit.stats.cpl_diagnostics.wz(w0_samples, wa_samples, z)

CPL equation of state w(z) = w0 + wa * z / (1+z), evaluated for every posterior sample at every redshift in z.

Returns:

ndarray, shape (n_samples, len(z))

CosmoFit.stats.cpl_diagnostics.wz_posterior_bands(w0_samples, wa_samples, z, quantiles=(2.5, 16, 50, 84, 97.5))

Percentile bands of w(z) over the posterior, at each redshift in z.

Returns:

• dict mapping each quantile (as given, e.g. 16) to an

• ndarray of length len(z).

CosmoFit.stats.cpl_diagnostics.crossing_redshift(w0_samples, wa_samples, z_min=0.0, z_max=2.5)

Posterior distribution of the redshift at which w(z) = -1,

z_cross = -(1 + w0) / (1 + w0 + wa)

Only samples with a physical crossing (finite, and inside [z_min, z_max]) are kept.

Returns:

• z_cross (ndarray) – Crossing redshifts for the physically-valid samples.

• fraction (float) – Fraction of all posterior samples with a physical crossing in range.

CosmoFit.stats.cpl_diagnostics.crossing_direction(w0_samples, wa_samples, z_min=0.0, z_max=2.5, z_ref=2.5)

Among the samples with a physical w(z)=-1 crossing, what fraction cross from quintessence-like (w>-1) at z=0 to phantom-like (w<-1) at high z, versus the reverse?

Returns:

• dict with keys ‘quintessence_to_phantom’,

• ’phantom_to_quintessence’ (fractions of the crossing

• sub-sample), and ‘n_crossing’.

CosmoFit.stats.cpl_diagnostics.REGIONS = ('phantom', 'quintessence', 'quintom-a', 'quintom-b')

The four regions the (w0, wa) plane splits into, in the order they are usually listed. See classify_region().

CosmoFit.stats.cpl_diagnostics.classify_region(w0, wa)

Which dark-energy region a (w0, wa) point belongs to.

A CPL equation of state runs from w(0) = w0 today to w(z -> inf) = w0 + wa in the distant past, monotonically. Whether it sits above or below the cosmological-constant value w = -1 at those two ends is what the standard classification is:

"phantom"

w < -1 at both ends – below -1 for all time.

"quintessence"

w > -1 at both ends; reachable by an ordinary canonical scalar field.

"quintom-a"

w > -1 in the past, w < -1 today: the equation of state crosses -1 going forward in time.

"quintom-b"

the reverse crossing, phantom in the past and quintessence-like today.

The crossing itself is the physically loaded part: a single canonical (or single phantom) scalar field cannot cross w = -1 at all, so a posterior sitting in either quintom region calls for something else – two fields, a non-canonical kinetic term, or modified gravity.

Parameters:

• w0 (float or array_like) – Scalars, or arrays of posterior samples.

• wa (float or array_like) – Scalars, or arrays of posterior samples.

Returns:

str, or ndarray of str – A scalar in, a scalar out; arrays in, an array of the same shape out.

CosmoFit.stats.cpl_diagnostics.region_fractions(w0_samples, wa_samples)

Posterior probability of each dark-energy region: the fraction of samples falling in each of REGIONS.

This is the number behind the picture that fit.plots.w0_wa_plane() draws – “the contours sit in the quintom-B region” stated as a probability rather than by eye. Because the regions tile the plane, the four fractions sum to 1.

Returns:

dict[str, float] – One entry per region, in REGIONS order.

Return type:

dict

CosmoFit.stats.cpl_diagnostics.mahalanobis_from_lcdm(w0_samples, wa_samples, lcdm_point=(-1.0, 0.0))

How far the LCDM point (w0, wa) = (-1, 0) sits from the 2D (w0, wa) posterior.

Returns both the raw Mahalanobis distance and – separately – the significance, because the two are not the same number in 2D and conflating them overstates the tension with LCDM.

distance (D) is the Mahalanobis distance, sqrt((x-mu)^T C^-1 (x-mu)). It is not a number of sigma. For a 2D Gaussian, D^2 follows a chi-square distribution with 2 degrees of freedom, not 1, so the probability enclosed at a given D is smaller than the familiar 1D intuition suggests:

D = 1.515  encloses 68.27%  (the "1 sigma" probability)
D = 2.486  encloses 95.45%  (the "2 sigma" probability)
D = 3.439  encloses 99.73%  (the "3 sigma" probability)

So a D of 2.20 is not 2.2 sigma – it corresponds to 91.1% exclusion, i.e. 1.70 sigma in the usual one-dimensional two-tailed sense. sigma below applies that conversion; quote it, not distance, when reporting a deviation from LCDM.

Returns:

dict –

With keys:

mean, covariance   the posterior mean and covariance used
distance_squared   D^2
distance           D, the Mahalanobis distance (NOT sigma)
p_value            P(chi2_2 > D^2): the probability of
                   being at least this far out under the
                   posterior
confidence_level   1 - p_value, i.e. the confidence at
                   which LCDM is excluded
sigma              the equivalent one-dimensional
                   two-tailed significance -- the number
                   to report

Notes

This is a Gaussian approximation to the posterior: it uses only the mean and covariance, so a strongly non-Gaussian (e.g. banana-shaped) w0-wa contour is not fully captured. Compare against the fraction of samples further out than the LCDM point if that matters.