A non-parametric consistency test of the $\Lambda$CDM model with Planck CMB data

TL;DR

This paper uses Gaussian process regression to perform a non-parametric, model-independent consistency test of the Planck CMB data against the base CDM model, finding no significant inconsistencies.

Contribution

It introduces a novel application of Gaussian process reconstruction as a model-independent consistency test for cosmological data against CDM predictions.

Findings

No significant inconsistencies found in Planck data with CDM model

Supports CDM as the standard cosmological model

Demonstrates the effectiveness of GP-based residual analysis

Abstract

Non-parametric reconstruction methods, such as Gaussian process (GP) regression, provide a model-independent way of estimating an underlying function and its uncertainty from noisy data. We demonstrate how GP-reconstruction can be used as a consistency test between a given data set and a specific model by looking for structures in the residuals of the data with respect to the model's best-fit. Applying this formalism to the Planck temperature and polarisation power spectrum measurements, we test their global consistency with the predictions of the base $\Lambda$CDM model. Our results do not show any serious inconsistencies, lending further support to the interpretation of the base $\Lambda$CDM model as cosmology's gold standard.

Figures

Tables

Table 1. Table 1 : Prior ranges for the hyperparameters.

	TT	TE	EE
${\log }_{10}{\sigma }_f$	$[-2,2]$	$[-2,1.2]$	$[-2,1.5]$
${\log }_{10}{l}_f$	$[0,4]$	$[0,4]$	$[0,4]$

Equations

$$\left[\begin{array}{c}\text{mathbfcal}{C}_{\mathbf{\ell }}\\ \text{f}\end{array}\right]=\mathcal{N}\left(\left[\begin{array}{c}\text{m}(\text{L})\\ \text{m}({\text{L}}^{\prime })\end{array}\right],\left[\begin{array}{cc}{K}_y(\text{L},\text{L})& K(\text{L},{\text{L}}^{\prime })\\ K({\text{L}}^{\prime },\text{L})& K({\text{L}}^{\prime },{\text{L}}^{\prime })\end{array}\right]\right),$$

$$\bar{\text{f}}=\text{m}({\text{L}}^{\prime })+K({\text{L}}^{\prime },\text{L})K_y^{-1}(\text{L},\text{L})(\text{y}-\text{m}(\text{L}))$$

$$\text{cov}(\text{f})=K({\text{L}}^{\prime },{\text{L}}^{\prime })-K({\text{L}}^{\prime },\text{L})K_y^{-1}(\text{L},\text{L})K(\text{L},{\text{L}}^{\prime }).$$

$$k(\ell ,{\ell }^{\prime })={\sigma }_f^2\text{exp}\left(-\frac{(\ell -{\ell }^{\prime }{)}^2}{2l_f^2}\right),$$

$$\mathcal{L}({\mathcal{C}}_{\ell }|\text{L},l_f,{\sigma }_f)\equiv \int \mathrm{d}\mathbf{f}\mathcal{L}(\text{mathbfcal}{C}_{\mathbf{\ell }}|\text{L},\mathbf{f},l_f,{\sigma }_f)p(\mathbf{f}|\text{L}),$$

$$\ln \mathcal{L}=-\frac{1}{2}(\text{y}-\text{m}(\text{L}){)}^T{K}_y^{-1}(\text{L},\text{L})(\text{y}-\text{m}(\text{L}))-\frac{1}{2}\ln \mid K_y(\text{L},\text{L})\mid -\frac{N}{2}\ln 2\pi .$$

$${\mathcal{E}}_{\mathrm{GP}}=\int \mathrm{d}{\sigma }_f\mathrm{d}{l}_f\tilde{p}({\sigma }_f,l_f)\mathcal{L}({\sigma }_f,l_f),$$

$$B\equiv \frac{{\hat{\mathcal{E}}}_{\mathrm{GP}}}{{\hat{\mathcal{E}}}_{bg}}$$

Full text

1 Introduction

The Planck mission’s [1] measurement of the cosmic microwave background’s anisotropies is currently probably the most powerful source of cosmological information. One of the main results of Planck is a confirmation of the base $\Lambda$CDM model as the simplest phenomenologically viable cosmological model. This conclusion is based on (i) extensive testing for known systematic effects at all stages of Planck data analysis from time-ordered data to parameter estimation [2], and (ii) a wide-ranging exploration of the space of physically motivated models [3, 4, 5].

Nonetheless, while the overall goodness-of-fit of the Planck data to base $\Lambda$CDM is reasonable, one cannot entirely rule out the possibility of (i) remaining unknown systematics, or (ii) the existence of a model that would be preferred over base $\Lambda$CDM by some measure such as the Bayesian model evidence.

