Bayesian power-spectrum inference with foreground and target contamination treatment

TL;DR

This paper introduces a Bayesian framework for jointly inferring cosmological power-spectra and 3D density fields from galaxy surveys, effectively handling foreground contamination and systematic uncertainties.

Contribution

It extends the ARES framework by incorporating block sampling for foreground and target contamination coefficients, enabling fully self-consistent and unbiased power-spectrum inference.

Findings

Successfully infers unbiased power-spectra from mock data with stellar contamination.

Accurately accounts for galaxy biases and noise levels in density field reconstruction.

Quantifies correlated uncertainties, revealing up to 10% correlations across Fourier space.

Abstract

This work presents a joint and self-consistent Bayesian treatment of various foreground and target contaminations when inferring cosmological power-spectra and three dimensional density fields from galaxy redshift surveys. This is achieved by introducing additional block sampling procedures for unknown coefficients of foreground and target contamination templates to the previously presented ARES framework for Bayesian large scale structure analyses. As a result the method infers jointly and fully self-consistently three dimensional density fields, cosmological power-spectra, luminosity dependent galaxy biases, noise levels of respective galaxy distributions and coefficients for a set of a priori specified foreground templates. In addition this fully Bayesian approach permits detailed quantification of correlated uncertainties amongst all inferred quantities and correctly marginalizes over observational systematic effects. We demonstrate the validity and efficiency of our approach in obtaining unbiased estimates of power-spectra via applications to realistic mock galaxy observations subject to stellar contamination and dust extinction. While simultaneously accounting for galaxy biases and unknown noise levels our method reliably and robustly infers three dimensional density fields and corresponding cosmological power-spectra from deep galaxy surveys. Further our approach correctly accounts for joint and correlated uncertainties between unknown coefficients of foreground templates and the amplitudes of the power-spectrum. An effect amounting up to $10$ percent correlations and anti-correlations across large ranges in Fourier space.

Figures

Equations

$$N_i=\bar{N}{R}_i(1+b{D}_i{\delta }_i)+{ϵ}_i,$$

$$⟨{ϵ}_i{ϵ}_j⟩=\bar{N}{R}_i{\delta }_{i,j}^K,$$

$$C_{i,\text{obs}}=\frac{N_{i,\text{spectro}}^g}{N_{i,\text{targets}}^g}=\frac{N_{i,\text{spectro}}^g}{N_{i,\text{photo}}^g}\frac{N_{i,\text{photo}}^g}{N_{i,\text{targets}}^g}=C_i{M}_i^{-1},$$

$$M_i=\prod _{n=0}^{N^{\text{fg}}}(1-{\alpha }_n{F}_{n,i}),$$

$$N_i^c={\bar{N}}^c{M}_i(\{{\alpha }_n^c\})R_i^c(1+b^c{D}_i{\delta }_i)+{ϵ}_i^c,$$

$$⟨{ϵ}_i^c{ϵ}_j^{c^{\prime }}⟩={\bar{N}}^c{M}_i(\{{\alpha }_n^c\})R_i^c{\delta }_{i,j}^K{\delta }_{c,c^{\prime }}^K.$$

$$\begin{gathered} \mathcal{L}(\{N_i^c\}|\{{\delta }_i\},\{{\bar{N}}^c\},\{b^c\})\propto \prod _{c=0}^{N^c}\prod _{i=0}^{N^v}{\left({\bar{N}}^c{M}_i(\{{\alpha }_n^c\})R_i^c\right)}^{-1/2} \\ \exp \{-\frac{1}{2}\frac{1}{{\bar{N}}^c{M}_i(\{{\alpha }_n^c\})R_i^c} \\ {\left[N_c^i-{\bar{N}}^c{M}_i(\{{\alpha }_n^c\})R_i^c(1+b^c{D}_i{\delta }_i)\right]}^2\}. \end{gathered}$$

$$\begin{gathered} \mathcal{P}\left(\{{\alpha }_n\}|\{N_i\},\{{\delta }_i\},\{\bar{N}\},\{b\}\right)\propto \\ \mathcal{P}\left(\{{\alpha }_n\}\right)\mathcal{L}(\{N_i\}|\{{\delta }_i\},\{\bar{N}\},\{b\},\{{\alpha }_n\}), \end{gathered}$$

$${\alpha }_0$$

$${\alpha }_1$$

$${\alpha }_2$$

$$⋮$$

$${\alpha }_{N^{fg}-1}$$

$$R_i^c=\frac{1}{|{\mathcal{V}}_i|}{\int }_{{\mathcal{V}}_i}{\text{d}}^3\text{@vec}x{C}^c(\hat{x})f^c(x),$$

$$N_i^c={\bar{N}}^c{R}_i^c{M}_i^c(1+b^c{\delta }_i)+\sqrt{{\bar{N}}^c{R}_i^c{M}_i^c}{ϵ}_i,$$

$$C(\theta {)}_n=\frac{1}{N-n}\sum _{i=0}^{N-n}\frac{{\theta }^i-⟨\theta ⟩}{\sqrt{\mathrm{Var}\left(\theta \right)}}\frac{{\theta }^{i+n}-⟨\theta ⟩}{\sqrt{\mathrm{Var}\left(\theta \right)}},$$

$$\log \mathcal{P}\left({\alpha }_k|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\}\right)=-\frac{1}{2}\sum _{i=0}^v\frac{{\left[N_i-\bar{N}{M}_i\left(\{{\alpha }_n\}\right)R_i(1+D_{i}b{\delta }_i)\right]}^2}{\bar{N}{M}_i\left(\{{\alpha }_n\}\right)R_i}-\frac{1}{2}\log \left(\bar{N}{M}_i\left(\{{\alpha }_n\}\right)R_i\right).$$

