Bayesian inference from photometric redshift surveys

Abstract

We show how to enhance the redshift accuracy of surveys consisting of tracers with highly uncertain positions along the line of sight. This increased redshift precision is achieved by imposing an isotropy and two-point correlation prior in a Bayesian analysis and is independent of the process that estimates the photometric redshift. In particular, our method can deal with arbitrary forms of redshift uncertainties for each galaxy. As a byproduct, the method also infers the three-dimensional density field, essentially super-resolving high-density regions in redshift space. Our method fully takes into account the survey mask and selection function. It uses a simplified Poissonian picture of galaxy formation, relating preferred locations of galaxies to regions of higher density in the matter field. The method quantifies the remaining uncertainties in the three-dimensional density field and the true radial locations of galaxies by generating samples that are constrained by the survey data. The exploration of this high-dimensional, non-Gaussian joint posterior is made feasible using multiple-block Metropolis–Hastings sampling. We demonstrate the performance of our implementation on a simulation containing 2 × 10⁷ galaxies. We further demonstrate the robustness of our method to prior misspecification by application to mock observations built from a large-scale structure simulation. In this test, initial Gaussian redshift uncertainties with δz ∼ 0.03 can yield final redshift uncertainties of δz_f ∼ 0.003 in high-density regions. These results bear out the promise of Bayesian analysis for upcoming photometric large-scale structure surveys with tens of millions of galaxies.

methods: numerical, methods: statistical, cosmology: observations, large-scale structure of Universe

Introduction

Understanding the origin and evolution of the Universe to high precision requires increasing amounts of data. In particular, understanding the nature of the recently observed accelerated expansion, measuring the Baryon acoustic scale in the matter distribution or establishing new standard rulers require probes of the three-dimensional (3D) galaxy distribution over large volumes, out to high redshifts. A particularly difficult problem in galaxy observations arises from the requirement of accurate redshift measurements as input to high-precision cosmology. Spectroscopy is the most accurate method for measuring the radial distances of galaxies along the line of sight. However, it is time-intensive and needs high photon fluxes. Therefore, spectroscopic redshifts will be unavailable for the vast majority of the (tens of) millions of galaxies to be observed by future galaxy surveys (Benjamin et al. ). These requirements have paved the way for photometric redshift (photo-z) measurements as valuable astronomical tools in future observations of the galaxy distribution. Photometric galaxy observations can not only measure redshifts for fainter objects than could be yielded with spectroscopy but also provide higher numbers of objects with redshift measurements per unit telescope time (Hildebrandt et al. ). Consequently, using photometry rather than spectroscopy permits exploring the galaxy distribution to higher redshifts at the expense of increased uncertainty in the determination of 3D galaxy positions (Baum ; Koo ; Connolly et al. ; Fernández-Soto et al. ; Blake & Bridle ).

Since major on-going and future galaxy surveys, such as Dark Energy Survey (DES), Physics of the Accelerating Universe (PAU), Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) and the space-based EUCLID mission, will be photometric, redshift uncertainties will play an increasingly important role in future observations of the large-scale structure (see e.g. Benítez et al. ; Coupon et al. ; Refregier et al. ; Frieman ). Nevertheless, precision cosmology requires accurate distance estimates. The impact of redshift uncertainties on the inference of cosmological signals such as the baryon acoustic oscillations, cluster counts or weak lensing, has been extensively discussed in the literature (see e.g. Bernstein & Jain ; Huterer et al. , ; Ishak ; Ma, Hu & Huterer ; Zhan & Knox ; Zhan ).

In this work we present a new Bayesian approach for redshift inference from galaxy redshift surveys with uncertain radial locations of galaxies. We emphasize that our method is independent and complementary to the methods inferring photo-zs from colour or morphology information, and is therefore intended to be used in conjunction with those. In order to increase the accuracy of galaxy positions our method exploits two physically motivated assumptions. The first assumption relies on the cosmological principle which states that the Universe should be homogeneous and isotropic. For this reason, the covariance matrix of the matter distribution in Fourier space is of diagonal shape, at least at the largest scales, with the diagonal being the cosmological power-spectrum. As a consequence, the sampler searches for galaxy redshift configurations that are compatible with an isotropic matter distribution. Secondly, we exploit a simplified physical picture of galaxy formation, which states that galaxies predominantly form in regions of higher density.

In practice, neither the 3D matter field nor the exact galaxy positions are known a priori. We must infer the 3D density field together with the locations of galaxies. In so doing, we account for the mutual dependences and correlated uncertainties of the inferred density field and galaxy positions, and show that the joint treatment can enhance their mutual inference.

In order to perform a full joint Bayesian analysis of the 3D density field and the redshifts of galaxies we implement a two-step Markov sampling process as depicted in Fig. ¹, which permits to explore the high-dimensional non-Gaussian joint posterior distribution in a numerically efficient way. In particular, within the framework of a multiple block Metropolis–Hastings sampler, we are able to break up the high-dimensional problem in a set of lower dimensional problems which can then be treated iteratively (Hastings ). This sampling approach permits us to produce data-constrained samples of the 3D density field and galaxy redshifts for photo-z surveys with tens of millions of galaxies in a numerically efficient way.

Figure 1.

Flow-chart depicting the two-step iterative block sampling procedure. First a set of galaxy redshift reconstructions is sampled conditional on the last density field. In a second step, a density field sample is generated conditional on the new galaxy redshifts, assuming them to be related to the density field through log-normal Poisson statistics.

In this work we are particularly interested in the final redshift uncertainties, which are output as sampled probability densities of galaxy redshifts constrained by the observations and incorporating the isotropy prior. These results improve estimates for the radial locations of galaxies and quantify the remaining uncertainties.

The paper is structured as follows. In Section 2 we give a brief overview and introduction to the method presented in this work. A brief overview of photo-z estimators that generate the input needed for our work is provided in Appendix . The density field sampling step and the corresponding log-normal Poissonian model for the density field is briefly discussed in Section 3. Since this sampling step has already been solved within the Hamiltonian Density Estimation and Sampling (HADES) algorithm, a more detailed description can be found in Jasche & Kitaura () and Jasche et al. (). The redshift sampling process and the derivation of the corresponding conditional redshift posterior distribution are discussed in Section 4. In Section 5 we describe the generation of the mock data used in our tests.