Here we argue that parametric methods are not amenable to finding such effects and suggest a test of the consistency of Planck data with base $\Lambda$CDM based on a non-parametric method which can be used to detect possible hidden patterns in the residuals of Planck’s angular power spectra. A statistically significant discovery of such a pattern would indicate problems with the model or our understanding of the data.

This test should be considered another item in the long list of challenges to base $\Lambda$CDM. Naturally, owing to any method’s inherent limitations, passing such a test cannot be an ultimate proof of base $\Lambda$CDM being correct, but would still bolster confidence in our present understanding of the physics governing the history of the Universe being adequate.

In particular, we employ Gaussian Process regression (see, e.g., [6, 7, 8, 9, 10, 11]) as a non-parametric regression tool to quantify the consistency of Planck 2015 temperature and polarisation angular power spectrum measurements with base $\Lambda$CDM.

Previous cosmology model-independent and non-parametric analyses of the Planck 2013 data in Refs. [12, 13] identified an inconsistency between the concordance model and the data at a $2$-$3\sigma$ level using different statistical approaches. However, further analysis showed that a systematic effect in the 217 GHz spectrum was mainly responsible for the inconsistency [14].

This paper is organized as follows: in Section 2 we explain the characteristics of the angular power spectrum data sets used in this work, elaborate the methodology of our Gaussian process-based consistency test and demonstrate its ability to find certain classes of patterns. In Section 3 we show and discuss the results of applying the test to Planck angular power spectrum data, and finally in Section 4 we present our conclusions.

2 Data and Methodology

In this section we introduce the CMB angular power spectrum data sets used in this analysis. We also briefly review the theory behind Gaussian process regression (see [6]) and explain the way we utilize it to test the consistency of the $\Lambda$CDM model with CMB angular power spectrum data.

2.1 Data

In this analysis we consider measurements of the CMB’s TT, TE and EE angular power spectra released by the Planck collaboration in 2015. We extract the data vectors and corresponding covariance matrices from the plik_lite likelihood available on the Planck Legacy Archive.111http://pla.esac.esa.int These data combine the information from Planck’s frequency channels and are already marginalised over foreground contributions.

The TT power spectrum data cover the multipole range $30\leqslant\ell\leqslant 2508$ and are given in $215$ bins with $\ell$-dependent widths, increasing from $\Delta\ell=5$ at $\ell_{\mathrm{min}}$ to $\Delta\ell=33$ at $\ell_{\mathrm{max}}$. The EE and TE spectra are provided in $199$ bins of widths $5\leq\Delta\ell\leq 17$ over a range of $30\leqslant\ell\leqslant 1996$. We plot the residuals of the data with respect to the $\Lambda$CDM model’s best-fit spectra in Figure 1. Even if the data were unbinned, in this multipole range the data follow to an excellent approximation a Gaussian likelihood distribution [2]. This motivates our use of Gaussian process regression, as described in the following subsection.

2.2 Gaussian process regression

A Gaussian Process (GP) $f(x)$ is a stochastic process which describes the properties of functions in function space using Gaussian distributions [6]. Given a set of $N$ data points $\{x_{i},y_{i}\},1\leq i\leq N$ and a covariance matrix of the data, $\Sigma$, this formalism can be used as a nonparametric regression tool to estimate the underlying function via the mean $\bar{f}(x)$ and its uncertainty via the covariance $\mathrm{cov}(f)$.

In the following, we will identify $\mathbf{y}$ with the vector $\mathbf{\mathcal{C}_{\ell}}$ formed by the binned $\mathcal{C}_{\ell}$ of the Planck plik_lite likelihood, $\mathbf{x}$ with the vector $\mathbf{L}$ of the weighted averages $\bar{\ell}_{i}$ of the multipoles included in bin $i$ and $\Sigma$ with the corresponding covariance matrix.

In addition to the data vectors and the covariance matrix, we need to specify the mean of the Gaussian process m, which can in general be a function of $\ell$, and a kernel or covariance function $\text{cov}(f(\ell),f(\ell^{\prime}))\equiv k(\ell,\ell^{\prime})$.

The Gaussian process evaluated at multipoles $\mathbf{L}^{\prime}$ follows the joint probability distribution given by [6, 10]

[TABLE]

where $\mathcal{N}$ is a multivariate normal distribution, $K_{y}(\textbf{L},\textbf{L})\equiv K(\textbf{L},\textbf{L})+\Sigma$ with the entries of $K(\textbf{L},\textbf{L})$ given by $[K(\textbf{L},\textbf{L})]_{ij}=k(\ell_{i},\ell_{j})$.