$${\tilde{R}}_i=\bar{N}{R}_i\prod _{n\setminus k}(1-{\alpha }_n{F}_{n,i}),$$

$$A_i={\tilde{R}}_i(1+D_{i}b{\delta }_i).$$

$$\log \mathcal{P}\left({\alpha }_k|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\}\right)=-\frac{1}{2}\sum _{i=0}^v\frac{{\left[N_i-(1-{\alpha }_k{F}_{k,i})A_i\right]}^2}{(1-{\alpha }_k{F}_{k,i}){\tilde{R}}_i}-\frac{1}{2}\log \left[(1-{\alpha }_k{F}_{k,i}){\tilde{R}}_i\right].$$

$$B_i$$

$$C_{k,i}$$

$${\gamma }_{k,i}$$

$$\log \left(\mathcal{P}\left({\alpha }_k|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\},\right)\right)=-\frac{1}{2}\sum _{i=0}^v\frac{{\left(B_i+{\alpha }_k{C}_{k,i}\right)}^2}{\left({\tilde{R}}_i-{\alpha }_k{\gamma }_{k,i}\right)}-\frac{1}{2}\log \left({\tilde{R}}_i-{\alpha }_k{\gamma }_{k,i}\right).$$

$$\forall i,{\tilde{R}}_i-{\alpha }_k{\gamma }_{k,i}>0$$

$$\forall i,{\tilde{R}}_i>0,$$

$${\tilde{R}}_i={\tilde{R}}_{i,k}^{\prime }+\omega {\gamma }_{k,i}.$$

$${\tilde{R}}_{i,k}^{\prime }$$

$$\omega$$

$$\omega =\underset{i}{\min }\left(\frac{{\tilde{R}}_i}{{\gamma }_{k,i}}\right).$$

$$\log \left(\mathcal{P}\left({\alpha }_k|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\},\right)\right)=-\frac{1}{2}\sum _{i=0}^v\frac{{\left(B_i+{\alpha }_k{C}_{k,i}\right)}^2}{{\tilde{R}}_{i,k}^{\prime }+\left(\omega -{\alpha }_k\right){\gamma }_{k,i}}-\frac{1}{2}\log \left({\tilde{R}}_{i,k}^{\prime }+\left(\omega -{\alpha }_k\right){\gamma }_{k,i}\right).$$

$$\mathcal{P}\left({\alpha }_k|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\},\right)\propto \prod _i{\int }_{-\infty }^{\infty }\mathrm{d}{t}_i\frac{{\mathrm{e}}^{-\frac{1}{2}{\sum }_{i=0}^v\frac{{\left(B_i+{\alpha }_k{C}_{k,i}-t_i\right)}^2}{\left(\omega -{\alpha }_k\right){\gamma }_{k,i}}}}{\sqrt{2\pi \left(\omega -{\alpha }_k\right){\gamma }_{k,i}}}\frac{{\mathrm{e}}^{-\frac{1}{2}\frac{t_i^2}{{\tilde{R}}_{i,k}^{\prime }}}}{\sqrt{2\pi {\tilde{R}}_{i,k}^{\prime }}}.$$

$$\xi =\omega -{\alpha }_k.$$

$$\mathcal{P}\left(\xi ,\{t_i\}|\{{\alpha }_n\}\setminus {\alpha }_k,\{N_i\},\{{\delta }_i\},\{\bar{N}\},\right)\propto \prod _i\frac{{\mathrm{e}}^{-\frac{1}{2}{\sum }_{i=0}^v\frac{{\left(B_i+\left(\xi -\omega \right)C_{ki}-t_i\right)}^2}{\xi {\gamma }_{k,i}}}}{\sqrt{2\pi \xi {\gamma }_{k,i}}}\frac{{\mathrm{e}}^{-\frac{1}{2}\frac{t_i^2}{{\tilde{R}}_{ik}^{\prime }}}}{\sqrt{2\pi {\tilde{R}}_{i,k}^{\prime }}},$$

Full text

1 Introduction

In recent years the cosmological community has witnessed great improvements in our understanding of the Universe. This progress is particularly due to the spectacular results of the Planck satellite mission and deep galaxy observations such as the ones provided by the Baryon Oscillation Sky Survey (Planck Collaboration et al. 2011; Dawson et al. 2013; Planck Collaboration et al. 2016). These results put high standards for future analyses of cosmological data with an ever increasing need to control uncertainties and systematic effects in observations in order not to misinterpret data when searching for cosmological signals. To address these needs, data science is challenged to provide ever more robust data models accounting for complex systematic effects and allowing for accurate marginalization over unknowns when interpreting cosmological data.

A particular challenge in existing and coming deep galaxy redshift surveys arises from the need to properly understand selection processes of galaxies from which cosmological surveys are constructed (see e.g. Huterer et al. 2013). Such identification was conducted for mitigating star-galaxy contamination of the first SDSS photometric galaxy catalogue (Scranton et al. 2002). The problem is further exacerbated by our lack of understanding of galaxies as tracers of the underlying dark matter field when performing cosmological inference. In particular all our indicators for completeness rely on the relative slow, homogeneous and isotropic evolution of galaxy densities relative to dark matter densities. If the observation is further hindered by instrumental and/or terrestrial effect, this leads to a complex and challenging analysis problem.

In particular, Ross et al. (2011) and Ho et al. (2012) have identified that contamination by bright stars alter significantly the intrinsic clustering signal of the observed photometric galaxy sample at large scales. The last SDSS release based on DR12 photometry still shows this problems in the measured correlation function (Ross et al. 2017). Effects due to foreground stars, dust, seeing, sky background intensity have the greatest potential to cause systematic deviations in the clustering signal (see e.g. Ho et al. 2015; Leistedt & Peiris 2014). Selection of spectroscopic galaxies is not affected immediately by the same effect, but it is subject to other systematics, such as fibre collisions, target priority conflicts and fibre plate fixations. All these data contaminations constitute a particular nuisance, since foreground affects are also affecting the noise properties of observed galaxy samples via varying attenuation or target contamination across the sky.

