Cosmology-marginalized approaches in Bayesian model comparison: the neutrino mass as a case study

TL;DR

This paper introduces a new Bayesian method to derive robust cosmological parameter constraints that account for model dependence, with a case study on neutrino mass bounds relevant to dark matter and future neutrino experiments.

Contribution

It presents a novel approach to marginalize over cosmological models in Bayesian analysis, improving the robustness of parameter constraints compared to traditional methods.

Findings

Preferred cosmological models are less favored when considering multiple models.

The method provides model-independent bounds on neutrino mass.

Application to neutrino mass yields more reliable limits for dark matter and experimental planning.

Abstract

We propose here a \emph{novel} method which singles out the \emph{a priori} unavoidable dependence on the underlying cosmological model when extracting parameter constraints, providing robust limits which only depend on the considered dataset. Interestingly, when dealing with several possible cosmologies and interpreting the Bayesian preference in terms of the Gaussian statistical evidence, the preferred model is much less favored than when only two cases are compared. As a working example, we apply our approach to the cosmological neutrino mass bounds, which play a fundamental role not only in establishing the contribution of relic neutrinos to the dark matter of the Universe, but also in the planning of future experimental searches of the neutrino character and of the neutrino mass ordering.

Tables

Table 1. Table 1: 95 % percent 95 95\%  CL upper limits on Σ ​ m ν Σ subscript 𝑚 𝜈 \Sigma m_{\nu} and Bayes factors in the different cosmological scenarios. Results are obtained either adopting Planck 2015 CMB temperature, low- ℓ ℓ \ell polarization and lensing data  Adam  et al. ( 2016 ); Ade  et al. ( 2016a , b ) plus BAO measurements  Anderson  et al. ( 2014 ); Beutler  et al. ( 2011 ); Ross  et al. ( 2015 ) (second and third columns) or the same data combination plus high-multipole CMB polarization measurements from the Planck 2015 data release (fourth and fifth columns). The different rows depict the bounds in different extensions of the Λ Λ \Lambda CDM model, while the last two rows illustrate, respectively, the model-marginalized 95 % percent 95 95\%  CL limit obtained via Eq. ( 7 ) and the posterior probability of the example model ℳ 0 subscript ℳ 0 \mathcal{M}_{0} (the preferred one, that is always the Λ Λ \Lambda CDM+ Σ ​ m ν Σ subscript 𝑚 𝜈 \Sigma m_{\nu}  scenario), see Eq. ( 5 ).

	CMB+lens+BAO		CMB+pol+lens+BAO
model	$\ln B_{i0}$	$\Sigma m_{\nu }$ [eV]	$\ln B_{i0}$	$\Sigma m_{\nu }$ [eV]
base=$\Lambda$CDM+$\Sigma m_{\nu }$	$0.0$		$0.0$
base+$A_{\mathrm{lens}}$	$-2.6$		$-2.4$
base+$N_{\mathrm{eff}}$	$-1.5$		$-2.3$
base+$w$	$-1.4$		$-0.1$
marginalized	$-$		$-$
$p_0$	0.65		0.48

Equations

$$p(\theta |d,{\mathcal{M}}_i)=\frac{\pi (\theta |{\mathcal{M}}_i)\mathcal{L}(\theta )}{Z_i},$$

$$p_i\equiv p({\mathcal{M}}_i|d)=\frac{{\pi }_i{Z}_i}{{\sum }_j{\pi }_j{Z}_j},$$

$$B_{ij}=Z_i/Z_j.$$

$$\frac{p_i}{p_j}=B_{ij}\frac{{\pi }_i}{{\pi }_j}.$$

$$p_0=\frac{Z_0}{{\sum }_i{Z}_i}={\left(\sum _i^N{B}_{i0}\right)}^{-1},$$

$$p(\theta |d)=\sum _i^{N}p(\theta |d,{\mathcal{M}}_i)p_i,$$

$$p(\theta |d)=\sum _i^{N}p(\theta |d,{\mathcal{M}}_i)Z_i/\sum _j^N{Z}_j.$$

Full text

Introduction— Bayesian parameter inference has been extremely successful in cosmology and astroparticle physics in the past two decades. This statistical technique is more powerful and adequate than traditional tools when dealing with large and complex data sets and with the impossibility to obtain different realizations of the object to study, our universe. In addition, the Bayesian probability theory has also been extensively exploited for model comparison purposes, offering not only the possibility of predicting but also of optimizing the most adequate theoretical frameworks to fit the cosmological observations, see e.g. Trotta (2008). However, despite the major accomplishments achieved by Bayesian parameter inference, both the role of parameterizations/priors and the possibility of different fiducial cosmologies (or models) may led to divergent predictions. The former have caused controversial arguments in the literature, particularly when extracting cosmological bounds on neutrino masses and on their ordering Simpson et al. (2017); Schwetz et al. (2017); Gariazzo et al. (2018); Long et al. (2018); Heavens and Sellentin (2018); Handley and Millea (2018).