The mean of f given the data can then be expressed as

[TABLE]

with covariance

[TABLE]

We take a Gaussian kernel function,

[TABLE]

characterised by two hyperparameters: firstly, the correlation length $l_{f}$, which determines the stiffness of the Gaussian process and sets the typical scale on which the resulting $\bar{\mathbf{f}}$ can vary significantly and secondly, the prior width $\sigma_{f}$, which limits the allowed deviation of $f$ from the mean function.222Note that in Ref. [6] the $l_{f}$ and $\sigma_{f}$ are called length-scale and signal variance, respectively. The choice of hyperparameters strongly influences the outcome of the regression; for example, a correlation length much shorter than the typical distance of two data points might lead to an overfitted function with too much structure, whereas a correlation length much larger than the total range covered by the data could result in an $\bar{f}$ that fits the data very poorly.

However, instead of simply arbitrarily fixing the hyperparameters, one can make use of the marginal likelihood to let the data decide on their most suitable values. The marginal likelihood is defined as the conditional probability of the data given a set of hyperparameters, marginalised over all realisations of the GP,

[TABLE]

with a Gaussian prior $p$, and can be written as

[TABLE]

Here, we perform the GP regression with the combination of $\sigma_{f}$ and $l_{f}$ that maximises $\mathcal{L}(\sigma_{f},l_{f})$333Alternatively, in a more strictly Bayesian spirit, one could also consider marginalising $f$ over the space of hyperparameters, weighted with $\mathcal{L}$, i.e., $f_{\mathrm{marg}}\propto\int\mathrm{d}\sigma_{f}\mathrm{d}l_{f}f(\sigma_{f},l_{f})\mathcal{L}(\sigma_{f},l_{f}).$.

2.3 Using GP regression as a consistency test

In fact, the marginal likelihood $\mathcal{L}(\sigma_{f},l_{f})$ holds additional useful information that allows us to construct a consistency test: the location of the maximum can tell us whether the data prefer extra structure that is not present in the assumed mean function.

Let us begin by discussing the expected form of $\mathcal{L}(\sigma_{f},l_{f})$ in the absence of extra structure: consider first the limiting case where one chooses the mean function to go exactly through all data points, i.e., $\mathbf{m(L)}=\mathbfcal{C}_{\mathbf{\ell}}$. In this case, any deviation from $\mathbf{m(L)}$ will decrease the likelihood, and thus for fixed $l_{f}$, the marginal likelihood will punish the addition of "extra covariance" and prefer a kernel function (Equation (2.4)) equal to zero. At the same time, for fixed $\sigma_{f}$, it will prefer a stiffer GP over a more pliable one. In other words, the global maximum of $\mathcal{L}$ will be at $\{\sigma_{f}\}\rightarrow\{0\}$. For data which do scatter around the background but are still consistent with the background function, the global maximum of $\mathcal{L}(\sigma_{f},l_{f})$ may be a shallow bump at finite values of $\sigma_{f}$ and $l_{f}$ instead.

Conversely however, if $\mathcal{L}$ has a sufficiently pronounced peak, this may hint at a discrepancy between the data and the mean function, with the best-fit $\{\sigma_{f},l_{f}\}$ indicating the typical amplitude and scale of the deviation, respectively. Thus, if we set the mean function to be the best-fit of a particular model to the data, an inspection of $\mathcal{L}(\sigma_{f},l_{f})$ can reveal the presence of hidden patterns in the data and hence serve as a consistency test between model and data [15].

Since our goal is to probe the consistency of Planck data with the $\Lambda$CDM model, it is natural to set the mean function to the respective base $\Lambda$CDM best-fit angular power spectra. In particular, we use the primary $\mathcal{C}_{\ell}$ that maximise the Planck TTTEEE+lowP likelihood, downloaded from the Planck Legacy Archive.444Since the 2015 Planck data release, results from a new low-$\ell$ polarisation likelihood (SimLow) have been presented by the Planck collaboration in Ref. [16], but the neither the likelihood nor the best-fit parameters or spectra are in the public domain at the time of writing. Despite the significant shift in the inferred optical depth to reionisation $\tau$, there appears to be no degradation to the quality of the fit to the high-$\ell$ Planck data [16]. We thus would not expect any qualitative changes to our results if we were to use the Planck TTTEEE+SimLow best-fit spectrum instead. Note that this is equivalent to performing a GP regression with mean zero on the residuals of the data with respect to the base $\Lambda$CDM best fit.