In the large scale structure community, foreground effects are traditionally treated by weighting observed galaxies to homogenize the distribution of the traced density field across the sky (see e.g. Ross et al. 2011; Sánchez et al. 2013; Ross et al. 2017). This approach neglects possible and sometimes counter intuitive correlations between contamination effects and power-spectrum amplitudes across large ranges in Fourier space. It also ignores the effects of modified observational noise properties due to target contamination. There have also been several methods proposed to account for additive contributions from unknown foregrounds in photometric and spectroscopic galaxy observations (see e.g. Tegmark et al. 1998; Leistedt & Peiris 2014; Ho et al. 2015). Also the literature on cosmic microwave background analyses provides a plenitude of approaches to account for linear additive foreground contributions (see e.g. Tegmark & Efstathiou 1996; Hinshaw et al. 2007; Eriksen et al. 2008; Sudevan et al. 2016; Vansyngel et al. 2016; Elsner et al. 2017).

However, all these approaches do not properly account for multiplicative foreground and target contaminations which also affect the noise. In this work we expand on the idea that foreground effects are more closely related to multiplicative rather than additive contributions. This is similar to the discussion presented in Huterer et al. (2013), who used a multiplicative correction in observed galaxy densities.

In this work we seek to account for these effects when inferring cosmological power-spectra from observations. Literature provides a plenitude of various statistically more or less rigorous approaches to measure power-spectra. Several of these methods rely on Fourier transform based methods or exploit Karhunen-Loève or spherical harmonics decompositions (see e.g. Feldman et al. 1994; Tegmark 1995; Hamilton 1997a; Yamamoto 2003; Percival et al. 2004; Tegmark et al. 1997, 2004; Pope et al. 2004; Fisher et al. 1994; Heavens & Taylor 1995; Tadros et al. 1999; Percival et al. 2004; Percival 2005). Other approaches aim at inferring the real space power spectrum via likelihood methods (Ballinger et al. 1995; Hamilton 1997a, b; Tadros et al. 1999; Percival 2005).

In the Bayesian community several approaches have been proposed to jointly infer three dimensional density fields and their corresponding cosmological power-spectra (see e.g. Jasche et al. 2010; Granett et al. 2012; Alsing et al. 2016). Also note, that similar approaches explored for analyses of cosmic microwave background data (see e.g. Wandelt et al. 2004; O’Dwyer et al. 2004; Eriksen et al. 2004; Jewell et al. 2004; Larson et al. 2007; Eriksen et al. 2007; Jewell et al. 2009).

To account for such effects of foreground and target contamination in a statistically rigorous fashion, we propose a hierarchical Bayesian approach to jointly and self-consistently infer three dimensional density fields, corresponding power-spectra and coefficients of a set of different foreground templates. In particular this work builds upon our previously developed Bayesian inference algorithm ARES (Algorithm for REconstruction and Sampling) (see e.g. Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015; Lavaux & Jasche 2016).

The manuscript is structured as follows. In Section 2 we give a brief overview of the statistical model that we propose. First we remind in Section 2.1 the hierarchical Bayesian inference approach on which our code, ARES, is based. Then we describe in Section 2.2 the necessary modifications of the model for foreground effects and in Section 2.3 the modification to the original inference algorithm. In Section 3 we describe the generation of artificial data used to test the performance of the sampling framework in section 4. Finally we discuss our results and give an outlook on future applications in Section 5.

2 The Bayesian inference model

This section provides a brief overview over our previously presented Bayesian inference framework ARES (see e.g. Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015). We will also give a detailed description of the modifications enabling us to account for foreground and target contamination in deep galaxy observations.

2.1 The ARES framework

As discussed in the introduction, the current work builds upon the previously developed Algorithm for Reconstruction and Sampling ARES (see e.g. Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015). This full Bayesian large scale structure inference method aims at precision inference of cosmological power-spectra from galaxy redshift surveys. Specifically it performs joint inferences of three dimensional density fields, cosmological power spectra as well as luminosity dependent galaxy biases and corresponding noise levels for different galaxy populations in the survey (Jasche et al. 2010; Jasche & Wandelt 2013; Lavaux & Jasche 2016).

In particular the ARES algorithm addresses a large scale data interpretation problem involving many millions of parameters. In the following, for the sake of clarity, we will describe the corresponding data model in case of a single galaxy population. The generalization of this data model to account for an arbitrary number of galaxy populations with respective stochastic and systematic uncertainties has been described in our previous works (Jasche & Wandelt 2013; Lavaux & Jasche 2016). In case of a single galaxy population the data model is given as:

[TABLE]

where $N_{i}$ is the number of galaxies in the $i$-th grid element, $\bar{N}$ is the mean density of the galaxy population, $R_{i}$ is the overall linear response operator of the survey, describing redshift and target completeness, $b$ is the galaxy population bias, $D_{i}$ is the cosmic growth factor at the position of $i$-th grid element, $\delta_{i}$ is the density contrast at a reference redshift in this same grid element and $\epsilon_{i}$ denotes random but structured instrumental noise. Also, as described in our previous work, the observational noise will be assumed to be a Gaussian approximation to Poissonian noise, neglecting the influence of the signal itself. This assumption yields the corresponding noise covariance matrix as:

[TABLE]

with $\delta^{K}_{i,j}$ being the Kronecker-Delta (i.e. equal to one if $i=j$ and zero otherwise). Finally, we assume a homogeneous and isotropic Gaussian prior for density contrast amplitudes $\delta_{i}$. For further details please consult our previous work Jasche & Wandelt (2013).