The performance of our multi-dimensional Metropolis–Hastings sampler in a realistic setting is studied in Section 6. Here we study the overall burn-in and convergence behaviour of the Markov sampler. In Section 7 we discuss the results for the inferred redshifts and 3D density fields. Finally, in Section 8, we discuss our results and give an outlook on future improvements and applications of our work.

Method

The method presented in this work aims at the joint inference of radial locations of galaxies as well as 3D density fields from observations with uncertain redshift measurements. In particular, we are interested in exploring the joint probability of the real-space density field {s_i} and the true redshift positions of galaxies {z_p} conditional on the observed redshifts and the positions in the sky {θ_n}, which is a set containing the declination and right ascension angles for all galaxies. The exploration of the joint probability distribution is numerically challenging partly due to high dimensionality but also because the problem requires non-linear and non-Gaussian inference methods. In situations like this, one usually has to rely on efficient implementations of a Metropolis–Hastings sampling scheme. However, it may be difficult to construct a single block Markov algorithm that converges rapidly to the target distribution.

Fortunately, it is possible to break up the parameter space into smaller blocks and to construct transition kernels for each individual block (Hastings ). Consequently, this multiple block sampling framework permits to build complex transition kernels for the joint problem. This approach greatly simplifies the search for suitable candidate-generating densities and enables us to use the previously developed non-linear log-normal Poissonian density sampler, presented in Jasche & Kitaura () and Jasche et al. (), in conjunction with a new redshift sampler described in this work. This block sampling approach was also at the heart of the Gibbs sampler for cosmic microwave background (CMB) power-spectrum inference originally proposed in Wandelt, Larson & Lakshminarayanan () and Jewell, Levin & Anderson (), which is now widely used.

We can explore the joint posterior distribution by iteratively sampling from the two conditional distributions

where in the first step we draw a random realization of radial galaxy positions {z_p} and in a second step we generate a 3D density field sample {s_i}. The iteration of these two random processes will provide us with samples of the joint posterior distribution. A schematic description of the sampling procedure is depicted in Fig. ¹.

The log-normal Poissonian model

The log-normal Poissonian posterior distribution is formulated in terms of galaxy number counts, rather than in terms of galaxy redshifts and positions in the sky. As demonstrated in Appendix , the density posterior distribution given in equation can be written as a function of galaxy number counts

where the number counts N_i({z_p}, {θ_p}) are calculated according to equation . Following Jasche & Kitaura (), the log-normal Poissonian posterior distribution can then be expressed as

where Q is the covariance matrix of the log-normal distribution, μ_i describes a constant mean field, R_k is a linear response operator, incorporating survey geometries and selection effects, is the mean number of galaxies in the volume and B(x)_k is a non-linear, non-local, bias operator at position x_k. In this fashion, the log-normal Poissonian distribution accounts for non-linearities in the density field, the local noise structure of the galaxy distribution as well as systematic uncertainties due to survey geometry, selection effects or galaxy biasing. Besides being physically motivated (see e.g. Coles & Jones ; Kayo, Taruya & Suto ), the log-normal distribution is also the maximum entropy prior for the density field on a logarithmic scale, when knowledge about the mean and covariance is available. It has further been shown that the reconstructed density field, employing a log-normal prior, agrees very well with the true density field, even at very high densities (see Kitaura, Jasche & Metcalf ). For these reasons, the log-normal distribution is the ideal least informative prior for the fully evolved density field in non-linear regimes.

The log-normal posterior distribution incorporates essentially two physically motivated assumptions that are important for the method described in this work. Due to the cosmological principle we will assume the matter distribution in real space to be statistically isotropic. This is incorporated in the posterior distribution in equation by noting that the Fourier transform of the covariance matrix Q for the logarithmic density field is assumed to have a diagonal form with the diagonal being the power-spectrum that depends on the |k| only. This assumption requires the individual density field samples to be statistically isotropic. In addition, the Poisson model can be considered as a simplified model for galaxy formation, which relates the number of galaxies found at a given position to the underlying density field, and hence singles out preferred locations for galaxies conditional on the underlying density field. Thus, in the absence of additional constraining information, these physical requirements enhance the inferred redshifts from observations with uncertain radial locations of galaxies.

For a detailed description on the numerical implementation of the large-scale structure sampler we refer to Jasche & Kitaura (). Unlike conventional Metropolis–Hastings algorithms, which move through the parameter space by a random walk, and therefore require prohibitive amounts of steps to explore high-dimensional spaces, the Hybrid Monte Carlo (HMC) sampler suppresses random walk behaviour by introducing persistent motions of the Markov chain through the parameter space (Duane et al. ; Neal , ). In this fashion, the Hamiltonian sampler maintains a reasonable efficiency even for high-dimensional problems (Hanson ).

The redshift posterior distribution

In the following we will focus on the derivation and implementation of the second sampling step in Fig. ¹, namely of generating random realizations of true galaxy redshifts {z_p} conditional on observations and a realization of the 3D density field. In particular, we are searching for a numerically efficient implementation of the random process

This problem is numerically challenging due to the high dimensionality, i.e. the number of galaxies in a survey. Additionally, the non-linearity of the problem generally prevents drawing samples directly from the redshift posterior distribution given in equation . Hence, in order to explore the space of possible galaxy redshift configurations we have to rely on a Metropolis–Hastings sampler. A naive implementation of the redshift sampling process might suggest drawing a realization for the whole set {z_p} at once. Such a scheme will lead to a high rejection rate, rendering such a method numerically unfeasible. Suppose, for example, that the sampler proposes a redshift configuration {z_p} which is acceptable except for one galaxy redshift z_n. Although most of the galaxy redshifts could be accepted in this scenario the complete set {z_p} would be rejected and the sampler stays at its current position in parameter space. For this reason, such a sampling approach would require to construct a very good proposal distribution for in order to yield efficient sampling rates. Generally, it is not trivial and may even be impossible to construct such a proposal distribution.

