A MARKOV CHAIN MONTE CARLO ALGORITHM FOR ANALYSIS OF LOW SIGNAL-TO-NOISE COSMIC MICROWAVE BACKGROUND DATA

We begin by assuming that the observed data may be modeled by a signal and a noise term,

$$\begin{equation} {\bf d}= {\bf A}{\bf s}+ {\bf n}, \end{equation} \tag{ 1 }$$

where d is a vector containing the data (at every pointing of the detectors), the matrix A involves both pointing and beam convolution (and where for this paper we will assume symmetric beams and neglect the details of this operation), and n is additive noise (here in the pixel domain). We assume both the CMB signal and noise to be Gaussian random fields with vanishing mean and covariance matrices S and N, respectively. In harmonic space, where s = ∑_ℓ,ma_ℓmY_ℓm, the CMB temperature covariance matrix is given by $\hbox{C}_{\ell m, \ell ^{\prime } m^{\prime }} = \langle a_{\ell m}^* a_{\ell ^{\prime } m^{\prime }}\rangle = C_{\ell } \delta _{\ell \ell ^{\prime }} \delta _{m m^{\prime }}$, C_ℓ being the angular power spectrum. A generalization to polarization merely requires the replacement of the signal matrix diagonal elements with 3 × 3 matrices of the form

$$\begin{equation} {\bf C}_{l} = \left[\begin{array}{ccc} C_{l}^{\rm TT} & C_{l}^{\rm TE} & C_{l}^{\rm TB} \\ C_{l}^{\rm ET} & C_{l}^{\rm EE} & C_{l}^{\rm EB} \\ C_{l}^{\rm BT} & C_{l}^{\rm BE} & C_{l}^{\rm BB}\end{array} \right]. \end{equation} \tag{ 2 }$$

For the discussion in this section, we focus on the temperature case, but note that the generalization to polarization is straightforward and discussed by Larson et al. (2007).

Given these assumptions, our goal is to quantify what has been learned about the underlying power spectrum of the CMB given the data, or how well the data constrain the cosmological parameters. One proceeds then, in a Bayesian framework, by writing down the posterior given the data,

$$\begin{equation} P(C_{\ell } | {\bf d}) \propto \mathcal {L}({\bf d}| C_{\ell }) P(C_{\ell }). \end{equation} \tag{ 3 }$$

Here, $\mathcal {L}({\bf d}| C_{\ell })$ is the likelihood and P(C_ℓ) is a prior on C_ℓ.

In order to derive the functional form of the likelihood, one imagines randomly choosing any relevant model (here a power spectrum drawn from P(C_ℓ)), and asks what sequence of effects needs to be modeled in order to simulate the data. Here, simulation is understood as conditioning on the chosen model, and leads to a joint density

$$\begin{eqnarray} P({\bf d},{\bf s},C_{\ell }) & = & P({\bf d},{\bf s}| C_{\ell }) P(C_{\ell }) \nonumber \\ & = & P({\bf d}| {\bf s}) P({\bf s}| C_{\ell }) P(C_{\ell }) \end{eqnarray} \tag{ 4 }$$

where the last line follows directly from our data model through the assumption of additive noise. Specifically, the factors in the above are

$$\begin{eqnarray} -2 \log P({\bf s}| C_{\ell }) & = & {\bf s}^{t} {\bf C}^{-1} {\bf s}- \log |{\bf C}| \nonumber \\ -2 \log P({\bf d}| {\bf s}) & = & -({\bf d}-{\bf s})^{t} {\bf N}^{-1}({\bf d}-{\bf s}) - \log |{\bf N}| \end{eqnarray} \tag{ 5 }$$

which follow from the assumption that both the signal and noise are independent Gaussian processes.

The idea of a "simulation chain" provides a conceptually clear approach to constructing a joint density, from which we immediately have the Bayesian posterior

$$\begin{equation} P(C_{\ell } | {\bf d}) = \int d{\bf s}\ P(C_{\ell }, {\bf s}| {\bf d}). \end{equation} \tag{ 6 }$$

The relevance of the above for this paper lies in relating what we refer to as the joint posterior, P(C_ℓ, s|d), and the more familiar likelihood $\mathcal {L}({\bf d}| C_{\ell }) \propto P(C_{\ell } | {\bf d}) / P(C_{\ell })$.