In order to provide full Bayesian uncertainty quantification the algorithm explores the joint parameter space of density fields, power-spectra, galaxy biases and noise parameters via an efficient block sampling scheme. We show in Figure 1 a visualization of this iterative sampling procedure including the foreground sampling method presented in this work. Iterative execution of these respective sampling steps provides us with a valid Markov chain and a numerical representation of the full joint target posterior distribution. Also note, we use here an upgraded version of the ARES algorithm for which we employ the messenger method discussed in Elsner & Wandelt (2013). This particular implementation of the Wiener posterior sampling method has been demonstrated to improve upon the statistical efficiency of previous implementations (Jasche & Lavaux 2015).

2.2 The foreground and target contamination model

Spectroscopic completeness is generally computed by the ratio of the number of observed spectra and the number of all photometric targets for a given area in the sky. This ratio is assumed to hold for any pointing of this area. However besides galaxies also a number of unknown contamination can contribute or affect observed photometric targets, artificially increasing or depleting their local number density. This contamination may include e.g. foreground stars, dust absorption or effects due to seeing. Naive estimate of the spectroscopic completeness from data therefore does not reflect the actual probability of obtaining a galaxy spectrum at a given position in the sky. From observations we can build an estimate of the completeness by calculating the ratio of the number of observed galaxy spectra $N^{g}_{i,\text{spectro}}$ and the number of target galaxies $N^{g}_{i,\text{targets}}$, both in the direction of the pixel $i$ in the sky:

[TABLE]

with $N^{g}_{i,\text{photo}}$ being the actual true number of galaxies that should have been identified by photometry and $M_{i}$ is the ratio of all the observed photometric targets $N^{g}_{i,\text{targets}}$ to the true sample of target galaxies $N^{g}_{i,\text{photo}}$. Equation (3) demonstrates the dilemma, that given some spectroscopic information (for which objects can clearly be identified as galaxies) there is no immediate way to decide how many galaxies were in the actual target sample. In Equation (3) this mismatch is quantified by the ratio $M_{i}=N^{g}_{i,\text{targets}}/N^{g}_{i,\text{photo}}$ between the number of real photometric targets that should have been chosen vs the actual but unknown number of all galaxy targets.

We expect two possible contributions to $M_{i}$: either there is an excess in the number of targets because the photometric information was insufficient to separate galaxies from stars or other objects in the sky, or there is a lack of galaxies that have not been detected due to, e.g., low surface brightness or dust absorption. Assuming such foreground contributions to be mild perturbations, the corresponding contamination map $M$ on the sky can be expressed as a product of small contaminants (greater or lesser than one). Individual contaminations can then be modelled via respective foreground templates. Note that this approach has also been adopted by the BOSS collaboration to correct their measurement of the cosmological power-spectrum (Ross et al. 2011; Ho et al. 2012; Anderson et al. 2012). Given this assumption we can express individual $M_{i}$ as:

[TABLE]

where $F_{n,i}$ is the foreground template of the $n$-th contribution at the $i$-th pixel of the map, $\alpha_{n}$ is the amplitude of respective foreground templates, $N^{\text{fg}}$ is the total number of foreground maps. We note that different surveys or even different sub-samples of observed galaxies may be subjected to different foreground effects. To consistently account for all these foreground effects when jointly analysing individual or several data sets the original data model implemented in ARES needs to be modified by a multiplicative correction of the survey response operator $R^{c}_{i}$. Specifically, we model the observed number of galaxies $N^{c}_{i}$ in a survey as:

[TABLE]

where the additive noise contribution $\epsilon^{c}_{i}$ drawn from a zero mean Gaussian distribution with covariance matrix given as:

[TABLE]

The superscript $c$ indicates different considered catalogues. These modifications render the posterior distribution of $\delta_{i}$ more complex. By construction the likelihood is given as:

[TABLE]

It should be remarked that foreground contributions, as modelled here, are not just mere additive contributions to the signal to infer, but they also have pronounced impact on the varying noise properties across the survey. Hence, as can be seen from Equation (7) inferring the foreground coefficients $\{\alpha_{n}\}$ is a highly non-linear analysis task.

2.3 Sampling foreground coefficients

As described by the likelihood distribution given in Equation (7) foregrounds of different catalogues as labelled by the superscript $c$ can be sampled independently. For this reason and without loss of generality here we provide the sampling procedure for a single galaxy catalogue only. Then the conditional posterior distribution for the coefficients $\{\alpha_{n}\}$ of respective foreground templates can be written as:

[TABLE]

where $\mathcal{P}\left(\{\alpha_{n}\}\right)$ is the prior of foreground coefficients and $\mathcal{L}(\{N_{i}\}|\{\delta_{i}\},\{\bar{N}\},\{b\},\{\alpha_{n}\})$ is the likelihood for a single catalogue as given by Equation (7). In the absence of any further information on the amplitudes of foreground coefficients $\alpha_{n}$ we follow the maximal agnostic approach by using a uniform prior $\mathcal{P}\left(\{\alpha_{n}\}\right)=1$. It can be seen from Equations (8) and (7) that conditional posterior distribution does not factorize in the coefficients $\alpha_{n}$ for respective foreground templates. To correctly account for the conditional dependencies between the coefficients $\{\alpha_{n}\}$, we propose to use a block sampling procedure to sequentially draw random variates of $\alpha_{n}$ with respect to all other values. This is achieved by introducing the following sequence of sequential sampling steps to the full ARES framework:

[TABLE]

where the symbol ”$\sim$” indicates a random draw from respective distributions. This sampling procedure integrates well into the ARES framework as indicated in Fig. 1. Despite the fact that drawing respective realizations of the foreground coefficients $\{\alpha_{n}\}$ is a non-linear process, there exists a direct sampling procedure. The detailed derivation of the foreground coefficient sampler is presented in Appendix A. The detailed algorithm for generating respective random variates is given in Algorithm 1.