If it were possible to find conditional redshift posterior distributions for each individual galaxy redshift z_n, one could build a multiple block Metropolis–Hastings sampler for the set of galaxy redshifts {z_p}.

Fortunately, this is possible in our simple, Poissonian galaxy formation model. As demonstrated in Appendix it is possible to find a conditional redshift posterior distribution such that a multiple block Metropolis–Hastings sampler can be constructed.

Galaxy redshifts can therefore be updated by sampling using the following random process

where for each transition we use a Metropolis–Hastings update step and is the set of galaxy redshifts containing all redshifts except the redshift z_n of the nth galaxy under consideration. The superscripts i at the {u_p}ⁱ indicate that the set of galaxy redshifts {u_p} have to be updated after each sampling step, since the next sampling step depends on the new position of the previously sampled galaxy. The conditional redshift posteriors describes the probability for the redshift z_n of the nth galaxy redshift conditional on the redshifts of all other galaxies, the density field, the observed redshifts and the positions in the sky.

At this point an important further simplification occurs. Given the density field s, the redshift of any one galaxy z_n is conditionally independent of the redshifts of all the other galaxies, so

which is also shown in Appendix B. This is particularly interesting for the numerical implementation of the random process given in equation , since it permits parallel, instead of sequential, sampling.

The computational advantage of sampling each galaxy redshift individually can be immediately realized. Even if some individual redshift proposals are rejected, other redshift samples will be accepted leading to an overall larger acceptance probability for the entire set {z_p}, compared to the scenario where {z_p} would be sampled as a whole. Additionally, the multi-dimensional sampling approach can be reduced to a one-dimensional sampling problem for each sampling step described in equation , as the redshift is explored along each line of sight.

The detailed form of the conditional redshift posterior is (Appendix B)

where r(z) is the radial comoving distance as a function of redshift z and W(x) is the voxel kernel function, as defined in Appendix . Since the Metropolis–Hastings algorithm does not require a probability distribution to be normalized in order to generate samples, we omitted all normalization constants in equation .

The first factor in equation is just the Jacobian of the coordinate transformation from spherical redshift to Cartesian comoving coordinates, which is required since the density field is defined in comoving Cartesian coordinates. The second factor is essentially the intensity function for the spatial Poisson process, as described in Section 3. It reflects the naive intuition that the probability of finding a galaxy at a given spatial location should be related to the underlying density field. Also note that the conditional redshift posterior distribution accounts for non-linear and non-local galaxy biasing B(s)_i, as well as for the systematics of the survey under consideration via the survey response operator R_i, which incorporates survey geometry and selection functions. The last factor in is the photo-z likelihood encoding the information about the redshift for the nth galaxy determined by a photo-z estimator. Since galaxy redshifts are determined separately for individual galaxies it is reasonable to assume that these likelihoods are independent of the observed redshifts of all other galaxies, meaning the redshift likelihood for an individual galaxy does not depend on the redshift observations of other galaxies.

We emphasize that so far we made no specific assumptions regarding the photo-z likelihood. The aim of this work is to find a physically meaningful way of increasing the accuracy of redshift inference once all other, colour or spatial information has been used. This is achieved by the two first factors in equation which could be considered as a redshift prior once the underlying density field is known. Thus, the joint sampling approach for the density field and the galaxy redshifts described in Section 2 provides us with samples of the joint distribution. Discarding the density fields then amounts to a Monte Carlo marginalization over the underlying 3D density fields and provides us with the redshift posteriors for all galaxies conditional only on the data and the isotropy prior.

It is important to remark that since no specific assumptions on the redshift likelihood are required, and since we are using a sampling approach, our method will in principle work with any specified redshift likelihood, or input photo-z probability distribution function (PDF). In particular, the Metropolis–Hastings algorithm employed to explore the redshift posterior distribution permits us to sample from arbitrary redshift likelihoods without any specific restrictions. Our method is therefore intended to be used as a complement to any photo-z estimator that uses galaxy-specific information such as a photometry or morphology to provide a photo-z PDF.

For this reason, our method is highly flexible and therefore broadly applicable to various different scenarios and galaxy observations.

Generating Mock Observations

In the previous sections we presented the derivation and the implementation of our method. Here we will describe the generation of artificial mock data sets, which will be used later to test our method.

Generating galaxy redshifts

The generation of an artificial mock galaxy survey will closely follow the description in Jasche & Kitaura (). Since galaxies trace the 3D matter distribution and in this work we are interested in recovering the true redshift locations of galaxies from redshift surveys with uncertain radial locations of galaxies, we prefer using a density field simulated from an N-body simulation rather than a log-normal realization. Such a density field permits us to test our method in a more realistic scenario, since the mock density field contains structures, such as voids, filaments and clusters. Furthermore, by using fully non-linear density fields of N-body simulations to construct our test cases, we will be able to demonstrate the applicability of the least informative log-normal prior in such a setting.

For this reason, we build our mock galaxy observation on a matter field realization calculated by the tree-PM code GADGET-2 (Springel ), starting from Gaussian initial conditions. The details of this simulation are described in Pillepich, Porciani & Hahn ().

We estimated the simulated matter density field on a 256³ Cartesian equidistant grid with side length of 1200 Mpc, using the Cloud in Cell (CIC) kernel and the super-sampling procedure described in Jasche & Kitaura () to account for aliasing effects. Note that we do not take peculiar velocities into account when mapping into redshift space, since any such distortion will be negligible for photometric survey data.

Given this density field we can then generate Poisson realizations, taking into account several survey properties such as noise, survey geometry and selection effects. In Fig. ² we show slices through the true underlying matter field realization and the simulated redshift survey.

Figure 2.

Slices through the mock matter field realization (left-hand panel), the mock galaxy survey generated from the mock matter field without distorted redshifts (middle panel) and the redshift distorted mock galaxy survey (right-hand panel). It is clearly visible that all small-scale structures are smeared out along the line of sight in a survey with photo-z errors of ∽0.03.

The survey properties are described by the galaxy selection function F_j and the observation Mask M_j where the product

