EVOLUTION OF THE COSMIC MICROWAVE BACKGROUND POWER SPECTRUM ACROSS WILKINSON MICROWAVE ANISOTROPY PROBE DATA RELEASES: A NONPARAMETRIC ANALYSIS

We are given CMB angular power spectrum data of the form

$$\begin{equation} Y_l = C_l + \epsilon _l \end{equation} \tag{ 1 }$$

consisting of N data points observed over multipole index range l_min ⩽ l ⩽ l_max. Here, C_l stands for the value of the true but unknown power spectrum at multipole index l. The noise variables _l are assumed to have a mean-0 normal distribution with known covariance matrix Σ. In practice, any reasonable estimate/approximation $\widehat{\Sigma }$ of this covariance matrix, such as an inverse Fisher matrix for a pilot parametric fit, can be used in place of Σ.

This nonparametric regression method is based on expanding the unknown regression function f, assumed to belong to an appropriate L₂ function space, in a complete orthonormal basis {ϕ_j(x)}, as

$$f(x) = \sum _{j=0}^\infty \beta _j \phi _j(x).$$

A basis that has proven useful in the CMB angular power spectrum context is the cosine basis defined over 0 ⩽ x ⩽ 1:

$$\begin{equation} \phi _j(x) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (j=0) \\ \sqrt{2} \cos (j\pi x) \; (j = 1, 2, \ldots)\end{array}\right.\!\!\!\! . \end{equation} \tag{ 2 }$$

Assuming that f is sufficiently smooth, we take

$$f(x) \approx \sum _{j=0}^{N-1} \beta _j \phi _j(x),$$

and estimate it as

$$\begin{equation} \widehat{f}_N(x) = \sum _{j=0}^{N-1} \widehat{\beta }_j \phi _j(x). \end{equation} \tag{ 3 }$$

While the method is asymptotically (i.e., as the data size N → ∞) basis-independent, choice of the basis may matter in any finite-N application; see Beran (2000a, 2000b) for a detailed discussion. This basis satisfies a discrete orthogonality property when the data Y_i are available over an equispaced grid {x_i = (2i + 1)/2N, 0 ⩽ i ⩽ N − 1} consisting of zeros of ϕ_N(x). In the CMB context, any contiguous range of N integer multipole indices l_min ⩽ l ⩽ l_max can be formally mapped onto this equispaced grid, hence we will not make any categorical distinction between data index i and the corresponding multipole index l.

The true angular power spectrum C_l ≡ f(x_l) is estimated as $\widehat{C}_l \equiv \widehat{f}_N(x_l)$ via coefficient estimates $\widehat{\beta }_j$, which are estimated as

$$\begin{equation} \widehat{\beta }_j = \lambda _j Z_j, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} Z_j = {1 \over N} \sum _{i=0}^{N-1} Y_i \phi _j(x_i) = {(U^T Y)_j \over \sqrt{N}}. \end{equation} \tag{ 5 }$$

Here, U is the orthonormal matrix with elements $U_{ij} = \phi _j(x_i) / \sqrt{N}$ and Y ≡ (Y₀, ..., Y_{N − 1})^T.

The task of obtaining coefficient estimates $\widehat{\beta } \equiv (\widehat{\beta }_0, \ldots, \widehat{\beta }_{N-1})^T$, and thereby the fit $\widehat{f}_N(x_i)$, is now relegated to determining the shrinkage parameters λ_j. Assuming smoothness for f (which implies a rapid decay of the true coefficients β_j with j), the shrinkage parameters λ_j are constrained to be monotonically decreasing

$$\begin{equation} 1 \ge \lambda _0 \ge \lambda _1 \ge \ldots \ge \lambda _{N-1} \ge 0 \;\; \mbox{(Monotone shrinkage)}. \end{equation} \tag{ 6 }$$

A special, discrete subset of the monotone shrinkage defined above is the nested subset selection (NSS) shrinkage, defined as

$$\begin{equation} \lambda _j = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 \quad\mbox{ for $0 \le j < J$} \\ \\[-6pt] 0 \quad\mbox{ for $J \le j < N$} \end{array}\right. \mbox{(NSS shrinkage)}. \end{equation} \tag{ 7 }$$

Either shrinkage type results in selective damping of high-frequency harmonics in the data Y_i, which results in smoothing of the fit $\widehat{f}_N$. A useful interpretation of shrinkage parameter λ_j is that it represents the effective degree of freedom for the jth coefficient estimate $\widehat{\beta }_j$. One can thus define the effective degrees of freedom (EDoF) for the entire fit $\widehat{f}_N$ as