Although we can analytically compute the integral of the joint posterior over the signal for the Gaussian signal and noise processes considered here, and therefore simply write down the functional form of the likelihood, it is too expensive to evaluate it for any specified C_l given high-resolution data. Furthermore, for more complicated data models (i.e., including foreground model uncertainties) we will not be able to perform the integrals over the additional degrees of freedom. Both situations then instead motivate sampling from the joint posterior, and thereby generating samples from P(C_ℓ|d) without ever evaluating P(C_ℓ|d). We now discuss the original Gibbs sampling approach proposed and implemented by Jewell et al. (2004), Wandelt et al. (2004), and Eriksen et al. (2004), and then introduce a new MCMC step which directly addresses the previously reported slow probabilistic convergence in the low signal-to-noise regime (Eriksen et al. 2004).

As stated above, our goal is to sample from the joint posterior

$$\begin{equation} - 2 \log P({\bf s}, C_{\ell }|{\bf d}) = \chi ^{2}({\bf d}, {\bf s}) + {\bf s}^{t} {\bf S}^{-1} {\bf s}+ \log | {\bf S}| + \log P(C_{\ell }). \end{equation} \tag{ 7 }$$

For notational convenience, we have here dropped constant factors of 2π, and also defined

$$\begin{equation} \chi ^{2}({\bf s}, {\bf d}) = ({\bf d}- {\bf s})^{t} {\bf N}^{-1}({\bf d}-{\bf s}). \end{equation} \tag{ 8 }$$

One approach to sample from this posterior is to use an algorithm known as Gibbs sampling, where we can alternately sample from the respective conditional densities,

$$\begin{eqnarray} {\bf s}^{i+1} &\leftarrow P({\bf s}| C_{\ell }^i, {\bf d}) \\ C_{\ell }^{i+1} &\leftarrow P(C_{\ell } | {\bf s}^{i+1}, {\bf d}). \end{eqnarray} \tag{ 9 }$$

Here, ← indicates sampling from the distribution on the right-hand side. After some burn-in period, during which all samples must be discarded, the joint samples (sⁱ, Cⁱ_ℓ) will be drawn from the desired density. Thus, the problem is reduced to that of sampling from the two conditional densities P(s|C_ℓ, d) and P(C_ℓ|s, d).

We now describe the sampling algorithms for each of these two conditional distributions, starting with P(C_ℓ|s, d). First, note that P(C_ℓ|s, d) = P(C_ℓ|s) which follows directly from the construction of the joint density of "everything" above. This is also intuitively easy to understand since if we already know the CMB sky signal, the data themselves tell us nothing new about the CMB power spectrum. Next, since the sky is assumed to be Gaussian and isotropic, the distribution reads

$$\begin{equation} P(C_{\ell } | {\bf s}) \propto P(C_{\ell }) \frac{e^{-\frac{1}{2} {\bf s}_{\ell }^{t}{\bf S}_{\ell }^{-1}{\bf s}_{\ell }}}{\sqrt{|{\bf S}_{\ell }|}} = P(C_{\ell }) \frac{e^{-\frac{2\ell +1}{2} \frac{\sigma _{\ell }}{C_{\ell }}}}{C_{\ell }^{\frac{2\ell +1}{2}}}, \end{equation} \tag{ 11 }$$

which, when interpreted as a function of C_ℓ, is known as the inverse Gamma distribution. In this expression, $\sigma _{\ell } = \frac{1}{2\ell +1} \sum _{m} |a_{\ell m}|^2$ denotes the observed power spectrum of s. Fortunately, there exists a simple textbook sampling algorithm for this distribution (e.g., Gupta & Nagar 2000), and we refer the interested reader to the previous papers for details. For an alternative, and more flexible, sampling algorithm, see Eriksen & Wehus (2009).

In order to describe the sky signal sampling step, we first define the mean-field map (or Wiener filtered data) to be $\hat{{\bf s}} = ({\bf S}^{-1} + {\bf N}^{-1})^{-1} {\bf N}^{-1} {\bf d}$, and note that the conditional sky signal density given the data and C_l can be written as

$$\begin{eqnarray} P({\bf s}| C_{\ell }, {\bf d}) &\propto e^{-\frac{1}{2} ({\bf s}-\hat{{\bf s}})^t ({\bf S}^{-1} + {\bf N}^{-1}) ({\bf s}-\hat{{\bf s}})}. \end{eqnarray} \tag{ 12 }$$