yields the linear response operator R_j of the galaxy survey at the jth volume element.

The toy selection function, used in our tests, is given by

where r_j is the comoving distance from the observer to the centre of the jth voxel. For our simulation we chose parameters β = 0.25 and r_th = 100 Mpc. The observation Mask M_j is the same as described in Jasche & Kitaura (). It was designed to approximate the most prominent features of the Sloan Digital Sky Survey geometry (Abazajian et al. ). As a result of the Poisson sampling of the density field we obtain galaxy number counts N_j at each grid-point in the 3D domain.

Next, in order to assign redshifts to the galaxies, all the N_j galaxies inside the jth volume element are uniformly distributed over the volume element to offset them from the centre position. This process does not destroy the Poissonian process by which the galaxy sample was realized, since it does not change the number of galaxies inside the volume element. Given an assumed cosmology we can then calculate the redshift z_p and the two angular positions θ_n for each galaxy. In particular, we assumed a standard Λcold dark matter cosmology corresponding to the following set of cosmological parameters ().

At the end of this procedure we obtained a set of real-space galaxy positions {θ_n, z_p} for N_gal = 21 340 276 galaxies.

Introducing redshift uncertainties

In the previous section we described how to obtain a set of unperturbed galaxy redshifts given a specific density field. Here we will describe how to generate the observed photo-z uncertainties.

Note that until now we have not specified a specific redshift likelihood. The method presented in this work is general, and does not need to make any specific assumptions about the functional shape of the redshift likelihood. In particular, the Metropolis–Hastings sampling approach permits us to use any desired redshift likelihood, whether given in analytic or numerical form.

In this work we will assume that the uncertainty in the observed redshifts is due to the observational strategy, i.e. due to the instrument. The redshift likelihood is therefore independent of the underlying density field

Although it is possible to use an arbitrary redshift likelihood, here, as an example, we will employ a truncated Gaussian distribution given as

where σ_z is the photo-z dispersion. This simple distribution permits us to test our method and to estimate its general performance. Given the true redshift z_n generated via the procedure described in Section 5.1 and by specifying a desired redshift dispersion σ_z, the observed redshift can be generated by drawing random samples from the truncated Gaussian distribution, given in equation , via rejection sampling. We distort the redshift of each individual galaxy following this recipe.

For the test case considered in this work we treat all galaxies with the same σ_z = 0.03. Since we are able to treat each galaxy individually, there is no problem in specifying a different photo-z likelihood for every galaxy. The possibility of modelling heterogenous data sets with photo-z PDFs that differ from galaxy to galaxy is very interesting since it allows incorporating detailed uncertainty information from photo-z estimators which could take into account colour or morphological information. It would also be possible to merge data sets with different redshift accuracies, including spectroscopic surveys. We will explore these promising ideas in future work.

The process of generating a galaxy distribution with redshift distortions corresponding to the likelihood given in equation is visualized in Fig. ². It can be seen that given a redshift dispersion of σ_z = 0.03 structures along the line of sight are smeared out on a scale of roughly ∼100 Mpc.

In the following we will apply our method to this mock survey.

Testing

In this section we will apply our method to our mock survey in order to estimate its performance. In particular, we focus on the convergence behaviour of the sampler which determines the efficiency of the method in a realistic setting.

Testing convergence and correlations

The Metropolis–Hastings Sampler is designed to have the target distribution, in our case the joint posterior distribution, as its stationary distribution (see e.g. Metropolis et al. ; Hastings ; Neal ). For this reason, the sampling process will provide us with samples from the specified joint posterior distribution after a sufficiently long burn-in phase. However, the theory of Metropolis–Hastings sampling by itself does not provide any burn-in criterion. Therefore, this initial burn-in behaviour has to be tested using numerical experiments.

Burn-in manifests itself as a systematic drift of the sampled parameters towards the true parameters from which the artificial data set was generated. We monitor it by following the evolution of parameters in subsequent samples.

As an example, in Fig. ³ we show slices through galaxy number counts calculated from successive redshift samples. While initially the galaxy distribution exhibits prominent redshift distortions, successive samples iteratively correct the galaxy positions towards their true real-space positions. Visually, already the 100th sample does not seem to be contaminated by redshift distortions. In order to demonstrate that galaxies are really moving towards their true real-space positions during the burn-in phase, we can study the correlation of the sampled redshifts with the true galaxy redshifts for successive samples. The results of this test are presented in Fig. ⁴. Here it can be seen that the correlation between the sampled and the true redshifts tightens throughout successive samples. Fig. ⁴ also clearly demonstrates the ability of our method to enhance significantly the constraints on the true galaxy redshifts from photo-z surveys.

Figure 3.

Slices through three successive galaxy count samples, first sample (left-hand panel), fifth sample (middle panel) and 100th sample (right-hand panel). The first sample shows prominent redshift smearing which quickly decays in subsequent samples.

Figure 4.

Correlation between the true redshift z_p and the redshifts of the ith sample for three successive samples, tenth sample (left-hand panel), 100th sample (middle panel) and 5000th sample (right-hand panel). It can be seen that in successive samples the correlation gets stronger as the galaxies move towards their true locations.

These aforementioned diagnostics particularly focus on monitoring the burn-in behaviour of the redshift samples. In order to probe the evolution of density samples we monitor power-spectra estimated from successive density samples. As the galaxies move towards their true locations these successive power-spectra will drift towards the power-spectrum of the underlying mock matter field realization. These successive power-spectra are depicted in Fig. ⁵. As can be clearly seen, successive power-spectra move towards the true power-spectrum of the underlying mock matter field realization. In summary these tests indicate that the initial burn-in phase requires of the order of 5000 samples.

Figure 5.

Power-spectra P(k) of successive samples during the burn-in phase. It can be seen that the samples move towards the true power-spectrum (black curve) of the matter field realization.

Next we will discuss the general convergence behaviour of the redshift samples in particular the correlation structure. In general, successive samples in a Markov chain are not independent, but are correlated with previous samples. In order to estimate how many independent samples are generated by the Markov chain one has to study this correlation effect. This can be done by following a similar approach as described in Eriksen et al. () or Jasche et al. () by assuming all individual galaxy redshifts in the Markov chain to be independent and estimating their correlation in the chain by the calculation of their autocorrelation function