$$\begin{equation} \mbox{EDoF}(\lambda) = \sum _{j=0}^{N-1} \lambda _j. \end{equation} \tag{ 8 }$$

This definition follows from the fact that for a linear smoother, EDoF of a fit is formally defined as $\mbox{tr}(\mathcal {H})$, where $\mathcal {H}$ is the matrix that connects the fitted values $\widehat{Y}$ to the data Y as $\widehat{Y} = \mathcal {H} Y$. For the present nonparametric regression method, $\mathcal {H}=U D U^T$, where U is the orthonormal basis matrix (Equation (5)), and D ≡ diag(λ₀, ..., λ_{N − 1}). This implies $\mbox{tr}(\mathcal {H}) = \mbox{tr}(D) = \sum _{i=0}^{N-1} \lambda _i$.

In the present formalism, the discrepancy between the (unknown) regression function f and its estimator $\widehat{f}_N$ is measured by the inverse-noise-weighted squared loss function L(λ), defined as

$$L(\lambda) = \int _0^1 { \left(f(x) - \widehat{f}_N(x) \over \sigma (x) \right)^2 } dx.$$

Here, σ²(x) is the (known) variance of the data Y at x, which accounts for the heteroskedasticity of the data Y. The loss L is considered a function of the vector of shrinkage parameters λ ≡ (λ₀, ..., λ_{N − 1})^T that determine the regression estimator $\widehat{f}_N$ via Equation (4). Risk R(λ), which is the expected value of L(λ), can be written as a sum of two non-negative terms; namely,

$$\begin{eqnarray*} R(\lambda) & = & \int _0^1 \left({ f(x) - \mathbb {E}(\widehat{f}_N(x)) \over \sigma (x) } \right)^2 dx\nonumber\\ && + \int _0^1 \mathbb {E}\left[ \left({ \widehat{f}_N(x) - \mathbb {E}(\widehat{f}_N(x) ) \over \sigma (x) } \right)^2 \right] dx. \end{eqnarray*}$$

These two terms represent, respectively, the integrated squared bias and the integrated variance of $\widehat{f}_N(x)$, both weighted by 1/σ²(x). Optimal smoothing is achieved, in principle, by minimizing R(λ) with respect to λ. Generally speaking, too little smoothing leads to a fit $\widehat{f}_N$ with low bias and high variance, and too much smoothing yields a fit with high bias and low variance. Minimal risk or optimal smoothing therefore can be thought of as a balance between the bias of $\widehat{f}_N$ and its variance.