Thus, P(s|C_ℓ, d) is a Gaussian distribution with mean equals $\hat{{\bf s}}$ and a covariance matrix equals (S⁻¹ + N⁻¹)⁻¹.

Sampling from this Gaussian distribution is straightforward, but computationally somewhat cumbersome. First, draw two random white noise maps ω₀ and ω₁ with zero mean and unit variance. Then solve the equation

$$\begin{equation} [{\bf S}^{-1} + {\bf N}^{-1}] {\bf s}= {\bf N}^{-1}{\bf d}+ {\bf S}^{-\frac{1}{2}} \omega _0 + {\bf N}^{-\frac{1}{2}} \omega _1 \end{equation} \tag{ 13 }$$

for s. Since the white noise maps have zero mean, one immediately sees that $\langle {\bf s}\rangle = \hat{{\bf s}}$, while a few more calculations show that 〈ss^t〉 = (S⁻¹ + N⁻¹)⁻¹.

The problematic part about this sampling step is the solution of the linear system in Equation (13). Since this a ∼10⁶ × 10⁶ system for current CMB data sets, it cannot be solved by brute force. Instead, one must use a method called Conjugate Gradients (CG), which only requires multiplication of the coefficient matrix on the left-hand side, not inversion. For details on these computations, together with some ideas on preconditioning, see Eriksen et al. (2004).

As originally applied to high-resolution CMB data, the Gibbs sampling algorithm as described above has very slow convergence at the high-ℓ, low signal-to-noise part of the spectrum. The reason for the slow convergence is easy to understand in light of the above: when sampling from P(C_ℓ|s), the typical step size is given by cosmic variance at all angular scales. In the high signal-to-noise regime, cosmic variance dominates the noise variance, and we are able to explore the full width of the posterior in only a few Gibbs iterations. However, in the low signal-to-noise end, cosmic variance is far smaller than the posterior variance, and it takes a prohibitively long time to converge probabilistically. This problem of "slow mixing" of the Gibbs sampler is illustrated in Figures 1 and 2. The long correlation length starting at signal to noise of unity leads to extremely long run times in order to produce a reasonable number of uncorrelated samples.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Assume that we have defined a deterministic sampling scheme for s, and that our new CMB map is given by some function

$$\begin{equation} {\bf s}_{n+1} = F\big({\bf s}_{n}, C_{\ell }^{(n+1)}, C_{\ell }^{(n)}\big). \end{equation} \tag{ 14 }$$

Then the condition of detailed balance for our MCMC sampler requires that

$$\begin{equation} F^{-1}\big({\bf s}_{n+1}, C_{\ell }^{(n+1)}, C_{\ell }^{(n)}\big) = F\big({\bf s}_{n+1}, C_{\ell }^{(n)}, C_{\ell }^{(n+1)}\big), \end{equation} \tag{ 15 }$$

or, in other words, that the inverse function is given by exchanging the order of the spectra in the function F. One simple function which has this property is

$$\begin{equation} {\bf s}_{n+1} = \left(\frac{C_{\ell }^{(n+1)}}{C_{\ell }^{(n)}}\right)^{\frac{1}{2}} {\bf s}_{n}. \end{equation} \tag{ 16 }$$

The total proposal matrix is then

$$\begin{eqnarray} w\big(C_{\ell }^{(n+1)}, {\bf s}_{n+1} | C_{\ell }^{(n)}, {\bf s}_{n}\big) & = & w\big(C_{\ell }^{(n+1)} | C_{\ell }^{(n)}, {\bf d}\big) \nonumber \\ & & \delta \left({\bf s}_{n+1} - \left(\frac{C_{\ell }^{(n+1)}}{C_{\ell }^{(n)}}\right)^{-\frac{1}{2}} {\bf s}_{n} \right), \nonumber \end{eqnarray}$$

and the "reverse" proposal is

$$\begin{eqnarray} w\big(C_{\ell }^{(n)}, {\bf s}_{n} | C_{\ell }^{(n+1)}, {\bf s}_{n+1}\big) & = & w\big(C_{\ell }^{(n)} | C_{\ell }^{(n+1)}, {\bf d}\big) \nonumber \\ & & \delta \left({\bf s}_{n} - \left(\frac{C_{\ell }^{(n)}}{C_{\ell }^{(n+1)}}\right)^{-\frac{1}{2}} {\bf s}_{n+1} \right). \nonumber \end{eqnarray}$$