In Figure 2 we illustrate the ability of this method to detect a hidden signal. In this particular example, we modulate the Planck TT data with a cosine function (see the lower panel of the top plot in Figure 1). Even for a relatively modest modulation amplitude of $8\,\mu\mathrm{K}^{2}$, the effective $\Delta\chi^{2}\equiv-2\Delta\ln\mathcal{L}\equiv-2\left(\ln\mathcal{L}(\sigma_{f},l_{f})-\ln\mathcal{L}_{\Lambda\mathrm{CDM}}\right)$ is smaller than $-25$, providing a clear signal. We can also observe the typical limiting behaviour of $\Delta\chi^{2}\rightarrow 0$ as $\sigma_{f}\rightarrow 0$ and $l_{f}\rightarrow\infty$.

2.4 Using the Bayesian evidence for model comparison

When it comes to deciding whether or not a GP-improved model is statistically preferred over the fit of the mean function to the data, one can take the idea behind Equation (2.5) one step further, and marginalise over the hyperparameter directions as well, leading to the Bayesian evidence,

[TABLE]

where one needs to specify a prior probability density $\tilde{p}$. This is to be compared to the evidence for the background model $\mathcal{E}_{\mathrm{bg}}$, which, in this case, is equivalent to the likelihood of the data for the base $\Lambda$CDM best-fit). Assuming equal model-priors, the Bayes factor

[TABLE]

quantifies the relative probabilities of the two models [17], with $\ln B>0$ indicating a preference for the GP-model.555Note that in our case we can exploit the fact that the models are nested, since the GP model reduces to the background model in the limit $\sigma_{f}\rightarrow 0$. It is thus not necessary to perform the integration over $\Lambda$CDM parameter space $\mathbf{\Theta}$ to calculate the full evidences $\hat{\mathcal{E}}_{\mathrm{GP}}\equiv\int\mathrm{d}\mathbf{\Theta}p(\Theta)\mathcal{E}_{\mathrm{GP}}(\Theta)$ and $\hat{\mathcal{E}}_{bg}\equiv\int\mathrm{d}\mathbf{\Theta}p(\Theta)\mathcal{L}(\mathbf{\Theta})$; one can instead evaluate $B$ through the Savage-Dickey density ratio [17], i.e., $B=\left.\frac{\mathcal{E}_{\mathrm{GP}}}{\mathcal{L}(\sigma_{f},l_{f})}\right|_{\sigma_{f}=0}$. We do emphasize that the Bayes factor can have a considerable dependence on the choice of $\tilde{p}$, and should be taken with a grain of salt. In the following, we will assume top-hat priors on the logarithms of the hyperparameters with the (data set dependent) prior ranges listed in Table 1.

In the fake-signal example of Section 2.3, these priors lead to Bayes factors of $\ln B=0.55$ for a modulation amplitude of $A=5\,\mu\mathrm{K}^{2}$ (a very mild preference for the GP-model), and $\ln B=8.3$ for the case of $A=8\,\mu\mathrm{K}^{2}$ (decisive evidence for the GP-model).

3 Results and Discussion

In Figure 3, we show the profiles of the logarithm of the marginal likelihood in the $\sigma_{f}$- and $l_{f}$-directions for the three sets of angular power spectrum residuals. For the best-fit hyperparameters, we plot the resulting mean and 1- and 2-$\sigma$ uncertainty bands of the Gaussian process, along with the residuals of the data, in Figure 1. As can be seen from Equation (2.3), if the covariance of the data $\Sigma$ dominates over the inherent covariance of the GP $K$, as is the case for the TT and EE data, the variance of the GP is given by $\sigma_{f}^{2}$. In other words, for these particular choices of hyperparameters, the data are not able to meaningfully constrain the GP beyond its original prior width, which might be a consequence of using a constant $\sigma_{f}$ (as opposed to a function $\sigma_{f}(\ell)$) for data with a strongly $\ell$-dependent variance. This is not the case for TE, where the size of the error bars is more uniform over the range of $\ell$ spanned by the data, and hence the variance of the GP is smaller than $\sigma_{f}^{2}$ here.