In this letter, we shall focus on the potential that Bayesian model comparison techniques offer for computing model-marginalized cosmological parameter limits, avoiding the biases due to the fiducial cosmology. We propose here a simple method to compute such solid and robust model-marginalized constraints.

In order to demonstrate the validity and robustness of this method, we shall illustrate a particular case and consider the sum of the neutrino mass $\Sigma m_{\nu}$ (see Refs. Lesgourgues et al. (2013); Lattanzi and Gerbino (2018); De Salas et al. (2018) for its key signatures on cosmology). Focusing exclusively on bounds from Cosmic Microwave Background (CMB) measurements, the final analyses from the Planck satellite set a $95\%$ CL limit of $\Sigma m_{\nu}<0.24$ eV Aghanim et al. (2018) after considering CMB temperature, polarization and lensing at all scales. Late-time observations of the large scale structure in the universe by means of the Baryon Acoustic Oscillation (BAO) method sharpen the limit above, as they help enormously in removing the degeneracies present in CMB data at the background level. Once BAO information is combined with Planck measurements, the limit is tightened to $\Sigma m_{\nu}<0.12$ eV at $95\%$ CL Aghanim et al. (2018), or even down to $\Sigma m_{\nu}<0.11$ eV when also considering Supernovae Ia luminosity distances. One obvious question is: how reliable and stable the cosmological neutrino mass limits quoted above are?

Even if not relying on the combination of potentially inconsistent data sets, for which the neutrino mass bounds become tighter 111As an example, adding the prior on the Hubble constant provided by Ref. Riess et al. (2018) one would get $\Sigma m_{\nu}<0.0970$ eV at $95\%$ CL Aghanim et al. (2018))., all of the aforementioned limits are based on the most economical $\Lambda$CDM scenario, which also leads to the tightest constraints on the neutrino mass. Surely, the bounds on $\Sigma m_{\nu}$ change when (a) new physics is added in the neutrino sector (for instance, changing the effective number of relativistic degrees of freedom, $N_{\rm{eff}}$ Hamann et al. (2007a, 2010); Giusarma et al. (2011); Hamann et al. (2011); Giusarma et al. (2013); Riemer-Sorensen et al. (2013); Archidiacono et al. (2013a); Di Valentino et al. (2013); Archidiacono et al. (2013b) or adding non-standard interactions Beacom et al. (2004); Bell et al. (2006); Hannestad (2005); Fardon et al. (2004); Afshordi et al. (2005); Brookfield et al. (2006a, b); Bjaelde et al. (2008); Mota et al. (2008); Ichiki and Keum (2008); Boehm et al. (2012); Archidiacono and Hannestad (2014); Dvali and Funcke (2016); Di Valentino et al. (2018)), (b) new physics appears in the early or late-time accelerating periods in the universe Hamann et al. (2007b); Joudaki (2013); Archidiacono et al. (2013b); de Putter et al. (2014); Di Valentino et al. (2016); Canac et al. (2016); Gerbino et al. (2017); Di Valentino et al. (2016); Gerbino et al. (2017); Huterer and Linder (2007); Baldi et al. (2014); Hu et al. (2015); Shim et al. (2014); Barreira et al. (2014); Bellomo et al. (2017); Peirone et al. (2018); Renk et al. (2017); Dirian (2017); Gavela et al. (2009); Reid et al. (2010); Lopez Honorez and Mena (2010); La Vacca et al. (2009); Guo et al. (2018) and/or, in general, (c) phenomenologically extended scenarios are considered Di Valentino et al. (2015). While one would naively expect that the neutrino mass limits within these more general cosmologies will always be relaxed, this has been shown to not to be the case for physical dark energy models Vagnozzi et al. (2018); Choudhury and Choubey (2018), for which the neutrino mass bounds get tighter than those obtained in the $\Lambda$CDM framework. It is therefore clear that one can artificially tune the cosmological neutrino mass limits in an optimistic or in a pessimistic manner.