3 Generation of Gaussian mock data

To test the validity and performance of the modified ARES sampling framework we follow a similar approach as discussed in our previous works (Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015). These mock catalogues are generated in accordance with the data model described in Equation (7), including various foreground effects. We generate artificial galaxy data on a cubic equidistant grid of side length $4000\,h^{-1}\mathrm{Mpc}$ consisting of $256^{3}$ grid nodes.

First a realization of a cosmic density contrast field $\delta_{i}$ is drawn from a zero-mean normal distribution with covariance matrix corresponding to a cosmological power-spectrum. This spectrum, including baryon acoustic oscillations, is calculated according to the prescription described in Eisenstein & Hu (1998) and Eisenstein & Hu (1999). For numerical evaluation we assume a $\Lambda$CDM cosmology with the set of parameters given as ($\Omega_{m}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{b}=0.0485$, $h=0.6774$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$), as determined by Cosmic Microwave Background observations of the Planck satellite mission (Planck Collaboration et al. 2016).

Following a similar description as discussed in Jasche & Wandelt (2013), we intent to create a realistic scenario to jointly analyse the SDSS DR7 main, the CMASS and the LOW-Z galaxy sample while accounting for their respective systematic effects. Corresponding artificial data sets are then drawn from the distribution given in Equation (7). These artificial data sets include the effects of noise, galaxy bias, survey geometries, selection and foreground effects. For the sake of this work we restrict our tests to accounting for the dominant foreground effects of dust extinction and stellar contamination. The templates for these two effects are presented in the upper panels of Fig. 2 the lower panels show the completeness masks for the respective surveys.

The map describing dust extinction has been generated straightforwardly from the SFD maps (Schlegel et al. 1998). In particular we constructed a HEALPix map of the reddening at $N_{\text{side}}=2048$ by linearly interpolating the values of the SFD map (Górski et al. 2005). The star map is built from different pieces of information. The first component consists in computing a MANGLE description111MANGLE software is originally provided by Swanson et al. (2008). of the geometry of the spectroscopic plates. We then count the number of stars with apparent magnitudes $20.3<i_{\text{PSF}}<20.6$ present in each single non-overlapping polygon. We convert these MANGLE description into an HEALPix map and divide the value in each pixel by the area of the overlapping polygon. This results in a map for which each pixel has an estimated star count by steradians, thus a star density. We have chosen to average by spectroscopic plates to reduce shot-noise in the estimate. A better estimate would have been obtained from the geometrical description of the photometric tiling but at the cost of increased noise.

Further, for the SDSS DR7 main sample component of the mock data we assume a radial selection function following from a standard Schechter luminosity function with standard r-band parameters ($\alpha=-1.05$, $M_{*}-5\mathrm{log}_{10}(h)=-20.44$), and we limit the survey to only include galaxies within an apparent Petrosian r-band magnitude range $13.5<r<17.6$ and within the absolute magnitude ranges $M_{min}=-17.0$ to $M_{max}=-23.0$. As usual, the radial selection function $f(z)$ is then given by the integral of the Schechter luminosity function over the range in absolute magnitude.

For the CMASS and LOW-Z component we have used numerical estimates of the selection functions by computing a histogram of the corresponding $N(d)$ in the actual data sets (e.g. for DR12 Ross et al. 2017) ($d$ being the co-moving distance from the observer). To account for the different selection effects in the northern and southern galactic plane we also split our mock data sets into the CMASS and LOW-Z catalogues correspondingly. The respective radial selection functions are presented in Fig. 3. The average of the product of the two dimensional survey geometry $C^{c}(\hat{x})=C^{c}_{i(\hat{x})}$ and the selection function $f(x)$ at each grid element in the three dimensional volume yields the survey response operator:

[TABLE]

with $|\mathcal{V}|$ the volume of the set $\mathcal{V}$, $\mathcal{V}_{i}$ indicating the volume represented by the $i$-th grid element.

Given these definitions and a realization of the three dimensional density field $\delta_{i}$, realizations of artificial galaxy observations for respective catalogues labelled with $c$ can be obtained by evaluating:

[TABLE]

where $\epsilon_{i}$ is a white-noise field drawn from zero-mean and unit variance normal distribution.

4 Results

In this section we discuss results obtained by applying the modified ARES algorithm to artificial mock data. In particular in this work we focus on the validity and statistical efficiency of the algorithm.

4.1 Statistical efficiency of the sampler

To test the statistical efficiency of our sampler we follow a standard test procedure as described in previous works (see e.g. Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015). In particular we test the initial burn-in phase by starting the sampler from an over dispersed state and monitoring transitions in parameter space as a sequence of steps in the Markov chain. Typically this test reveals a coordinated drift of inference parameters towards their target values. Once the chain moved to a preferred region in parameter space it starts to correctly explore the target posterior distribution via the random walk Gibbs sampling approach. At this stage it is assumed that the sampler has passed the initial burn-in phase and we start recording samples of the Markov chain. In Fig. 4 we show the sequence of sampled posterior power-spectra during the burn-in phase. For this test we started from an over dispersed state by multiplying the initial guess of a Gaussian random density contrast field by a factor of $0.1$. Fig. 4 nicely demonstrates the initial drift towards the preferred region in parameter space. As can be seen the initial burn in phase consists of about $\sim 2000$ Markov steps. As a side remark we note that generation of individual samples require an investment of $\sim 1.87$ CPU-hours per sample for the present scenario of dealing with five different galaxy sub-catalogues and two foregrounds.

To further test the statistical efficiency of the Gibbs sampling procedure we estimated its efficiency in generating independent Markov samples. Generally subsequent samples in a Markov chain are correlated and do not qualify for independent samples of the target posterior distribution. To estimate how many independent samples can be drawn from a Markov chain with a given number of transition steps one has to determine the length over which sequential samples are correlated. This correlation length characterizes the statistical efficiency of generating independent realizations of a parameter $\theta$ as follows:

[TABLE]

where $n$ is the distance in the chain measured in iterations, $\left\langle\theta\right\rangle=1/N\sum_{i}\theta^{i}$ and $\mathrm{Var}\left(\theta\right)=1/N\sum_{i}\left(\theta^{i}-\left\langle\theta\right\rangle\right)^{2}$ and $N$ is the total number of samples in the chain.

As an illustration in Fig. 5 we show the correlation length for the power-spectrum amplitudes of different modes in Fourier space, as indicated in the plot. It can be seen that the typical correlation length for a BOSS like survey analysis is on the order of $\sim 100$ Markov transitions. These results demonstrate the numerical feasibility of complex full Bayesian analyses of present and next generation surveys.

4.2 Inference of template coefficients

As discussed above the proposed hierarchical Bayesian inference machine aims to account for the uncertainties of systematic effects arising from foreground effects. In particular the ARES framework correctly accounts for all joint and correlated uncertainties of different inference parameters and even across the five different mock galaxy surveys as used in this work. To illustrate this fact in Fig. 6 we present two dimensional marginal posterior distributions for the corresponding foreground template coefficients. As can be seen, different panels indicate various degrees of correlation between different foreground coefficients across the five mock catalogues. We would also like to point out that generally, distributions for foreground template coefficients are highly non-Gaussian, and can have sharp transitions due to the requirement that effective contamination templates are required to have a positive sign.

As an interesting by-product of our sampling procedure we are able to provide effective survey response maps which account for the a priori unknown systematics due to foreground and target contamination. In particular mean and standard deviations of such maps can be estimated by evaluating Equation (4) for every foreground template parameter coefficient in the Markov chain and multiplying it with the estimated completeness map $C_{i,obs}$. The result is demonstrated in Fig. 7. It can be seen that most corrections appear at the boundary of the mask. These are the regions most affected by foreground stars in our test scenario. Note, that corresponding standard deviations, as shown in the right panel of Fig. 7, also show increased uncertainty in these regions. This demonstrates that the algorithm accounts for larger uncertainty for more unreliable portions of the data, and optimally extracts cosmological information from observations, as discussed in the following.

4.3 Inferred three dimensional density fields

Although this work focusses primarily on the inference of cosmological power-spectra it is instructive to also look at inferred three dimensional density fields. In particular we would like to highlight the impact of foreground contaminations on the inference of density fields from galaxy redshift surveys. To do so we compare two ARES analyses of the generated mock galaxy catalogue, with and without foreground treatment. From these two Markov chains we can calculate the ensemble mean density and corresponding variance fields. Results are presented in Fig. 8. As can be seen, the analysis without a detailed treatment of foreground contaminations shows residual large scale features and erroneous power particularly close to the survey boundaries. These regions are affected the most by stellar contamination as indicated by the corresponding foreground map shown in Fig. 2.

In contrast our ARES run with detailed Bayesian foreground treatment shows a homogeneous density distribution throughout the entire observed domain. Also note that variance maps of the corrected and uncorrected ARES run look very similar. Apparently the erroneous features in the uncorrected ARES analysis do not affect the corresponding variance map. Since for the uncorrected run, the data model does not account for the systematic uncertainties associated to foreground contaminations, the reconstructed erroneous large scale power in the field will be fully attributed to the inferred large scale structure. From Fig. 8 it is visually evident that our data model including the treatment of foreground contaminations is much more robust against such misinterpretations. In the following we will have a closer look at inferred cosmological power-spectra.

4.4 Inferred cosmological power-spectra

One of the most important features of the ARES algorithm is its ability to jointly infer three dimensional density fields, corresponding cosmological power-spectra, galaxy biases, noise levels and coefficients of several foreground templates including a detailed treatment of all joint and correlated uncertainties. Since the ARES framework yields proper Markov chains we are able to correctly marginalize over all joint uncertainties, when focussing on the analysis of specific target quantities such as the cosmological power-spectrum. Specifically in our previous work Jasche & Wandelt (2013) we have already demonstrated that the ARES algorithm reveals and correctly treats the anti-correlation between bias amplitudes and power spectrum, a 20% effect across large ranges in Fourier space. Here we also take into account the unknown coefficients of several foreground templates.

To study the impact of foreground contamination on the analyses of cosmological power-spectra in deep galaxy surveys we compare inference results obtained from our two ARES runs with and without corresponding corrections. The results for inferred power-spectra are presented in Fig. 9 where we show the univariate marginal posterior distribution for power-spectrum amplitudes at different modes in Fourier space. For the Markov chain without foreground treatment one can clearly observe excessive power at the largest scales. This observation corresponds to the excessive large scale power observed in corresponding inferred three dimensional density fields, as discussed above. In addition one can observe a slight bias with respect to the true underlying power-spectrum from which mock observations were generated. This can easily be understood by inspecting the data model described in Equation (5), where one can see that there is a certain potential for a degeneracy between foreground coefficients and galaxy biases. If foreground effects are not treated correctly, then some of the foreground contributions will erroneously be compensated by sampled galaxy bias amplitudes, introducing the offset between true and recovered power-spectra shown in Fig. 9. In contrast inferred power-spectra for the run with foreground treatment are unbiased with respect to the true underlying power-spectrum over the full domain of Fourier modes considered in this work. In particular the shape of the recovered power-spectrum at the largest scales is in excellent agreement with the true fiducial model.

We further studied the impact of different foreground coefficients on power-spectrum amplitudes at different scales in Fourier space. In particular we calculated the cross correlation matrix between the foreground coefficients of the five mock catalogues and power-spectrum amplitudes from the posterior samples of the corresponding Markov chain. Results are presented in Fig. 10. It can be seen that correlations and anti-correlations can amount of up to ten percent across all modes in Fourier space.

