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

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.

We propose here a novel method which singles out the 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.
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. [1]. 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 [2][3][4][5][6][7].
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 Σm ν (see Refs. [8][9][10] 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 Σm ν < 0.24 eV [11] 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 Σm ν < 0.12 eV at 95% CL [11], or even down to Σm ν < 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?
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 Σm ν when studying a number of possible cosmological scenarios, i.e. the minimal ΛCDM universe with massive neutrinos and its extensions. Adopting Planck 2015 data [59], the tightest bound we obtain within a ΛCDM universe is Σm ν < 0.23 eV at 95% CL, which relaxes to Σm ν < 0.35 eV when the uncertainty on the cosmological model is taken into account using our model-marginalization method. arXiv:1812.05449v2 [astro-ph.CO] 10 Jan 2019 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 ΛCDM+Σm ν model, which arises when it is compared with each of its extensions individually, will not correspond to a global 1σ level strength when considering the entire ensemble of extended scenarios.
Bayesian statistics-The Bayes theorem, which represents the foundation of Bayesian statistics, reads: where π(θ|M i ) and p(θ|d, M i ) are the prior and posterior probabilities for the parameters θ within a model M i , L(θ) is the likelihood as a function of the parameters θ, given the data d and the model M i , and . The Bayes theorem can also be written in a slightly different form to obtain the model posterior probability [60]: where π i ≡ π(M i ) refers to the model prior probability. In the Bayesian model comparison framework, the socalled Bayes factor provides a measure of whether the data have increased or decreased the odds of model M i relative to a second model M j : 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: 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 Adopting the commonly exploited Jeffreys' scale [61], the strength of the posterior odds can be ranked as inconclu- . 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 1σ (inconclusive) to 3σ (strong) level when considering a Gaussian variable.
Using Eq. (2), and selecting one among the available models, labelled M 0 without loss of generality, one can write, provided all priors are identical for all models: where we have used the definition of the Bayes factor. Notice that the posterior probability of the selected model M 0 depends on the Bayes factors with respect to all the possible models. For each data combination we will choose M 0 to be the preferred model. In practice, this is the one that has more influence on the modelmarginalized posterior. Since the model M 0 is the preferred one, we will always have 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 | 5 for i = 0, thus indicating apparently strong results according to the usually adopted Jeffreys' scale, the posterior probability of M 0 is p 0 0.955, which would indicate a mild 2σ significance for a Gaussian measure. In the same way, having N = 7 and | ln B i0 | 2.5 for i = 0, which usually indicates a weak preference, would give p 0 0.67, which would correspond to less than 1σ preference for M 0 . The tools of model comparison also allow us to compute a model-marginalized posterior distribution for the parameter θ, taking into account the posterior probability of each model M i resulting from the data d [1]: where the posterior probabilities of θ within each model M i are weighted according to the model posterior probabilities p i . These can be written using Eq. (2) to obtain the fundamental formula This is the expression that we will use to obtain modelmarginalized 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 Λ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 [62], 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 [63,64]. We consider two possibilities: a) both temperature and low-polarization (CMB), or b) temperature and polarization at all multipoles (CMB+pol). In both cases we also include the Planck CMB lensing determination (lens) [65]. BAO geometrical information from the SDSS BOSS DR11 [66], the 6DF [67] and the SDSS DR7 MGS [68] 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 modelmarginalized 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 [69] together with CosmoMC [70], with PolyChord [60, 71] (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 ΛCDM model with freely varying neutrino masses and some of its oneparameter extensions. In particular, we consider the ΛCDM+Σm ν , ΛCDM+Σm ν +A lens , ΛCDM+Σm ν +N eff and ΛCDM+Σm ν +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 lens varies in [0,5], N eff in [1,5] and w in [−3, 0].
Results: the neutrino mass as a case study- Table I 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 Σm ν . As aforementioned, a robust model-marginalized limit on Σm ν 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. I, for the two data combinations considered here.
In order to compute the model-marginalized result we consider, together with the simplest ΛCDM model with freely varying neutrino masses, some cosmological scenarios which are usually explored in the literature, see e.g. Planck 2015 data analyses [64]. These models contain extra parameters which are either partially or significantly degenerate with the neutrino mass. For instance, Σm ν has a correlation, among others, with the phenomenological parameter A lens , which rescales the lensing amplitude in the CMB spectra. Since current CMB constraints on the neutrino mass are mostly due to the reduction in the lensing potential induced by a larger neutrino mass, there is a degeneracy between Σm ν and A lens : a value A lens < 1 (> 1) would allow for a lower (higher) value of Σm ν . Notice from the results depicted in Tab. I that the neutrino mass bounds, in absence of high-multipole polarization, are worsened when the A lens parameter is allowed to vary. When going from a ΛCDM to a ΛCDM+A lens scenario the 95% CL limit changes from Σm ν < 0.28 eV to Σm ν < 0.38 eV [72]. Another parameter potentially degenerate with Σm ν is the number of relativistic degrees of freedom N eff , albeit the latest Planck analyses have shown that data are able to disentangle between the different physical effects induced by Σm ν and N eff [11] on temperature and polarization anisotropies. While the 95% CL limit without high-polarization is Σm ν < 0.37 eV, information from high multipoles brings the neutrino mass constraint extremely close to the bound obtained within the ΛCDM model [73]. Finally, a freely-varying constant dark energy equation of state w can also affect the bounds on Σm ν . If w is allowed to vary, the matter energy density can take very high values, compensating for the suppression induced in large scale structure due to an increased value of Σm ν , and therefore these two parameters will be correlated in a significant way. As a result, the limit is relaxed to Σm ν < 0.42 eV at 95% CL, both without and with high-multipole polarization data, see Tab. I.
From all the limits above and using the Bayes factors also listed in Tab. I, by means of Eq. (7) it is possible to obtain the marginalized limits on Σm ν shown in the second-to-last row of Tab. I. Notice that the 95% CL upper limits obtained within the most economical ΛCDM picture are significantly relaxed (they are increased up to 50%) when considering extended scenarios. For a visual  [63][64][65] plus BAO measurements [66][67][68] (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 ΛCDM model, while the last two rows illustrate, respectively, the model-marginalized 95% CL limit obtained via Eq. (7) and the posterior probability of the example model M0 (the preferred one, that is always the ΛCDM+Σmν scenario), see Eq. (5).
comparison of the one-dimensional posterior probabilities of Σm ν in the various models considered here and of the model-marginalized one, we provide Figures 1 and 2, where we also show the sampled prior distribution [74] Notice that the method following Eq. (7) allows a proper weighing of the information from each model, building a robust estimate for the neutrino mass that can be used as an input in neutrino particle physics. The possible applications of the method, however, are significantly wider than what explored here. The last row of Tab. I shows the posterior probabilities p 0 for the example model M 0 , computed from Eq. (5). M 0 is chosen to be the preferred one by each of the two data combinations, and it turns out to be the minimal ΛCDM scenario with free neutrino masses in both cases. The posterior probability, which depends on the Bayesian evidences of various models, is shown in the second and fourth columns for the two possible data combinations. Here one should clarify an important aspect of Bayesian model comparison. While the Bayes factors with respect to the extended models, if considered separately, indicate a weak-to-moderate [75] Bayesian preference for the ΛCDM model accordingly to the Jeffreys' scale [61], and therefore individually corresponding to a 1.1−2.7σ probability (in Gaussian terms) in favor of the ΛCDM framework, it is clear from Eq. (5) that such naive expectations are no longer true when more than one model is accesible. The values of the posterior probabilities for the example ΛCDM model never reach the 1σ level strength in terms of a Gaussian variable. Based on these results, therefore, it is possible to say that the ΛCDM model, despite being more likely than its extensions, is not strongly preferred by data. This is a crucial result of our analyses, with strong implications in many other early universe fundamental physics searches, as for example in the case of the inflationary landscape, where many models arise and are usually ranked by means of Bayesian comparison techniques, see e.g. Ref. [76]. One-dimensional posterior probabilities for Σmν for different cosmological models, arising from Planck 2015 CMB temperature, low-polarization and lensing data [63][64][65] plus BAO measurements [66][67][68]. We also depict the model-marginalized bound obtained using Eq. (7) and the prior sampled on Σmν , see text for details.

Discussion
robust machinery to compute model marginalized limits. We have proposed a method which allows to minimize the uncertainty related to multiple model choices on the determination of parameter constraints. We have applied our novel method to the neutrino mass case, exploiting current publicly available cosmological data. We show that the limits on the neutrino masses can significantly change when one realizes that present measurements are not able to unambiguously tell us the cosmological model that nature has chosen. The statistical Gaussian preference for the favoured model, indeed, always becomes inconclusive when there are a number of other possible models, even if equally disfavored by observations. An updated and extended analysis using the proposed method will come after the release of Planck 2018 likelihoods. As Fig. 1, but including also CMB high multipole polarization measurements from the Planck 2015 data release [63,64].