These a priori harmless uncertainties translate into very serious dilemmas for neutrino particle physics searches. The near and far future neutrinoless double beta decay roadmap provides a very important example. It seems therefore mandatory to build a method to extract model-independent cosmological neutrino mass bounds. It is among our major goals to apply our novel model-marginalized method to $\Sigma m_{\nu}$ when studying a number of possible cosmological scenarios, i.e. the minimal $\Lambda$CDM universe with massive neutrinos and its extensions. Adopting Planck 2015 data 222Although the final data analyses have already been presented by the Planck collaboration, the data and likelihood codes are not public at the time of writing., the tightest bound we obtain within a $\Lambda$CDM universe is $\Sigma m_{\nu}<0.23$ eV at $95\%$ CL, which relaxes to $\Sigma m_{\nu}<0.35$ eV when the uncertainty on the cosmological model is taken into account using our model-marginalization method.

At the same time, we can use Bayesian tools in order to compare the models we are studying and obtain which one is preferred by data. Noticeably, even if the best scenario is strongly favoured over its competitors when comparing pairs of models with a Bayes factor analysis, its global statistical evidence falls abruptly when all the models are considered simultaneously, making this preferred model less likely. In the scenarios explored here, this will imply that the weak-to-moderate Bayesian preference for the minimal $\Lambda$CDM+$\Sigma m_{\nu}$ model, which arises when it is compared with each of its extensions individually, will not correspond to a global $1\sigma$ level strength when considering the entire ensemble of extended scenarios.

Bayesian statistics— The Bayes theorem, which represents the foundation of Bayesian statistics, reads:

[TABLE]

where $\pi(\theta|\mathcal{M}_{i})$ and $p(\theta|d,\mathcal{M}_{i})$ are the prior and posterior probabilities for the parameters $\theta$ within a model $\mathcal{M}_{i}$, $\mathcal{L}(\theta)$ is the likelihood as a function of the parameters $\theta$, given the data $d$ and the model $\mathcal{M}_{i}$, and $Z_{i}=\int d\theta\,\pi(\theta|\mathcal{M}_{i})\,\mathcal{L}(\theta)$ is the Bayesian evidence of $\mathcal{M}_{i}$ Trotta (2008). The Bayes theorem can also be written in a slightly different form to obtain the model posterior probability Handley et al. (2015a):

[TABLE]

where $\pi_{i}\equiv\pi(\mathcal{M}_{i})$ refers to the model prior probability. In the Bayesian model comparison framework, the so-called Bayes factor provides a measure of whether the data have increased or decreased the odds of model $\mathcal{M}_{i}$ relative to a second model $\mathcal{M}_{j}$:

[TABLE]

The Bayes factor enters the definition of the posterior probability ratio between two models, which indicates how much one of the two is preferred over the other, after using the information provided by data:

[TABLE]

If the two models are equivalent according to our initial knowledge, i.e. the model priors are the same, the final preference driven by data is determined by the Bayes factor. In terms of posterior odds, the preference for the favored model is $B_{ij}:1$, if $\mathcal{M}_{i}$ is preferred over $\mathcal{M}_{j}$. Adopting the commonly exploited Jeffreys’ scale Jeffreys (1961), the strength of the posterior odds can be ranked as inconclusive ($|\ln B_{ij}|<1$), weak ($1<|\ln B_{ij}|<2.5$), moderate ($2.5<|\ln B_{ij}|<5$), or strong ($|\ln B_{ij}|>5$). Very importantly, this arises from the fact that when comparing two mutually exclusive models, the mentioned ranks correspond roughly to what is usually indicated as $\lesssim 1\sigma$ (inconclusive) to $\gtrsim 3\sigma$ (strong) level when considering a Gaussian variable.

Using Eq. (2), and selecting one among the available models, labelled $\mathcal{M}_{0}$ without loss of generality, one can write, provided all priors are identical for all models:

[TABLE]

where we have used the definition of the Bayes factor. Notice that the posterior probability of the selected model $\mathcal{M}_{0}$ depends on the Bayes factors with respect to all the possible models. For each data combination we will choose $\mathcal{M}_{0}$ to be the preferred model. In practice, this is the one that has more influence on the model-marginalized posterior. Since the model $\mathcal{M}_{0}$ is the preferred one, we will always have $B_{i0}=Z_{i}/Z_{0}<1$ (or $\ln B_{i0}<0$) for $i\neq 0$.

Assuming that (i) more than two models are possible; and (ii) all the models have the same prior probabilities, then Eq. (5) implies that the posterior probability of the preferred model is smaller than what the single Bayes factors would suggest in a one-to-one comparison. For example, if $N=8$ and all the Bayes factors are $|\ln B_{i0}|\simeq 5$ for $i\neq 0$, thus indicating apparently strong results according to the usually adopted Jeffreys’ scale, the posterior probability of $\mathcal{M}_{0}$ is $p_{0}\simeq 0.955$, which would indicate a mild $2\sigma$ significance for a Gaussian measure. In the same way, having $N=7$ and $|\ln B_{i0}|\simeq 2.5$ for $i\neq 0$, which usually indicates a weak preference, would give $p_{0}\simeq 0.67$, which would correspond to less than $1\sigma$ preference for $\mathcal{M}_{0}$.