The risk function R(λ), unfortunately, depends on the unknown regression function f, and therefore needs to be estimated. A particular estimator of this risk, which is of the SURE (Stein's unbiased risk estimator; Stein 1981) kind, takes the following form:

$$\begin{equation} \widehat{R}(\lambda) = Z^T \bar{D} W \bar{D} Z + \mbox{tr}(D W D B) - \mbox{tr}(\bar{D} W \bar{D} B), \end{equation} \tag{ 9 }$$

where Z ≡ (Z₀, ..., Z_{N − 1})^T, D ≡ diag(λ₀, ..., λ_{N − 1}), $\bar{D} = I_N - D$, B = U^TΣU/N is the covariance of Z, and I_N is the N × N identity matrix. The positive (semi)definite weight matrix W is defined as

$$\begin{equation} W_{jk} = \sum _l \Delta _{{jkl}} w_l, \end{equation} \tag{ 10 }$$

where w_l is the lth coefficient in the expansion (1/σ²(x)) ≈ ∑^{N − 1}_{j = 0}w_jϕ_j(x) and, for the cosine basis (Equation (2)),

$$\begin{eqnarray*} \Delta _{{jkl}} & = & \int _{0}^{1} \phi _j(x) \phi _k(x) \phi _l(x) \hspace{2.84526pt} dx\nonumber\\ & = & \left\lbrace \begin{array}{@{}ll}1, & \mbox{if }\#\lbrace j,k,l=0\rbrace =3, \\\\[-6pt] 0, & \mbox{if }\#\lbrace j,k,l=0\rbrace =2, \\\\[-6pt] \delta _{jk}\delta _{0l} + \delta _{jl}\delta _{0k} + \delta _{kl}\delta _{0j}, & \mbox{if }\#\lbrace j,k,l=0\rbrace =1, \\\\[-6pt] \frac{1}{\sqrt{2}} (\delta _{l,j+k}+\delta _{l,|j-k|}) & \mbox{if }\#\lbrace j,k,l=0\rbrace =0. \\ \end{array} \right. \end{eqnarray*}$$

We denote the optimal shrinkage obtained by minimizing risk $\widehat{R}(\lambda)$ by $\widehat{\lambda }$. The best NSS shrinkage $\widehat{\lambda }_{{\rm NSS}}$ is obtained simply by evaluating risk $\widehat{R}(\lambda)$ for each of the (N + 1) NSS shrinkage vectors and choosing the one with the least risk. Monotone shrinkage usually results in a lower risk than the NSS shrinkage because of the greater freedom available in the monotone set of shrinkages. We will discuss risk minimization subject to monotonicity constraints (Equation (6)) in Section 2.3.

Figure 1 illustrates the contrasting behavior of the nonparametric risk (red curve) and the WMAP 1-year likelihood function (green curve; Verde et al. 2003) for the WMAP 1-year data (Hinshaw et al. 2003b), as a function of the EDoF of all NSS fits. This figure is motivated by the fact that cosmologists, by and large, are better-acquainted with parametric likelihood-based methods. Each integer value on the horizontal axis represents one NSS fit, from the zero function at EDoF = 0 to the fit that simply interpolates through the data (EDoF = N). Optimal smoothing for NSS shrinkage occurs at EDoF = 12 where the nonparametric risk function attains its minimum over the NSS set of fits. The likelihood function, on the other hand, keeps on improving with the EDoF indefinitely.

Zoom In Zoom Out Reset image size

Conventional regression methods provide a confidence band around the fit that quantifies the uncertainty in the fit. In contrast, this nonparametric methodology quantifies the uncertainty surrounding the nonparametric fit in the form of an elegant construct, namely, a (1 − α) confidence set at a pre-specified confidence level 0 ⩽ (1 − α) ⩽ 1. The (1 − α) confidence set for the coefficient vector β is defined as

$$\begin{equation} \mathcal {D}_{N,\alpha } = \big\lbrace \beta : (\beta -\widehat{\beta })^T W (\beta -\widehat{\beta }) \le {r}^2_\alpha \big\rbrace, \end{equation} \tag{ 11 }$$

which is centered at the vector of estimated coefficients $\widehat{\beta }$, and the confidence radius r_α is given by

$$\begin{equation} {r}^2_\alpha = {\widehat{\tau } z_\alpha \over \sqrt{N}} + \widehat{R}(\widehat{\lambda }). \end{equation} \tag{ 12 }$$

Here, z_α is the upper α quantile of the standard normal distribution, and

$$\begin{equation} \widehat{\tau }^2 / N = 2 \mbox{tr}(ABAB) + Z^T Q Z - \mbox{tr}(QB), \end{equation} \tag{ 13 }$$

with Q = 4(ABA + WDBDW − 2ABDW) and A = DW + WD − W. In practice, the risk estimator (Equation (9)) and/or $\widehat{\tau }^2$ (Equation (13)) may turn out to be negative for particular data/covariance matrix realizations. In such cases, the squared confidence radius (Equation (12)) may be negative (or may not be 0 for α = 0). For minimization purposes, the risk estimator (Equation (9)) is adequate and appropriate (Beran & Dümbgen 1998). For confidence radius purposes, we suggest the following modifications to avoid the negativity problem:

$$\begin{eqnarray} \widehat{R}_+ & = & Z^T \bar{D} W \bar{D} Z + \max \lbrace 0, \mbox{tr}(D W D B) - \mbox{tr}(\bar{D} W \bar{D} B) \rbrace \nonumber \\ \widehat{\tau }_+^2 / N & = & 2 \mbox{tr}(ABAB) + \max \lbrace 0, Z^T Q Z - \mbox{tr}(QB) \rbrace \\ {r}^2_{\alpha +} & = & \max \left\lbrace 0, {\widehat{\tau }_+ z_\alpha \over \sqrt{N}} + \widehat{R}_+ \right\rbrace. \nonumber \end{eqnarray} \tag{ 14 }$$

At worst, this adjustment will make the confidence radius bigger, resulting in, e.g., wider confidence intervals, but more conservative inferences. Similar modifications have been suggested in Beran & Dümbgen (1998), Beran (2000a), and Genovese et al. (2004).

The corresponding confidence set on the true regression function f is given by

$$\begin{equation} \mathcal {B}_{N,\alpha } = \left\lbrace f(x) = \sum _{j=0}^{N-1} \beta _j \phi _j(x) : \beta \in \mathcal {D}_{N,\alpha } \right\rbrace. \end{equation} \tag{ 15 }$$

The quadratic form of the inverse noise-weighted loss function and the fact that the weight matrix W is positive (semi)definite implies that both confidence sets $\mathcal {D}_{N,\alpha }$ and $\mathcal {B}_{N,\alpha }$ are ellipsoidal in shape. For any functional T of the spectrum f, such as location or height of a peak or a dip, the (1 − α) confidence interval is defined as

$$\begin{equation} \mathcal {I}_{N,\alpha } = \bigg(\min _{f \in \mathcal {B}_{N,\alpha }} T(f), \max _{f \in \mathcal {B}_{N,\alpha }} T(f) \bigg). \end{equation} \tag{ 16 }$$

Moreover, prior information that the true regression function f belongs to a subset $\mathcal {P}_{N,\alpha }$ of the confidence set $\mathcal {D}_{N,\alpha }$ (e.g., f has k peaks over the range of x-values represented in the data) can be incorporated in the analysis by replacing $\mathcal {D}_{N,\alpha }$ with $\mathcal {P}_{N,\alpha } \cap \mathcal {D}_{N,\alpha }$. This methodology further provides the formal assurance that, asymptotically,

• 1.  
  $\mathcal {B}_{N,\alpha }$ ($\mathcal {D}_{N,\alpha }$) will contain the true spectrum f (true coefficient vector β) with probability ⩾(1 − α), and
• 2.  
  confidence intervals $\mathcal {I}_{N,\alpha }$ on any number of functionals T(f), computed from the confidence set $\mathcal {B}_{N,\alpha }$, will be simultaneously valid at the same confidence level (1 − α), and that these will trap their corresponding true but unknown values with probability ⩾(1 − α).

In this section, we show how the risk function $\widehat{R}(\lambda)$ (Equation (9)) can be minimized subject to the monotonicity constraints 1 ⩾ λ₀ ⩾ λ₁ ⩾ ... ⩾ λ_{N − 1} ⩾ 0. The risk function corresponding to the unweighted loss function (W = I_N) has a simple weighted-sum-of-squares form, and can be minimized exactly and efficiently using the pooled adjacent violators (PAV) algorithm (Robertson et al. 1988). While the risk function corresponding to the inverse-noise-weighted loss function (W ≠ I_N) is still quadratic in λ, it can no longer be expressed as a weighted sum-of-squares and the PAV algorithm cannot be used to minimize it.

It can be shown that, disregarding terms that do not depend on λ, the risk function $\widehat{R}(\lambda)$ (Equation (9)) can be written as

$$\widehat{R}(\lambda) = {1 \over 2} \lambda ^T H \lambda - \lambda ^T h,$$

where H_jk = 2z_jz_kW_jk, h = (H − V)(1, 1, ..., 1)^T, and V_jk = 2W_jkB_kj. H and V are both manifestly symmetric. Positive (semi)definiteness of W implies that H is a positive (semi)definite matrix, implying that $\widehat{R}(\lambda)$ is a convex function. The system of linear inequality constraints 1 ⩾ λ₀ ⩾ λ₁ ⩾ ... ⩾ λ_{N − 1} ⩾ 0 implies that the constrained region (Figure 2) has a convex trianguloidal shape determined by flat surfaces. The original risk minimization problem can therefore be formulated as the following equivalent convex quadratic minimization problem:

$$\begin{eqnarray} \mbox{Minimize}\quad & & \widehat{R}(\lambda) = {1 \over 2} \lambda ^T H \lambda - \lambda ^T h \nonumber \\ \mbox{subject to}\quad & & C \lambda \le (0,0,\ldots,0)^T \\ \mbox{and}\quad & & 0 \le \lambda _i \le 1 \;\; \mbox{for all $i$,} \nonumber \end{eqnarray} \tag{ 17 }$$

where C is the (N − 1) × N matrix

$$C = \left[ \begin{array}{rrrrrrr}-1 & 1 & 0 & \ldots & 0 & 0 & 0 \\ 0 & -1 & 1 & \ldots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & -1 & 1 & 0 \\ 0 & 0 & 0 & \ldots & 0 & -1 & 1 \\ \end{array} \right].$$

An additional linear inequality or equality constraint of the form

$$\begin{equation} \sum _{i=0}^{N-1} \lambda _i \le q \;\;\; \mbox{ or } \;\;\; \sum _{i=0}^{N-1} \lambda _i = q, \end{equation} \tag{ 18 }$$

where q constrains the EDoF of the fit (0 < q ⩽ N), can easily be accommodated in this formulation. This reformulation of the monotone risk minimization problem makes it possible to use standard minimization methods (Pshenichny & Danilin 1978; Powell 1985; Goldfarb & Idnani 1983) and computational tools (Schittkowski 2007; Vanderbei 1999) for convex quadratic minimization problems subject to linear constraints and simple bounds.

Zoom In Zoom Out Reset image size

To motivate the discussion in this section, consider the full-freedom monotone fit to the WMAP data sets obtained as above (green curves in Figures 3 and 4). By full-freedom monotone fit, we mean the fit that minimizes risk over the entire monotone-admissible region (Figure 2) without any restriction on the EDoF of the fit. This fit turns out to be quite wiggly especially at the high-l end because of the high noise variance here. Such wiggliness implies the presence of high-frequency components in the fit, which in turn implies a large number of EDoF in the fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Without cosmological pre-conditioning (i.e., from a completely agnostic and data-driven perspective) and when viewed in relation to the data, it is clear that this fit is not at all unreasonable, given the high noise levels at the high-l end. However, all cosmological models suggest far smoother shapes for the angular power spectrum. In the context of the present methodology, one candidate for a smoother fit is the NSS fit (i.e., one that has the minimal risk over the set of (N + 1) NSS fits). Indeed, this possibility has been exploited, e.g., in Genovese et al. (2004). However, the NSS fit (see the blue curve in Figures 3 and 4) may also turn out to be somewhat unsatisfactory with respect to cosmological expectations (and, some times, also with respect to trends reflected in the full-freedom monotone fit). This is primarily because of the limited freedom available in the NSS class. Note again that the NSS fit is not entirely unreasonable from an agnostic viewpoint.

The monotone set, on the other hand, offers the possibility of harnessing local minima in the risk function that are constrained to lie in appropriate "smoother" subsets of the full monotone set. This may be achieved in two distinct ways that may be combined for greater effect.

• 1.  
  By imposing one of the additional constraints (Equation (18)) on the EDoF of the fit. Examples of such restricted-freedom monotone fits are the red curves in Figures 3 and 4.
• 2.  
  By truncating the expansion (Equation (3)) to p number of coefficients (p < N) and then performing monotone risk minimization over this subset of the full monotone set.

In practice, such a smoother restricted-freedom monotone fit can be obtained by gradually reducing the value of q (Equation (18)) starting from the EDoF of the NSS fit until all low-amplitude, high-frequency wiggles in the fit disappear. Generally, the resulting fit has a lower risk than the NSS fit with EDoF = q and is manifestly consistent with trends captured by the full-freedom monotone fit. We find it useful to present (or consider) all three fits (NSS, full-freedom monotone, and restricted-freedom smoother monotone) together; this helps build a realistic picture about estimated trends in the data and thereby about the shape of the underlying true spectrum. Like the NSS fit, the smoother restricted-freedom monotone fit will generally be more biased than the full-freedom monotone fit. This greater bias, however, is partially compensated for by a larger risk which results in a larger confidence radius value (Equation (12)), a larger confidence set, and therefore more conservative inferences.

In this paper, we need to probe the confidence set for a fit for two purposes: (1) for validating cosmological models against a nonparametric fit (see Section 3.2) and (2) for finding the uncertainties on specific features of the fit such as peak heights and locations, and below we describe our method to scan and sample the confidence set for the latter. Our particular method for probing the confidence set for determining uncertainties on features of the fit is based on the following observations.

• 1.  
  The confidence set $\mathcal {D}_{N,\alpha }$ on the vector of coefficients (β₀, ..., β_{N − 1}), by construction, is centered at the vector of estimates $(\widehat{\beta }_0, \ldots, \widehat{\beta }_{N-1})$.
• 2.  
  The confidence interval defined in Equation (16) requires locating extreme variations in any functional T of the power spectrum f; e.g., location or height of a peak. The largest possible variations in T will be located as far away from the center of the confidence set as possible, i.e., on its surface.
• 3.  
  Cosmologically meaningful and sufficiently smooth variations in the fitted spectrum $\widehat{C}_l$ are most likely to be located in the projection of the full confidence set onto the lowest M dimensions.

We therefore generate a uniform sample from the projection of the full confidence set surface onto the lowest d dimensions, where 2 ⩽ d ⩽ M, with M ≲ 23. For convenience, we use the smoother restricted-freedom monotone fit for this purpose, with the justification that the confidence set corresponding to the full-freedom monotone fit happens to be nested inside that for this fit, for all four data realizations. The Cholesky factorization W = u^Tu is used to transform the original confidence ellipsoid $\mathcal {D}_{N,\alpha }$ (whose principal axes may not be aligned with coordinate directions in the β-space) into a sphere of the form $\lbrace \psi : \Vert \psi - \widehat{\psi } \Vert ^2 \le r^2_\alpha \rbrace$, where ψ = uβ and $\widehat{\psi } = u \widehat{\beta }$. The surface of this ψ-sphere can be efficiently sampled with uniform density using a standard algorithm (Section 3.4.1.E.6 on p.130 of Knuth 1981), and then transformed back into β-space, preserving uniformity of density because of the linearity of the transformation. From a sufficiently large sample of such variations of the power spectrum, we further selected those functions for which successive peaks and dips are separated by at least 50 multipole moments l. This cosmologically motivated selection criterion ensures that (1) the sampled functions are sufficiently smooth, and (2) high-frequency wiggles are not counted as peaks or dips when estimating uncertainties on locations and heights of peaks and dips (see Section 3.3). Based on cosmological considerations, we restricted the search to functions with three peaks (WMAP 1-, 3-, and 5-year data) or four peaks (7-year data). The set of functions thus sampled is used to estimate uncertainties on specific features of the fit. As an aside, we note that the confidence set construct and the formal guarantees related to confidence intervals (Equation (14)) do not necessarily imply a uniform density over the confidence set; uniform sampling is used here as a convenient computational device for scanning the confidence set surface in an unbiased fashion.

Considering that the noise in all four data sets is highly heteroskedastic and noise levels are especially high for large l, it would be useful to make an assessment of how such noise in the data affects local uncertainties in the fitted spectrum and to quantify how well the angular power spectrum value at each l is determined by the data.

To this end, we compute, for each data set, the approximate 95% confidence interval on each $\widehat{C}_l$ using 5000 function variations from the confidence set as outlined in Section 2.5. The length of this vertical confidence interval at given l, divided by the absolute value of the fit, $\vert \widehat{C}_l \vert$, provides an approximate indication of how well each $\widehat{C}_l$ is determined via the following interpretation: a value ≪ 1 indicates that the fit is well determined by the data, whereas values ≳ 1 imply that the data contain little or no information about the height of the power spectrum. This approach, which is inspired by the boxcar probe approach of Genovese et al. (2004), has the practical advantage of not having adjustable parameters (i.e., the boxcar width) in the procedure.

In Figure 5, we plot this height, scaled by the value of the fit, as a function of the multipole index l, for all four data realizations. We see that the range of l-values over which the fit is well determined expands consistently between 1-, 3-, and 5-year data realizations, from l ≈ 546 (1-year), to l ≈ 667 (3-year), to l ≈ 804 (5-year). On the other hand, while the l-range of the data expanded substantially between WMAP 5 and 7, the information contained in the data does not appear to have grown proportionately beyond l ≈ 842 for the 7-year data.

Zoom In Zoom Out Reset image size

In each of the four plots in Figures 3 and 4, we have also included parametric fits based on the ΛCDM (Hinshaw et al. 2003b, 2007; Nolta et al. 2009; Larson et al. 2011) and HΛCDM (Primack et al. 1995; Primack & Gross 2001) models (see the figure caption for details). The parametric ΛCDM-based fits turn out to be quite close to the respective nonparametric fits wherever the data are precise. This is remarkable considering that our nonparametric method does not rely on any cosmologically motivated prior information whatsoever. Moreover, the parametric ΛCDM-based fit appears to get closer to the respective nonparametric fit across the four WMAP data releases. The closeness of a parametric fit $(C_{l_{{\rm min}}}, \ldots, C_{l_{{\rm max}}})$ to the corresponding nonparametric fit $(\widehat{C}_{l_{{\rm min}}}, \ldots, \widehat{C}_{l_{{\rm max}}})$ can be measured through the distance

$$d(C,\widehat{C}) \approx \sqrt{ {1 \over N} \sum _{l=l_{{\rm min}}}^{l_{{\rm max}}} \left(C_l - \widehat{C}_l \over \sigma _l \right)^2 }.$$

Using Equation (12), this distance can be further expressed as the smallest confidence level α beyond which the parametric fit is rejected as a candidate for the true spectrum.

Table 1 lists distances of parametric fits based on the two cosmological models (ΛCDM and HΛCDM) from the nonparametric full-freedom monotone fit for the corresponding data realization (see the caption for specific details and description). The progression of distance values between a parametric ΛCDM fit and the corresponding nonparametric fit clearly shows that the two are getting closer as the data become precise. In contrast to this, the angular power spectrum generated by the particular HΛCDM model considered, which is almost as strong a contender for the true power spectrum as the ΛCDM fit with respect to the 1-year confidence set, is progressively pushed away to the boundary of the 95% confidence set for the 7-year confidence set. Visually, this trend can be understood on the basis of the differences between the HΛCDM power spectrum and the nonparametric fit (e.g., differences in the heights of the first peak; see Figures 3 and 4) that result in pushing this particular HΛCDM model out of the confidence set. Given the formal guarantees of this methodology (Section 2.2), the WMAP 7-year data thus rules out, at ≈95% confidence level, the particular HΛCDM model considered.

We now consider the problem of determining uncertainties on the locations and heights of peaks and dips in the nonparametrically fitted spectrum. The motivation for this exercise comes from the fact that peak locations and heights contain valuable information about cosmological models and parameters (Doran & Lilley 2002; Durrera et al. 2003).

Following the prescription outlined in Section 2.5, we sampled a set of 5000 function variations from the confidence set for each data realization. Peak and dip locations and heights were recorded for each peak and dip over this set of functions. This results in an empirical scatterplot that is indicative of the joint distribution of location and height for each peak or dip, under the assumption of uniform surface density on the confidence set $\mathcal {D}_{N,\alpha }$.

Figure 6 shows the results of probing the 95% confidence sets for uncertainties on peaks and dips, as outlined above, with 5000 acceptable function variations for each data realization. The box around a peak or a dip represents the largest horizontal and vertical variations in the scatter. In accordance with the confidence interval defined in Equation (16), these form the 95% confidence intervals on the location and height of a peak/dip. Table 2 lists these confidence intervals together with 95% confidence intervals on peak height ratios.

Zoom In Zoom Out Reset image size

Table 2. 95% Confidence Intervals on Several Features of the Angular Power Spectrum

Data	Peak Location	Peak Height	Dip Location	Dip Height	Peak Location Ratio	Peak Height Ratio
1-year	$l_1: \hfill (186,252)$	$h_1: \hfill (4968, 6133)$	$l_{1+{1 \over 2}}: \hfill (378,507)$	$h_{1+{1 \over 2}}: \hfill (1152, 2029)$	⋅⋅⋅	⋅⋅⋅
	$l_2: \hfill (440,680)$	$h_2: \hfill (1879, 3933)$	$l_{2+{1 \over 2}}: \hfill (494,772)$	$h_{2+{1 \over 2}}: \hfill (-2929, 2868)$	$l_2/l_1: \hfill (1.953,3.217)$	$h_2/h_1: \hfill (0.329,0.724)$
	$l_3: \hfill (559,897)$	$h_3: \hfill (41, 10733)$	$l_{3+{1 \over 2}}: \hfill (639,900)$	$h_{3+{1 \over 2}}: \hfill (-11376,9962)$	$l_3/l_1: \hfill (2.495,4.385)$	$h_3/h_1: \hfill (0.0075,1.906)$
3-year	$l_1: \hfill (184,254)$	$h_1: \hfill (5130, 6199)$	$l_{1+{1 \over 2}}: \hfill (384,456)$	$h_{1+{1 \over 2}}: \hfill (1402, 1917)$	⋅⋅⋅	⋅⋅⋅
	$l_2: \hfill (479,606)$	$h_2: \hfill (2172, 2997)$	$l_{2+{1 \over 2}}: \hfill (593,856)$	$h_{2+{1 \over 2}}: \hfill (80, 2145)$	$l_2/l_1: \hfill (2.028,3.000)$	$h_2/h_1: \hfill (0.380,0.555)$
	$l_3: \hfill (646,977)$	$h_3: \hfill (1506, 6045)$	$l_{3+{1 \over 2}}: \hfill (729,1000)$	$h_{3+{1 \over 2}}: \hfill (-4501, 5113)$	$l_3/l_1: \hfill (2.761,4.801)$	$h_3/h_1: \hfill (0.264,1.082)$
5-year	$l_1: \hfill (187,251)$	$h_1: \hfill (5249, 6301)$	$l_{1+{1 \over 2}}: \hfill (387,445)$	$h_{1+{1 \over 2}}: \hfill (1489, 1934)$	⋅⋅⋅	⋅⋅⋅
	$l_2: \hfill (485,596)$	$h_2: \hfill (2306, 2955)$	$l_{2+{1 \over 2}}: \hfill (608,793)$	$h_{2+{1 \over 2}}: \hfill (971, 2095)$	$l_2/l_1: \hfill (2.040,2.963)$	$h_2/h_1: \hfill (0.381,0.525)$
	$l_3: \hfill (663,982)$	$h_3: \hfill (1863, 4635)$	$l_{3+{1 \over 2}}: \hfill (731,1000)$	$h_{3+{1 \over 2}}: \hfill (-2978, 3157)$	$l_3/l_1: \hfill (2.883,4.672)$	$h_3/h_1: \hfill (0.321,0.833)$
7-year	$l_1: \hfill (187,252)$	$h_1: \hfill (5177, 6377)$	$l_{1+{1 \over 2}}: \hfill (390,442)$	$h_{1+{1 \over 2}}: \hfill (1512,1931)$	⋅⋅⋅	⋅⋅⋅
	$l_2: \hfill (492,589)$	$h_2: \hfill (2328, 3015)$	$l_{2+{1 \over 2}}: \hfill (623,803)$	$h_{2+{1 \over 2}}: \hfill (1074,2063)$	$l_2/l_1: \hfill (2.060,2.887)$	$h_2/h_1: \hfill (0.386, 0.534)$
	$l_3: \hfill (681,965)$	$h_3: \hfill (1871, 4119)$	$l_{3+{1 \over 2}}: \hfill (733,1104)$	$h_{3+{1 \over 2}}: \hfill (-3111, 2709)$	$l_3/l_1: \hfill (3.000,4.553)$	$h_3/h_1: \hfill (0.323, 0.732)$
	$l_4: \hfill (815,1193)$	$h_4: \hfill (-1391,7726)$	$l_{4+{1 \over 2}}: \hfill (951,1200)$	$h_{4+{1 \over 2}}: \hfill (-11102,6364)$	$l_4/l_1: \hfill (3.606,6.047)$	$h_4/h_1: \hfill (-0.244, 1.340)$

Notes. Negative values for heights and height ratios are a reflection of the high noise at the high-l end. Here, l_m (h_m) stands for the location (height) of the mth peak, and $l_{m+{1 \over 2}}$ ($h_{m+{1 \over 2}}$) denotes the location (depth) of the mth dip.

Download table as:  ASCIITypeset image

The following features of these results are worth pointing out. As is well known, the first peak was very clearly resolved in the 1-year data itself. Our results are manifestly consistent with this observation, in the sense that its box does not overlap with any other box. However, the uncertainty on the first peak does not shrink appreciably across the four data realizations. Further, our results clearly indicate that the second peak is resolved cleanly only in the 5-year data, whereas the third and fourth peaks are not resolved completely even in the 7-year data.

Consider the following relationship (Hu et al. 2001; Doran & Lilley 2002) between the location l_m of the mth peak, the acoustic scale l_A, and the shift parameter ϕ_m:

$$\begin{equation} l_m = l_A(m - \phi _m). \end{equation} \tag{ 19 }$$

Substituting the end-points of the 95% confidence interval for the mth peak location, this relationship results in a hyperbolic band of allowed values in the l_A–ϕ_m plane. Such bands, derived from 95% confidence intervals on the first three peaks (Table 2), are shown in Figure 7. Additional information from other sources is required to constrain these bands to physically meaningful regions in the l_A–ϕ_m plane. For example, if we assume l_A = 300 (Page et al. 2003) then, based on the 7-year data, the 95% confidence intervals for ϕ_m will be ϕ₁: (0.1600, 0.3767), ϕ₂: (0.0367, 0.3600), ϕ₃: (− 0.2167, 0.7300). Conversely, additional constraints on ϕ_m could be used to generate a confidence interval on l_A. From a model-independent point of view, we note that the (l_A, ϕ_m) bands for different peaks m appear to overlap around ϕ_m ≈ 0 and 200 ≲ l_A ≲ 400. We interpret this as a nonparametric realization of the nearly harmonic structure of peaks in the CMB power spectrum.

Zoom In Zoom Out Reset image size

Another interesting feature in Figure 6 is the tiny but clearly observable scatter for the very first dip at the low-l end. This scatter corresponds to extreme power spectrum variations that reside on the surface of the 95% confidence set and have an upturn at low l values. In the ΛCDM cosmology, such upturn at the low-l end is primarily the result of the integrated Sachs–Wolfe (ISW) effect and is seen in all parametric ΛCDM fits in Figures 3 and 4. It would therefore be interesting to see what could be said about the low-l upturn (and thereby about the ISW effect) based on the nonparametric confidence set.

Note that our nonparametric fits, which are at the center of their respective confidence sets, do not show a low-l upturn. However, the 7-year parametric ΛCDM fit, e.g., does show a clear upturn at the low-l end. This parametric fit is at a distance corresponding to a confidence level of about 10% (≈0.12σ; see Table 1) from our 7-year nonparametric full-freedom monotone fit. This means that the confidence set for the 7-year nonparametric fit contains spectra with a low-l upturn at most as far away as the 7-year parametric fit. We therefore conclude conservatively that the low-l upturn as a feature of the CMB angular power spectrum cannot be ruled out at any confidence level in excess of about 10% Actually, there are indications in our results (not shown) that such upturned variations of the spectrum may be much closer to the center of the confidence set for the 7-year full-freedom monotone fit; this needs further investigation.