The condition of detailed balance including deterministic moves requires the consideration of some technical points which we leave to Appendix A. There we show that the full Metropolis–Hastings accept probability reads

$$\begin{eqnarray} A & = & \min \left[ 1, \frac{e^{- \chi ^{2}({\bf s}_{n+1}, {\bf d})}}{e^{-\chi ^{2}({\bf s}_{n}, {\bf d})}} \frac{w\big(C_{\ell }^{(n)} | C_{\ell }^{(n+1)}, {\bf d}\big)}{w\big(C_{\ell }^{(n+1)} | C_{\ell }^{(n)}, {\bf d}\big)} \right]. \end{eqnarray} \tag{ 17 }$$

The significance of the above is that we can make relatively large changes to the power spectrum in the low signal-to-noise regime, where N⁻¹ is getting small, since the χ² is affected only very mildly by changes in any low signal-to-noise mode.

We note the interesting point (discussed more completely in Appendix B) that if one changes variables in the joint posterior from CMB maps, s, to whitened maps, ${\bf x}= {\bf C}_{\ell }^{-\frac{1}{2}} {\bf s}$, and then Gibbs sample in the new variables (C_ℓ, x), the resulting accept probability for changes in C_l conditioned on x is numerically identical to the above. While at first there might seem to be a conceptual distinction between MCMC schemes in a change of variables or those involving partial deterministic proposals, it turns out that this numerical equivalence holds more generally, and in Appendix B we establish a relation between both approaches. We bring this point to the attention of the reader as there could be other deterministic proposal schemes or another change of variables which lead to improvements over the approach presented in this paper.

For the numerical demonstration of the MCMC algorithm presented in this paper, we use a simple symmetric Gaussian proposal, truncated at C_ℓ > 0 (or, for polarization, the region where the resulting CMB covariance matrix is positive definite), for the power spectrum,

$$\begin{equation} w\big(C_{\ell }^{(n+1)} | C_{\ell }^{(n)}, {\bf d}\big) \propto e^{-\frac{1}{2} \left(\frac{C_{\ell }^{n+1}-C_{\ell }^{n}}{\tau _{\ell }}\right)^2} I(C_{\ell }>0), \end{equation} \tag{ 18 }$$

where τ_ℓ is a measure of the typical step size taken between two samples. Note that because this proposal density is symmetric, the ratio of C_ℓ proposals cancels, and the acceptance probability is entirely determined by the change in χ².

It should be noted that while the above MCMC step satisfies detailed balance, it is not irreducible, in the sense that there is not a nonvanishing probability in reaching any state from any other state in a finite number of MCMC steps; the phases are unchanged in each MCMC step. However, alternating these steps with a traditional Gibbs sampling step gives a combined "two-step" MCMC algorithm which indeed is irreducible, and therefore provably converges to the joint posterior. Once again, the details are left to the Appendices for the interested reader.

A general advantage of the Gibbs sampler is the fact that it is free of tunable efficiency parameters. The same is not true for the Metropolis–Hastings MCMC algorithm; for satisfactory sampling performance, it typically has to be tuned quite extensively. In this section, we describe three specific features that helps in this task, namely (1) step size tuning, (2) slice sampling, and (3) binning.

First, we have to ensure that the step size of our Gaussian proposal density roughly matches the width of the target distribution, in order to maintain both a reasonable acceptance rate and high mobility (for a general discussion of these issues and technical considerations of how to tune proposals, we refer the interested reader to Roberts & Rosenthal 2001). We do this by performing an initial test run, producing typically a few hundreds C_ℓ samples, and compute the standard deviation of these samples for each ℓ. These are then adopted as the proposal widths for the main run, scaled by some number less than unity, typically between 0.05 and 0.5. For the initial test run, we approximate the posterior width by the noise variance alone,

$$\begin{equation} \tau _{\ell }^{2} = \frac{2}{2\ell +1} \frac{N_{\ell }}{b_{\ell }^{2}}, \end{equation} \tag{ 19 }$$

because the MCMC sampler is used only in the low signal-to-noise regime. In this expression, N_ℓ is the power spectrum of the instrumental noise alone, and b_ℓ is the product of the Legendre transform of the beam and the HEALPix window function.