The tools of model comparison also allow us to compute a model-marginalized posterior distribution for the parameter $\theta$, taking into account the posterior probability of each model $\mathcal{M}_{i}$ resulting from the data $d$ Trotta (2008):

[TABLE]

where the posterior probabilities of $\theta$ within each model $\mathcal{M}_{i}$ are weighted according to the model posterior probabilities $p_{i}$. These can be written using Eq. (2) to obtain the fundamental formula

[TABLE]

This is the expression that we will use to obtain model-marginalized limits in the following, under the assumption that all the models have the same priors.

Some final comments are due. To obtain the most robust model-marginalized estimate one should in principle consider the largest number of possible models. In the cosmological context, these should include the $\Lambda$CDM and all its possible extensions, plus scenarios with any possible modified gravity paradigm and their extensions: this is clearly computationally impossible. From an Occam’s razor perspective, however, the models with an unnecessarily large number of parameters will be generally penalized by the Bayesian evidence calculation 333This is true as long as the additional parameters are constrained by the data, as no penalty applies to unconstrained parameters. The selection of models considered in the analysis should take into account equal priors only for reasonable extensions of the simplest model., so that their final weight in Eq. (7) will be negligible, while most of the contribution will be given by the most economical models that better fit the data. While our method allows to marginalize over the freedom related to different models or additional parameters, since it is based on the comparison of Bayesian evidences obtained in the different models, it still has a residual dependence on the shape and the width of the adopted priors.

Cosmological data analyses— The data we shall exploit to derive model-marginalized constraints from cosmological observations include measurements of the CMB angular power spectrum and of the BAO signature in the matter power spectrum. Awaiting for the final release from the Planck collaboration, we use here their 2015 data release Adam et al. (2016); Ade et al. (2016a). We consider two possibilities: a) both temperature and low-$\ell$ polarization (CMB), or b) temperature and polarization at all multipoles (CMB+pol). In both cases we also include the Planck CMB lensing determination (lens) Ade et al. (2016b). BAO geometrical information from the SDSS BOSS DR11 Anderson et al. (2014), the 6DF Beutler et al. (2011) and the SDSS DR7 MGS Ross et al. (2015) surveys complements the data sets used in our numerical analyses. We are aware that this combination may not provide the strongest cosmological constraints. However, it is not our main goal here to outperform the current cosmological constraints, but to exemplify the novel model-marginalized approach here proposed. After the Planck final public release, our method will be applied to an extended set of cases with respect to those considered here.

In our numerical calculations we use the Boltzmann solver CAMB Lewis et al. (2000) together with CosmoMC Lewis and Bridle (2002), with PolyChord Handley et al. (2015b, a) (version 1.9) as the algorithm devoted to extract the Bayesian evidences.

In our demonstrative analysis, we restrict our set of models to the simplest $\Lambda$CDM model with freely varying neutrino masses and some of its one-parameter extensions. In particular, we consider the $\Lambda$CDM+$\Sigma m_{\nu}$, $\Lambda$CDM+$\Sigma m_{\nu}$+$A_{\rm{lens}}$, $\Lambda$CDM+$\Sigma m_{\nu}$+$N_{\rm{eff}}$ and $\Lambda$CDM+$\Sigma m_{\nu}$+$w$ models, as discussed more in detail in the next paragraphs. In the numerical calculations, all the parameters that are shared among the different models are sampled adopting the same linear priors as in the default PolyChord settings, except for the sum of the neutrino masses which is varied in the range $[0.06,\,5]$ eV. For the additional parameters we adopt linear priors in the following ranges: $A_{\rm{lens}}$ varies in $[0,\,5]$, $N_{\rm{eff}}$ in $[1,\,5]$ and $w$ in $[-3,\,0]$.

Results: the neutrino mass as a case study— Table Cosmology-marginalized approaches in Bayesian model comparison: the neutrino mass as a case study summarizes the results from our novel method applied to a particular physics case that is usually constrained by cosmological observations: the sum of the neutrino masses $\Sigma m_{\nu}$. As aforementioned, a robust model-marginalized limit on $\Sigma m_{\nu}$ is absolutely required, as it is crucial for a number of issues. In particular, it is a very important input when deciding the experimental strategy for neutrino character (Dirac versus Majorana) searches. We show such model-marginalized limit in the second-to-last row of Tab. Cosmology-marginalized approaches in Bayesian model comparison: the neutrino mass as a case study, for the two data combinations considered here.