where n is the distance in the chain measured in iterations. In Fig. ⁶ we show the autocorrelations of redshift samples for several randomly selected galaxies. We can further define the correlation length of the redshift sampler as the distance in the chain n_C beyond which the correlation coefficient C_l(n) has dropped below 0.1. As can be seen in Fig. ⁶ the correlation length is less than 200 samples, which demonstrates the high mixing rate of the sampler.

Figure 6.

Correlation functions C_l(n) as a function of distance in the Markov chain n for 10 randomly selected galaxies. It can be seen that the correlation length between the samples in the chain is of the order of 200 samples.

Despite the high dimensionality of the problem considered here, these tests demonstrate that the joint posterior of galaxy redshifts and 3D density field can be explored in a numerically feasible way by our method.

Results

In the previous section we demonstrated the ability of our method to efficiently explore the joint parameter space of the 3D density field and the true redshifts of galaxies from redshift surveys with uncertain radial location of galaxies.

The algorithm will identify high-density regions in volume elements where the number counts of galaxies are high. In this fashion, the information of the entire galaxy sample is used to infer the 3D density field. In turn, the inference of individual galaxy redshifts uses this reconstructed 3D density field as described in Section 4. Therefore information propagates from the entire galaxy sample to reduce the redshift uncertainties of each individual galaxy.

In successive Markov samples individual galaxies will preferentially move to regions of higher density within the range allowed by the photo-z inference. For this reason, a worry might be that the galaxy may get stuck and not move any further once it finds a density peak. We find that this does not occur since our burn-in and autocorrelation tests show that galaxies quickly find and explore their redshift range.

It is illustrative to follow the trajectory of an individual galaxy along its line of sight throughout successive samples. Generally, in each sample step both, the density field along the line of sight and the galaxy redshift, will change. In Fig. ⁷ we show the trajectory of a random galaxy and the evolution of the corresponding density field along the line of sight for 1000 samples. It can be seen that the galaxy preferentially moves to regions with densities higher than the mean density and avoids underdense regions. Also the galaxy stays longer in regions of higher density, which correspond to regions of higher probability, than in other regions.

Figure 7.

The plot shows the trajectory of a randomly selected galaxy along the line of sight as a function of the sampling step n in the Markov chain (blue dashed line). The background displays the density field along the line of sight in each sampling step. It can be seen that the galaxy jumps preferentially between different regions of higher density in the range allowed by the redshift likelihood from a photo-z estimator.

This is the anticipated behaviour of the Markov sampler, which is supposed to provide more samples of high probability regions and less samples in regions of lower probability in order to construct a probability density.

Next, we want to discuss the final scientific results of the sampling procedure. In particular, we are interested in the accuracy of the inferred galaxy redshifts. Intuitively, the accuracy of the inferred redshifts should depend on the environment at the galaxies’ original location. Galaxies born in high-density regions will be attracted by a high-density peak along their line of sight, while galaxies born in low-density regions have fewer constraints. In the limit of very low density the galaxy will explore the full range allowed by the input photo-z likelihood along the line of sight.

High-density regions are better sampled by galaxies, which further decreases the uncertainty for the inferred density field at these positions, which results in an overall higher joint probability for the density field and galaxy redshifts at these locations. For this reason, redshifts of galaxies located in high-density regions such as clusters or super clusters will generally be recovered more precisely than those of galaxies in regions with lower density.

Redshift samples

To demonstrate this effect and to see which regions can be best recovered from a redshift survey with uncertain radial locations of galaxies, we study the average redshift standard deviation as a function of the true underlying density field δ_true

with being the true density contrast at true location of the pth galaxy, σ_p is the galaxy redshift standard deviation calculated from the samples and K(x) is the binning kernel as defined in equation . Equation yields the redshift standard deviation measured from the Markov chain averaged over density bins of the true underlying mock matter-field realization. This quantity is a measure of how well the true redshift can be recovered from an uncertain galaxy redshift survey.

Fig. ⁸ clearly shows a strong dependence of the average redshift standard deviation on the amplitude of the true density field. Very dense structures can be recovered very well. As expected, Fig. ⁸ demonstrates that high-density regions such as clusters can be much better recovered from uncertain redshift surveys than regions of lower density such as filaments or voids. However, on average the redshift precision has improved for all galaxies compared to the initial uncertainty.

Figure 8.

Mean standard deviation as a function of density contrast δ.

The final redshift posterior distributions are generally highly non-Gaussian. We therefore estimate the conditional probability distribution of redshift departures from the true redshift Δz = z_true − z conditional on the true underlying density field δ_true

which is the probability of finding the departure Δz in a region with density contrast δ_true after our method has been applied. Initially, this distribution coincides with the redshift likelihood, which is independent of the density field and is basically a Gaussian distribution with redshift dispersion σ_z = 0.03. The comparison between the initial distribution and that obtained by our method is shown in Fig. ⁹. This plot demonstrates that although the initial redshift likelihood is independent of the true underlying density field, the final redshift posterior distribution strongly depends on it.

Figure 9.

Probability distribution for the difference Δz between measured and true redshift conditional on the true underlying density field δ_true. The left-hand panel shows the initial Gaussian redshift likelihood, as measured from the initial particle distribution, which is conditionally independent of the true underlying density field. The right-hand panel shows the a posteriori conditional probability density distribution for Δz estimated from the Markov samples. It can be seen that the accuracy of inferred redshifts depends on the true underlying density field. Also note that the resulting probability distribution is non-Gaussian.

Beyond measures of accuracy, it is instructive to inspect the redshift posteriors we derive for individual galaxies. In Fig. ¹⁰ we show some randomly selected redshift posterior distributions. In general the inferred redshift posterior distributions are non-Gaussian and multi-modal. Each mode corresponds to a possible association of the galaxy with peaks in the density distribution. The area under each mode measures the probability for the galaxy to be associated with that density peak. When these modes sit on top of a broader distribution this accounts for the residual probability that the galaxy is not associated with any of these peaks. It is clear that by any reasonable measure of information, such as Shannon's negentropy, these PDFs are far more informative regarding the redshift of the galaxy than the original photo-z likelihoods.