Additionally we tested whether the sampler correctly accounts for the combined effects of foreground contaminations, galaxy biases and unknown noise amplitudes by estimating the co-variance matrix of inferred power-spectra from the ARES runs. As can be seen in Fig. 11 the co-variance matrix for both runs exhibit strong diagonal shape indicating that the algorithm correctly accounted for the otherwise erroneous mode coupling introduced by survey geometry and foreground effects. Residual off-diagonal contributions amount to less than ten percent.

These results clearly demonstrate the feasibility of dealing with strong but unknown foreground contaminations when inferring cosmological power-spectra from deep galaxy observations.

5 Summary and Conclusion

Major challenges for the analysis of next generation deep galaxy redshift surveys arise from the requirement to account for an increasing amount of systematic and stochastic uncertainties. In particular foreground effects and target contaminations due to e.g. stars and dust, can greatly affect the observation of galaxies. If not accounted for properly, these effects can yield erroneous modulations of galaxy number counts across the sky, which hinders the immediate inference of power-spectra and three dimensional density fields from such galaxy samples.

To address this issue, in this work we have described a fully Bayesian treatment of unknown foreground contamination in the process of inferring cosmological power-spectra from deep galaxy surveys. In particularly we build upon the previously presented Bayesian inference framework ARES (for details see Jasche et al. 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2015). The ARES algorithm aims to jointly infer three dimensional density fields, corresponding cosmological power-spectrum, luminosity dependent galaxy biases and unknown noise levels. Being a full Bayesian inference engine ARES further correctly provides joint and correlated uncertainties for all target quantities by performing an efficient Markov chain Monte Carlo sampling within a block sampling scheme as indicated in Fig. 1.

In this work we extend this hierarchical Bayesian framework by also including an additional sampling procedure to account for foreground and target contamination effects. As discussed in Section 2.2 such contaminating effects particularly affect the estimation of the spectroscopic completeness of a given galaxy survey. Naive estimation of the probability of acquisition of galaxy spectra at given positions in the sky, by calculating ratios of observed galaxy spectra and all observed photometric targets, ignores the possibility that photometric targets are contaminated by foreground stars or dust extinction. Such effects are likely to artificially increase or decrease the estimated number of galaxy targets in observations. In consequence these effects introduce an artificial modulation of observed galaxy number densities across the sky, which in turn yields erroneous large scale power in inferred cosmological power-spectra.

As demonstrated in Section 2.2, foreground and target contaminations can be accounted for by describing the mismatch between actual galaxies and observed photometric targets as multiplicative correction factors to estimated spectroscopic completeness maps. These correction factors can be described as combinations of various templates for foreground effects, such as introduced by stars or dust, and corresponding unknown template coefficients. The aim of this work is to jointly infer these template coefficients together with the three dimensional density field, corresponding cosmological power-spectra, galaxy biases and unknown noise levels. This goal can be achieved by simply adding an additional sampling scheme for foreground coefficients to the block sampling framework of the ARES algorithm. Interestingly, as discussed in section 2.3 and Appendix A, we were able to derive a direct sampling scheme by introducing auxiliary random fields and marginalizing over them. The corresponding algorithm to generate random realizations of foreground template coefficients is given in Algorithm 1.

We test the performance of improved ARES algorithms via applications to mock galaxy observations. To further evaluate the impact of foreground and target contamination effects on the inference of cosmological power-spectra we compare two test scenarios, with and without treatment of foreground effects. Corresponding artificial galaxy observations were self-consistently generated according to the data model described in Section 3. In particular, these artificial observations seek to emulate realistic features of the SDSS DR7 main, the LOW-Z and the CMASS galaxy sample. Effectively this results in five artificial galaxy surveys that are jointly handled by the ARES framework, while self-consistently accounting for their respective systematic and stochastic uncertainties, including survey geometries, selection effects, galaxy biasing and foregrounds.

This artificial data was then used to test the statistical performance of our sampling framework. In particular we tested the burn-in behaviour of the algorithm by starting the Markov chain from an over-dispersed state. As described in Section 4.1, we start the chain with an initial power-spectrum scaled by a factor 0.1. The following burn-in behaviour then manifested itself by a coherent drift of sequential power-spectrum samples towards preferred regions in parameter spaces. We estimated the burn-in phase to be completed after $\sim 2000$ sampling transitions. The statistical efficiency of the sampler was estimated by measuring the correlation length between subsequent posterior power-spectrum samples. As demonstrated in section 4.1 the sampler exhibits a correlation length of a few hundred samples. This leaves us with a numerically efficient sampling framework to explore cosmological power-spectra in deep galaxy surveys. It should be remarked that there exist various possibilities to further improve statistical efficiencies but details are left to future publications.

Results for inferring foreground and target contamination template coefficients are described in section 4.2. In particular we have used two realistic foreground templates describing foreground stars in the galaxy and dust extinction. Furthermore, our implementation is general enough to account for an arbitrary amount of foreground templates. To demonstrate the feasibility of inferring these parameters jointly with the cosmological power-spectrum, three dimensional density fields, galaxy biases and noise levels, we presented two dimensional marginalized distributions for template coefficients. These results show that different contamination contributions can be recovered within $\leq$ 2.3 sigma of their input values. Given inferred foreground coefficients it is further possible to reconstruct an effective completeness mask and corresponding uncertainties. As expected, in the example tested in this work, uncertainties in recovered completeness masks are largest in regions where stellar contaminations of the target distribution are also large.