Next, the Metropolis–Hastings MCMC algorithm is inefficient in spaces with too many free parameters. For this reason, we divide the power spectrum coefficients, C_ℓ, into subsets, each containing typically only 10–20 multipoles. Then we propose changes to one subset at a time, while keeping all other multipoles fixed. Finally, we loop over subsets, and thus effectively implement a multipole slice Gibbs sampler for the full power spectrum.

This is computationally feasible, because a single MCMC proposal only requires a single χ² evaluation, which has a computational cost of a single spherical harmonic transform. Since drawing a full sky map from P(s|C_ℓ, d) in the classical Gibbs sampling step requires $\mathcal {O}(10^{2})$ spherical harmonic transforms, we can indeed afford to perform many MCMC proposals for each Gibbs step, without dominating the total cost.

Nevertheless, for very high resolution analysis, it is often beneficial to bin several C_ℓ together, both in order to increase the signal to noise of the joint coefficient, and to decrease the number of parameters that needs to be sampled by MCMC. We implement this by defining a new binned spectrum, weighted by ℓ(ℓ + 1)/2π, as follows:

$$\begin{equation} C_{b} = \frac{1}{N_b}\sum _{\ell \in b} \frac{\ell (\ell +1)}{2\pi } C_{\ell }. \end{equation} \tag{ 20 }$$

Here, b = [ℓ_min, ℓ_max] denotes the current bin, and N_b = ℓ_max − ℓ_min + 1 is the number of multipoles within the bin. These new (and fewer) coefficients are then sampled with the above MCMC sampler, after which the original spectrum coefficients are given by

$$\begin{equation} C_{\ell } = \frac{2\pi }{\ell (\ell +1)} C_{b}. \end{equation} \tag{ 21 }$$

The high-ℓ temperature simulation is designed to mimic the five-year WMAP temperature data (Hinshaw et al. 2009) with one exception, namely that the noise is assumed spatially uniform, in order to facilitate analytic comparison. Specifically, the CMB realization was drawn from the best-fit ΛCDM model derived from WMAP alone (Komatsu et al. 2009), including multipoles up to ℓ_max = 1000, and then smoothed with the instrumental beam of the WMAP V1 differencing assembly, and pixelized at HEALPix⁸ resolution N_side = 512. Finally, uniform noise of σ₀ = 40 μK RMS was added to each pixel. This corresponds to an S/N of unity at ℓ ∼ 550, roughly similar to the five-year WMAP data. We analyze this simulation both with and without the WMAP KQ85 sky cut (Gold et al. 2009).

In both analyses, we adopted the Gaussian proposal density with tuned variances, as described above. We also bin the power spectrum in progressively wide bins, starting at ℓ = 600, to maintain a reasonable signal to noise per sampled power spectrum parameter. Ten bins were sampled jointly per proposal, while all others were kept fixed.

In the full-sky case, we produced a total of 31,800 samples over 60 chains, and in the cut-sky case a total of 6800 samples. The cost for producing one sample in the latter, and by far most expensive, set was 2.5 CPU hours, for a total of 17,000 CPU hours. The number of MCMC steps per Gibbs step was one in the former and 20 in the latter. (Since the signal sampler dominates the cut-sky Gibbs chain, one can perform more low S/N steps without slowing down the overall code significantly.) In addition to these two main sample sets, we also produced two longer chains with each 3500 samples for the full-sky case, both with and without the new MCMC step turned on, in order to compare the Markov chain correlation lengths before and after including the MCMC sampler.

We first consider the full-sky data set, and in Figure 1 we show a segment of each of the two longer chains for three selected multipole bins. The top panel shows ℓ = 600, which is the first bin to be sampled by MCMC, the middle panel shows ℓ = 732–742, where there is still some signal in the data, and, finally, the bottom panel shows ℓ = 855–1000, which is strongly noise dominated. Starting with the top panel, we see that the red curve (Gibbs+MCMC) scatters significantly faster than the black curve (Gibbs only), implying more efficient sampling. This trend becomes even stronger with lower signal to noise, until the last case, where the Gibbs-only chain essentially does not move at all, while the MCMC sampler does probe the full range. Note, however, that even the MCMC sampler has a significant correlation length in this range, and this implies that there is still some room for improvement to be made in defining our proposals.

Next, these considerations are quantified in Figure 2, where we plot the Markov chain correlation length as a function of distance in the chain, for six bins with and without the MCMC sampler. As first reported by Eriksen et al. (2004), we see that the Gibbs-only correlation length increases dramatically with decreasing signal to noise, rendering the algorithm essentially useless in this regime. However, we also see that the new MCMC step effectively resolves this issue, as the correlation length (here defined by having a correlation less than 0.2) now is less than ∼40 steps. This is a dramatic improvement, and makes the algorithm useful even in this range. Nevertheless, we once again point out that it is possible to make further improvements by establishing better proposal densities.