Figure 10.

Randomly selected marginalized redshift posterior distributions estimated from the Markov samples (black curves). The red vertical lines denote the true redshift of this specific galaxy and the red curve shows the distribution from which the observed redshift was drawn. The green vertical lines are the observed redshifts and the blue vertical lines denote the ensemble mean redshift estimated from the samples.

The density field

In this section, we will discuss the inferred 3D density field. As already described in the previous section, the accuracy of the inferred redshifts depends on the amplitude of the true underlying density field. In particular, the higher the underlying density the more accurately we can identify the true location of galaxies. For this reason, we expect a similar behaviour for the recovered 3D density field. For example, we expect to recover high-density regions more accurately than low-density regions. However, the corresponding joint uncertainty of redshift and density estimates can be accurately estimated from the Markov samples. As an example we calculated the ensemble mean and variance for the density field. The results are presented in Fig. ¹¹ where we show three slices through the ensemble mean density and the ensemble variance field from different sides. The results demonstrate that high-density regions are much better localized than lower density or even underdense regions. In particular, the ensemble mean density field shows elongated structures along the line of sight for underdense objects such as voids. This result was anticipated. It reflects the fact that these objects are poorly constrained by the data.

Figure 11.

Three different slices from different sides through ensemble mean density (left-hand panels), ensemble variance (middle panels) and the 3D response operator R_i (right-hand panels). Especially the variance plots demonstrate that the method accounted for the full Poissonian noise structure introduced by the galaxy sample. One can also see the correlation between high-density regions and high variance regions, as expected for Poissonian noise.

To further quantify this visual impression and to estimate the accuracy of the inferred density field we study point-to-point comparisons between the true underlying density field and the recovered ensemble mean density field smoothed on different scales. The density fields are smoothed with a spherical top hat filter in Fourier space with different filter radii k_th. We can then calculate the corresponding correlation coefficient for the different filter thresholds k_th

where the superscript denotes that the field has been filtered with a top hat filter corresponding to k_th. The result of this calculation is presented in Fig. ¹². To further study the dependence of the accuracy of the inferred density field on the true underlying density field, we will also calculate the correlation coefficients only between voxels in the true and the inferred ensemble mean density for which δ_true > δ_th. We can therefore study the correlation between true and recovered density fields as a function of smoothing scale and density amplitude. As demonstrated by Fig. ¹², the accuracy for the inferred structures in the density field strongly depends on the amplitude of the true underlying density field. As expected, high-density objects can generally be much better constrained than low-density or even underdense regions. It is interesting to point out that high-density objects still show a reasonably high correlation even at the shortest scales resolved by the grid (∼ 4.69 Mpc) although the data correspond to a density field smeared out on a scale of about 100 Mpc due to a redshift variance of σ_z = 0.03. These results are therefore in agreement with the results for the estimated redshifts presented in the previous section. In the absence of additional information or data the inference of redshifts and 3D density field strongly depends on the amplitude of the true underlying density field. For further illustration we also plot the correlation coefficient between the true underlying density field and the density contrast δ^data estimated directly from the galaxy sample as

As can be seen in Fig. ¹², our method improves the correlation between the true underlying density field and the estimated ensemble mean field on all scales. Also note that the naive estimator given in equation suffers from observational effects such as survey geometry and selection effects at large scales, while these effects have been accounted for in our method.

Figure 12.

Correlation coefficient r(k_th) as a function of filter scale k_th for different density amplitude thresholds δ_th as indicated in the plot. The left-hand panel corresponds to the correlations between the initial density contrast, estimated from the galaxy counts, and the true density field, while the right-hand panel depicts the correlation between the amplitudes of the true underlying and the ensemble mean density field. It can be seen that the correlation between amplitudes in the true underlying density field and the inferred ensemble mean density field increases on all filter scales with increasing density thresholds δ_true > δ_th. Also note the improved correlation at the large scales due to the correction of observational effects such as survey geometry and selection effects.

Discussion and conclusion

Present and upcoming photometric galaxy surveys provide us with (tens of) millions of galaxies sampling the density field out to high redshifts. These enormous numbers of galaxies can only be obtained at the expense of low accuracy in the measurement of the radial locations of galaxies along the line of sight. It has been demonstrated that these redshift uncertainties will generally influence the inference of cosmological signals such as the baryon acoustic oscillations, cluster counts or weak lensing. For this reason, accurate treatment of redshift uncertainties is required in order to use these surveys to full effect.

In this work, we presented a new, fully Bayesian approach that improves the accuracy of inferred redshifts. Specifically our method performs a joint inference of the 3D density field and the true radial locations of galaxies from surveys with highly uncertain redshift measurements. This joint approach permits us to impose physically motivated constraints on the redshift and density inference, which increase the accuracy of the inferred locations of galaxies, essentially super-resolving high-density regions in redshift space.

Where does the additional information come from? The isotropy and 2D correlation prior allow improving the inference of 3D structure from the highly precise angular information in the radial direction. If the projected density of tracers on the sky shows a peak then it is likely that this corresponds to one or more 3D peaks along the line of sight. Many galaxies with imprecise redshifts populate the line of sight. Galaxies in high-density regions can be associated with structures that are consistent with their photo-z likelihoods and the observed two-dimensional density. Each such structure has a relatively accurate redshift derived from the averaged redshift of the galaxies it contains.

In this work we employed a simplified picture of galaxy formation. The Poisson model of galaxy formation serves to connect the isotropy prior on the underlying density field and the 3D locations of galaxies. As demonstrated in Section 4, in this picture, the probability of finding a galaxy at a given radial position is proportional to the intensity of the spatial Poisson process. As can be seen from equation , ignoring survey systematics and galaxy biasing, the probability of finding a galaxy at a given location is essentially proportional to the density, as would be expected intuitively. Further, equation also demonstrates that our method can accurately account for survey systematics and can in principle be combined with a non-linear and non-local galaxy biasing model.