We note that the minima of the projection of $-2\ln\mathcal{L}$ in the $l_{f}$-direction do not lie at $\sigma_{l}=0$, but at $(\sigma_{f},l_{f})=(2.4,6.0)$ for TT, $(1.4,40)$ for TE and $(0.019,\infty)$ for EE (see also Figure 4 for the joint 2-dimensional likelihood). Taken at face value, this might suggest some preference for additional covariance in the residuals. However, when dealing with real data which scatter around the true model, this should not come as a surprise. We find that effective $\Delta\chi^{2}$ of the optimal Gaussian process regression with respect to the base $\Lambda$CDM best fit is $\Delta\chi^{2}_{\mathrm{eff}}=-0.4$ for the TT data, $\Delta\chi^{2}_{\mathrm{eff}}=-3.2$ for TE and $\Delta\chi^{2}_{\mathrm{eff}}=-0.2$ for EE, respectively. Such mild improvements in the marginal likelihood cannot be interpreted as evidence for extra structure and are likely to be due to the stochastic spread of the data.

And indeed, for our choice of prior ranges, the Bayesian evidence does not show any preference for the GP-model either, with $\ln B^{\mathrm{TT}}=-0.5$ for the TT data, $\ln B^{\mathrm{TE}}=-0.3$ for the TE data and $\ln B^{\mathrm{EE}}=-1.6$ for the EE data.

In the case of the TT data, this conclusion is further supported by three observations that are related to the smallness of $\Delta\chi_{\mathrm{eff}}^{2}$: firstly, the optimal values of $l_{f}$ are $\mathcal{O}(10)$, comparable to the average bin size. Secondly, the amplitude of the modulation is much smaller than the mean standard deviation of the residuals (cf. Figure 1). Thirdly, the maximum local deviation of the GP from the mean does not exceed one standard deviation, with the global significance being even smaller due to the look-elsewhere effect. The latter two points are also true for the EE case, but here, instead of adding structure, the GP prefers a nearly constant offset, albeit with an amplitude that is completely negligible when compared to both the amplitude of the EE power spectrum and the standard deviation of the residuals. We should note that a similar constant amplitude offset was reported recently by [18]. While the mean amplitude of the offset in this analysis seems to be smaller than the one found in [18], results of the two analyses are consistent considering the uncertainties. This effect could potentially be caused by a permille-sized mismatch in the calibration of the polarisation spectra with respect to the temperature spectra, or might also be related to the absence of the low-$\ell$ Planck data in our approach (note that the base $\Lambda$CDM best-fit spectra we assume are the best-fits with respect to the full Planck likelihood).

Things are slightly more interesting for the TE data: the preferred correlation length is clearly larger than the average bin size and the optimal prior width reaches about 22% of the average standard deviation of the residuals. Nonetheless, neither the local deviation of the GP from the base $\Lambda$CDM best-fit nor the absolute $\Delta\chi^{2}_{\mathrm{eff}}$ indicate a serious discrepancy. A likely reason for this could be the fact that the base $\Lambda$CDM model’s best-fit to TE only does not completely coincide with the full TTTEEE data fit, which is still dominated by the TT data (explaining why the TT data GP does not display the same behaviour). But at the same time, the TE data by themselves are starting to become competitive with the TT data for parameter estimation purposes [3], so unlike the EE data, they are sensitive enough for this small difference in best-fits to actually be picked up by the GP regression test.

4 Conclusion

We have tested the consistency of Planck TT, TE and EE angular power spectrum data with the base $\Lambda$CDM model using a non-parametric test based on Gaussian process regression. Being non-parametric and applied directly to the observable (rather than a fundamental quantity, e.g., the primordial power spectrum), this method is sensitive to general inconsistencies between data and model, no matter whether they are caused by assuming a wrong physical model or unknown systematics in the data.

It should be pointed out that since we used binned data, the scope of our analysis is limited to structures with a characteristic variation length of the order of the respective bin size. It would be interesting to repeat the analysis on unbinned data in order to increase sensitivity to finer structures with $\Delta\ell\sim 1$, but unfortunately, Planck’s frequency-combined, foreground-marginalised power spectra are not publically available in unbinned form at present.

Nevertheless, our results do not show any serious inconsistencies; besides a statistically non-significant deviation in the TE spectrum, we find excellent agreement between Planck data and the best-fit base $\Lambda$CDM spectra: especially the TT residuals display a perhaps remarkable lack of hidden patterns. We therefore conclude that the $\Lambda$CDM model passes yet another challenge thrown at it, and confirm its consistency with Planck data.

5 Acknowledgments

J.H. gratefully acknowledges support of an Australian Research Council Future Fellowship (FT140100481) and would like to thank Aarhus University and the Korea Astronomy and Space Science Institute for hospitality during the preparation of this work. A.S. would like to acknowledge the support of the National Research Foundation of Korea (NRF-2016R1C1B2016478) and Sydney Institute for Astronomy, University of Sydney for hospitality during the preparation of this work.