In Figure 3, we consider the convergence properties of the ∼30 k samples set, by computing the Gelman–Rubin statistic R (Gelman & Rubin 1992) as a function of ℓ. Typically, one recommends that R should be less than, say, 1.2 in order to claim convergence. We see that this holds everywhere for this sample set, and typically it is even less than 1.05. Note also the step at ℓ = 600, showing clearly the beneficial effect of the MCMC sampler.

Zoom In Zoom Out Reset image size

Next, in Figure 4, we compare the marginal distributions derived from this sample set with the analytic result,

Zoom In Zoom Out Reset image size

$$\begin{equation} P(C_{\ell }|{\bf d}) \propto \prod _{\ell \in b} \frac{e^{-\frac{2\ell +1}{2} \frac{\sigma _{\ell }^{{\rm S+N}}}{b_{\ell }^2 C_{\ell }+N_{\ell }}}}{\big(b_{\ell }^2 C_{\ell }+N_{\ell }\big)^{\frac{2\ell +1}{2}}}. \end{equation} \tag{ 22 }$$

Here, b = [ℓ_min, ℓ_max] indicates a given multipole bin, b_ℓ denotes the product of the instrumental beam and the HEALPix pixel window, and σ^S+N_ℓ is the power spectrum of the noisy data map. We see that the new algorithm reproduces the analytic distributions very well, and this verifies the overall method.

Finally, the cut-sky power spectrum with one-sigma confidence regions is shown in three panels in Figure 5, focusing on different ℓ-ranges, namely all ℓ, the S/N ∼ 1 transition region, and the low S/N region. This completes the high-ℓ temperature analysis validation.

Zoom In Zoom Out Reset image size

We now consider polarization analysis, and construct a new low-ℓ simulation for this purpose. This simulation does not mimic any planned experiment, but is rather designed to highlight the analysis method itself. Specifically, we drew a new CMB realization from the best-fit WMAP ΛCDM spectrum that includes a nonzero tensor contribution, including multipoles up to ℓ_max = 150, and convolved this with a 3° FWHM Gaussian beam, and pixelized it at N_side = 64. Uniform noise of 5 μK RMS was added to the temperature component, and 1 μK RMS to the polarization components. The five-year WMAP polarization sky mask was imposed on the data.

We allowed for nonzero C^TT_ℓ, C^TE_ℓ, C^EE_ℓ, and C^BB_ℓ spectra, but fixed C^TB_ℓ = C^EB_ℓ = 0. These spectra were then individually binned to maintain a reasonable signal to noise per bin. (Details on how to introduce individual binning of each power spectrum were recently described by Eriksen & Wehus 2009.) Again, a tuned Gaussian proposal density was used in the MCMC step. A total of 12,000 samples were produced over 12 chains, and the CPU time per sample was 55 s, for a total of ∼200 CPU hours.

In Figure 6, we show one C_ℓ chain for each of the four sampled spectra, for the last (and therefore most difficult) bin in each case. Note that the C^EE_ℓ and C^BB_ℓ spectra have essentially vanishing signal to noise, and therefore these chains reach zero values. Clearly, we see that mixing properties of these chains are satisfactory, and the correlation lengths are quite short.

Zoom In Zoom Out Reset image size

In Figure 7, we show the Gelman–Rubin statistics for each of the four power spectra, and with the single exception of the very last bin of C^EE_ℓ, all R values are well below 1.1. Thus, all spectra have converged well everywhere.

Zoom In Zoom Out Reset image size

Finally, in Figure 8, we show the reconstructed marginal power spectra for each polarization component, overplotted on the input spectrum. The agreement is very good. Note, however, that these spectra are direct marginals, and not a joint maximum likelihood estimate. They are therefore not individual unbiased estimators. In particular, the marginal C^EE_ℓ power spectrum is biased slightly high because of the combination of the C^TT_ℓC^EE_ℓ − (C^TE_ℓ)² > 0 positivity constraint and relatively low signal to noise. Consideration of the joint polarization posterior, which is an unbiased estimator, is postponed to a future publication.

Zoom In Zoom Out Reset image size