The joint inference of radial galaxy positions and the 3D density field as well as the quantification of joint uncertainties requires exploring the posterior density conditional on observations. Doing so via Markov chain Monte Carlo methods is a numerically challenging problem due to the high dimensionality, non-Gaussianity and non-linearity of the problem.

We show that the complex sampling procedure can be split up into many simpler sub-sampling steps within the framework of a multiple block Metropolis–Hastings algorithm. In particular, this permits us to split the problem into the two tasks of individually sampling 3D density fields, conditional on the true galaxy locations, and sampling galaxy redshifts, conditional on the density field. This sampling approach is depicted in Fig. ¹. We explore the log-normal Poissonian posterior distribution for the 3D density field based on the galaxy distribution using the HADES algorithm (Jasche & Kitaura ). In this work we add a sampling procedure for the redshift posterior distribution as shown as the first step in equation .

In the Poissonian picture of galaxy formation the locations of individual galaxies are conditionally independent of the locations of all other galaxies once the underlying density field is given. We demonstrate this intuitively plausible fact in Section 4. This fact greatly simplifies the numerical implementation of the redshift sampler, since it allows treating each galaxy individually, leading to a numerically efficient, parallel sampling procedure despite the high dimensionality and non-linearity of the problem.

We extensively tested our method with artificial mock galaxy observations, basing on the results of N-body simulations, to evaluate the overall performance of the algorithm in a realistic scenario with tens of millions of galaxies. As a remark we should mention that the method has been successfully tested previously with a Poisson sample from a correlated log-normal realization. However, artificial observations built from large-scale structure simulations represent a more realistic test scenario. The generation of the artificial galaxy survey from a N-body simulation is described in Section 5. In this test all galaxy positions have been distorted corresponding to a truncated Gaussian with equal redshift dispersion for all galaxies. In particular, we used a redshift dispersion σ_z = 0.03 which amounts to an uncertainty along the line of sight of about ∼100 Mpc. However, it should be noted that since our method treats each galaxy individually, it is generally possible to take into account the photo-z information and uncertainties appropriate for each individual galaxy.

The primary goal of the tests we conducted was to establish the numerical feasibility of our method. In particular, we focused on evaluating the burn-in and convergence behaviour of the multiple block Metropolis–Hastings sampler. We studied the burn-in behaviour of the sampler by monitoring the systematic drift of various inferred quantities towards their true values throughout subsequent sampling steps. As described in Section 6.1, following the evolution of inferred power-spectra, shown in Fig. ⁵, and the flow of galaxies towards their true locations, as demonstrated by Fig. ⁴, indicates a burn-in time of about 5000 samples. Further, we studied the correlation length of individual galaxy redshifts across the Markov chain. In particular, Fig. ⁶ demonstrates that the correlation length is typically less than 200 samples. These tests clearly demonstrate the feasibility of the method proposed in this work.

As a scientific result the method provides improved redshifts for each galaxy and the 3D density field as well as the corresponding joint uncertainties. In Section 7.1 we discuss the results for the inferred galaxy locations. As was anticipated, the accuracy of inferred redshifts depends on the density amplitude of the true underlying density field. Fig. ⁸ demonstrates that on average all galaxy redshifts have been improved with respect to the original redshift dispersion, up to a factor of 10 in high-density regions.

It can clearly be seen, from our test cases, that the accuracy of inferred redshifts for galaxies located in high-density regions is generally larger than for galaxies in regions of lower density. This result is further clarified by Fig. ⁹ which shows the probability distribution of departures from the true redshifts conditional on the underlying density field. Redshifts of galaxies in high-density regions can be particularly well recovered from a survey with uncertain radial galaxy positions. These results correspond to our expectation, and simply reflect the fact that there are usually many galaxies in high-density objects. The joint information of these many observed redshifts permits to accurately infer the position of the high-density object and hence the locations of galaxies present in this object. The lower the density the smaller is the number of galaxies contributing to a feature and hence the number of observed redshifts. As a result it is generally more difficult to constrain the locations of low-density features than the location of high-density features.

This effect is also reflected in the inferred density field (see Section 7.2). Fig. ¹¹ shows the ensemble mean density field and the corresponding ensemble variance estimated from the set of Markov samples. It is visually clear that high-density regions can be much more accurately located than low density or even underdense regions. To quantify this impression, we studied the point-to-point correlation of the estimated ensemble mean with the true underlying density field for different smoothing scales and different density thresholds. The results shown in Fig. ¹² clearly demonstrate that high-density regions can be accurately recovered for all smoothing scales, while the correlation between ensemble mean and true density field worsens with decreasing density. For this reason, the results of inferred density fields and galaxy redshifts are consistent and strongly depend on the density amplitudes of the true underlying density field. In addition, by employing mock observations generated from large-scale structure simulations, which do not follow the assumed prior distribution, the method proves to be robust to the choice of a correlated log-normal prior, as anticipated from maximum entropy arguments. However, in spite of using a least informative prior the method proves to achieve significant improvement in the recovery of galaxy positions and the density field.

We emphasize that the scientific output of our method is not a single estimate but a sampled representation of the joint posterior of 3D density field and galaxy positions. Any desired statistical summary as well as corresponding uncertainties can be calculated from this set of samples.

For example, by marginalizing over the density field we are able to provide the marginalized redshift posterior distributions, as depicted in Fig. ¹⁰, which permit easy propagation of uncertainties through subsequent analysis stages.

Our method does not rely on specific assumptions on the provenance of the redshift likelihood. It is therefore applicable generally to surveys of isotropic distributions of tracers with uncertain radial positions.

In summary, the presented method is a flexible and numerically efficient addition to the analysis toolbox for large-scale galaxy surveys and has the potential to add substantial value to their scientific output.

Acknowledgments

We thank Francisco S. Kitaura, Torsten A. Enßlin and Simon D. White for useful discussions and encouraging us to pursue the project described in this work. Further, we thank Cristiano Porciani for providing us with the simulated density field and Nina Roth for providing us with required reading routines and information on how to handle the simulation data. Particular thanks go to Rainer Moll and Björn M. Schäfer for useful discussions and support with many valuable numerical gadgets. This work has been supported by the Deutsche Forschungsgemeinschaft within the Priority Programme 1177 under the project PO 1454/1-1, NSF grant AST 07-08849 and by NFS grant NSF AST 09-08693 ARRA.