In Section 4.3, we have studied the impact of foreground and target contaminations on the inference of three dimensional density fields from galaxy surveys. To do that we have contrasted a run with and without foreground treatment. Ignoring foreground effects when inferring density fields yields excessive large scale power particularly in regions most affected by such contaminations. In contrast detailed Bayesian treatment of foreground systematics yields inferred density fields showing a homogeneous distribution of power across the inference domain. It must be remarked that if foreground contaminations are not explicitly modelled within the data model, then their effects will be attributed to a real signal.

This result is also in agreement with inferred cosmological posterior power-spectra, as presented in section 4.4. The inferred power-spectrum of the run without foreground treatment agrees with the visual impression obtained from corresponding three dimensional density fields. In particular it reflects the observed excess in large scale power. Ignoring foreground effects may further lead to an overall bias across the entire Fourier domain with respect to the true underlying power-spectrum. In contrast detailed Bayesian foreground treatment yields inferred power-spectra in agreement with respect to the underlying fiducial truth over the entire range of Fourier modes, as considered in this work.

We have further tested the impact of different foreground effects on the inference of power-spectrum amplitudes by estimating their correlations with foreground template coefficients throughout the full Fourier range. These results show that correlations and anti-correlations can amount up to ten percent throughout large ranges in Fourier space.

The ARES algorithm also accounts for artificial mode coupling between power-spectrum amplitudes as introduced by survey geometries and completeness masks. To demonstrate that fact we have estimated the correlation matrix of power-spectrum amplitudes from our Markov runs. These tests show that residual artificial mode coupling is typically much less than ten percent. These results indicate the validity of the algorithm in scenarios with heavily masked data. We are currently using our method to process the actual BOSS data, which will presented in our companion paper (Lavaux & Jasche 2017).

As demonstrated in this work detailed treatment of foreground and target contamination is essential to recover unbiased estimates of three dimensional density fields and corresponding cosmological power-spectra from present and next generation surveys. The proposed ARES algorithm provides a joint and statistically rigorous Bayesian inference framework to achieve this goal and prevent misinterpretation of observations.

Acknowledgements.

This research was supported by the DFG cluster of excellence ”Origin and Structure of the Universe” (www.universe-cluster.de). This work made in the ILP LABEX (under reference ANR-10-LABX-63) was supported by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. The Parkes telescope is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. We acknowledge financial support from ”Programme National de Cosmologie and Galaxies” (PNCG) of CNRS/INSU, France. This work was supported in part by the ANR grant ANR-16-CE23-0002.

Appendix A Sampling foreground coefficients

In this appendix, we derive the sampling procedure for a single foreground template coefficient $\alpha_{k}$. As described in Section 2.3, the full joint distribution of all template coefficients may be sampled by performing a sequential iterative block sampling procedure. In the absence of any a priori information on the foreground coefficient $\alpha_{k}$ we follow a maximally conservative approach by assuming a uniform prior distribution. Using Equation (8)), the logarithm of the conditional posterior distribution for a foreground coefficient in a single galaxy catalogue can then be written as:

[TABLE]

We can simplify notation by noting that different foreground templates contribute multiplicatively an effective survey response operator. We can therefore collapse all multiplicative foreground contributions, except the one currently under consideration, into the effective survey response operator given as:

[TABLE]

where $k$ labels the currently considered foreground template and we have used Equation (4) to factorize foreground contributions. We further introduce the vector:

[TABLE]

Given this notation the conditional posterior distribution for $\alpha_{k}$ simplifies to:

[TABLE]

The expression can be further compressed by introducing the following indexed quantities:

[TABLE]

Consequently the conditional posterior distribution can be expressed as:

[TABLE]

In order for Equation (19) to represent a proper probability distribution the following positivity requirement for the variances needs to hold:

[TABLE]

Similarly it is required that:

[TABLE]

which states that the survey response operator should be positive definite. To ensure these requirements we split the effective survey response operator $\tilde{R}_{i}$ as follows:

[TABLE]

By also requiring $\tilde{R}^{\prime}_{i,k}>0$ we yield the following requirement for the scalar quantity $\omega$:

[TABLE]

We therefore choose $\omega$ to be:

[TABLE]

Given these definitions one can express the conditional posterior distribution as:

[TABLE]

Note that this distribution can be described conveniently as the marginalization over a set of auxiliary fields $\{t_{i}\}$:

[TABLE]

This approach therefore follows a similar line of reasoning as discussed in our previous work when presenting a messenger field Gibbs sampler (Jasche & Lavaux 2015). Here we propose to jointly sample the template coefficient parameter $\alpha_{k}$ with the auxiliary messenger field $t_{i}$ via a two step block sampling procedure. First we generate realizations of the messenger field $t_{i}$ conditional on the current value of $\alpha_{k}$ and then we draw a new value of $\alpha_{k}$ conditional on the given realization of $t_{i}$. We loop over this small block ten times to ensure a minimal mixing of the variables. To simplify notation we introduce the following change of variable:

[TABLE]

This yields the joint distribution:

[TABLE]

As can be easily confirmed, sampling the messenger field $t_{i}$ amounts to simply generating normal random variates with following means and variances:

[TABLE]

and

[TABLE]

To generate realizations of the $\xi$ values, we introduce the following quantities:

[TABLE]

and

[TABLE]

With this definition the conditional distribution to sample $\xi$ turns into a generalized inverse Gaussian (GIG) distribution given as:

[TABLE]

where $N_{v}$ is the number of observed grid elements. The GIG distribution can be conveniently sampled with standard approaches as described in the literature (see e.g. Dagpunar 1988). Finally to obtain a sample for the foreground coefficient $\alpha_{k}$ we invert the transformation in Equation (28):

[TABLE]

The respective realizations of the $t_{i}$ field are not longer required and are immediately discarded, which amounts to a marginalization over the $t_{i}$ values. An efficient sampling algorithm that avoids storing the full $t_{i}$ vector is proposed in Algorithm 1.