Methods for Obtaining Redshifts from Photometry

The approach we present in this paper uses the output of photo-z estimators as its input. We give a short summary of these methods and references to relevant work below.

A variety of different methods to infer accurate redshifts from photometric galaxy observations have been proposed. These techniques can be roughly categorized as kernel regression methods, artificial neural networks and χ² model testing (Wolf ). The class of kernel regression methods utilize sets of model realizations or spectroscopic galaxy samples to approximate the conditional probability density of a galaxy redshift conditional on photometric properties via specified kernel functions. Given this approximate probability density the non-linear mapping between a galaxy's photometric properties and its redshift can be established via a kernel regression that provides the conditional expectation value of a galaxy redshift (Nadaraya ; Watson ; Boris et al. ; Wang et al. ; Wolf ).

Artificial neural networks aim at fitting the non-linear relation between observed photometric parameters and the redshift of galaxies (Firth, Lahav & Somerville ; Collister & Lahav ). In this approach, the non-linear relation between photometry and galaxy redshifts is represented by a multi-layer perceptron neural network, trained by a sufficiently large training set for which both photometry and accurate redshifts are known (Collister & Lahav ). Generally, artificial neural networks provide a single unique parameter estimate. However, if the solution is degenerate, they tend to return the most likely one (Wolf ). It is also possible to estimate uncertainties of the inferred quantities by re sampling the input parameters from their error distributions. The outputs of an artificial neural network will then provide a probability distribution characterizing the uncertainty in the inferred quantity (Collister & Lahav ; Wolf ).

Additionally, Bayesian redshift inference methods have been proposed in the literature (Benítez ; Wolf ). These models rely on χ² model testing, assuming a parametrized model between photometric properties and redshifts of galaxies (Wolf ). Uncertainty information on the observations then permits us to determine a probability density distribution for the inferred quantity and to provide expectation values and parameter uncertainties (Wolf ).

Aside from the photometric problem, literature also proposes Bayesian methods to correct for redshift distortions introduced by peculiar motions of the observed objects (see e.g. Kitaura & Enßlin ; Weig & Enßlin ). Note that the method presented in this work is very general and can easily be adapted to also solve these kinds of problems.

In any case, improvements in photo-z determinations help us to improve the technique we describe here. We emphasize that our work begins where photo-z estimation leaves off.

Density Posterior

The density posterior distribution given in equation can be written as

by noting that the posterior distribution is conditionally independent of the observed redshifts once the true redshifts {z_p} are known. In order to relate this distribution to a log-normal Poissonian distribution we would like to express it in terms of galaxy counts rather than in terms of galaxy redshifts and positions in the sky. We can therefore write equation as the marginalization over all possible galaxy number count configurations {N_i}

where we assumed the conditional independence . Further, the number counts are defined in comoving Cartesian coordinates; we will therefore introduce the set of corresponding comoving Cartesian galaxy positions and marginalize over it

Since there exists a unique relation between the comoving Cartesian and the true redshift coordinates, we can write

with α_p and δ_p being the right ascension and declination of the pth galaxy, respectively. By introducing the vectors y_p as

the integral in equation can be solved to yield

Also note that the probability distribution can be written as

where δ^K(x) is the Kronecker delta, x_i is the position of the ith volume and is the counting kernel

with xⁱ being the components of the Cartesian vector x and K(x) being

and Δx is the length of the volume element. Therefore, introducing the galaxy number counts

we can solve equation to yield

Hence, the density posterior distribution can be described in terms of number counts.

Redshift Posterior Distribution

In this section we will derive the expression for the redshift posterior distribution of an individual galaxy. Making frequent use of Bayes rule we can write the conditional redshift posterior distribution described in Section 4 as

We will now consider the probability distribution , which basically describes a random draw of a particle out of an ensemble of particles. Changing the coordinates from spherical redshift coordinates to comoving Cartesian coordinates yields

where is the derivative of the comoving radial distance with respect to the redshift z. The Poissonian assumption for the galaxy distribution assumes that the number counts, and hence the number of particles inside individual cells, are independent. Further, the Poissonian model assumes the particles inside a volume element to be homogeneously distributed. We therefore yield

In the following we will omit the arguments of the number counts {M_i({ξ_p}, {u_p})} = {M_i}. We can then write

By noting that we can write

The probability is merely the probability of randomly picking a particle at a specific position conditional on the set of number counts. Hence we can write

where N_tot is the total number of particles and ΔV is the volume of the volume element. Next we will consider the fraction

where we used the result, obtained in Appendix , that the density posterior solely depends through the number counts on the set of galaxy coordinates. As can be easily seen by applying Bayes law we can write this fraction as the ratio of the Poissonian likelihoods

with being the intensity of the Poisson process. Since the number densities {N_i} and {M_i} differ only by one at the position x_n of the galaxy under consideration, the ratios of the product in equation are all equal to one except at this particular galaxy position. We can therefore write

Inserting all results into equation yields

Since the product in the third line of equation is only different from one at the galaxy position x_n we can write

Also note, since in this work we will use Metropolis–Hastings algorithms to explore the conditional redshift posterior distribution, we do not require to explicitly calculate the normalization constant of the probability distribution given in equation . Up to an overall normalization constant we can therefore write

This final result demonstrates that the conditional redshift posterior is essentially proportional to the intensity function of the Poissonian likelihood. It is also interesting to note that the right-hand hand side of equation does not depend on the galaxy coordinates {u_p} and {ξ_p} of all galaxies. This shows that given the underlying density field {s_i} the individual redshift posterior distribution is conditionally independent of all other galaxies

If we make the additional assumption that the statistical process by which the observed redshifts for individual galaxies is independent of the observed redshifts of all other galaxies we may write the redshift likelihood as

Using this definition in equation and omitting all normalization constants yield

This result permits us to sample each galaxy redshift z_n individually.