Bayesian photometric redshifts of blended sources

ABSTRACT

Photometric redshifts are necessary for enabling large-scale multicolour galaxy surveys to interpret their data and constrain cosmological parameters. While the increased depth of future surveys such as the Large Synoptic Survey Telescope (LSST) will produce higher precision constraints, it will also increase the fraction of sources that are blended. In this paper, we present a Bayesian photometric redshift (BPZ) method for blended sources with an arbitrary number of intrinsic components. This method generalizes existing template-based BPZ methods, and produces joint posterior distributions for the component redshifts that allow uncertainties to be propagated in a principled way. Using Bayesian model comparison, we infer the probability that a source is blended and the number of components that it contains. We extend our formalism to the case where sources are blended in some bands and resolved in others. Applying this to the combination of LSST- and Euclid-like surveys, we find that the addition of resolved photometry results in a significant improvement in the reduction of outliers over the fully blended case. We make available blendz, a Python implementation of our method.

methods: statistical, galaxies: distances and redshifts, cosmology: observations

1 INTRODUCTION

Current photometric surveys such as CFHTLens (Heymans et al. 2012), KiDS (de Jong et al. 2013), and DES (Dark Energy Survey Collaboration 2016) image galaxies over large volumes of the Universe to probe the growth of structure and the distribution of matter on large scales. Through techniques such as galaxy clustering and cosmic shear, these surveys are able to constrain cosmological parameters and conduct tests of the standard ΛCDM cosmological model (e.g. Kilbinger et al. 2013; Dark Energy Survey Collaboration 2018).

These observational tests require the redshift distribution of the sample to make model predictions for comparison. Additional information can also be obtained by considering the redshift dependence using tomography, where galaxies are placed into one of several redshift bins (e.g. Hu 1999; Petri, May & Haiman 2016; Joudaki et al. 2018). However, the large number of sources required to constrain cosmological parameters to high precision makes obtaining spectroscopic redshifts for the entire sample unfeasible. As a result, photometric redshifts are a vital part of the cosmological analysis pipeline of galaxy surveys.

Photometric redshift methods seek to infer the redshift of galaxies from noisy observations of their flux in several broad-band filters. They provide an alternative to spectroscopic redshifts that requires less telescope time, at the expense of a reduction in precision. As a result, photometric redshifts can be applied to galaxies too faint and samples too large for spectroscopic observations. There are two general classifications for photometric redshift methods that utilize flux information; template-based and empirical methods.

Template-based methods use a set of galaxy spectra that are assumed to be representative of every galaxy they are applied to. These templates are redshifted and integrated over the response function of each filter to produce predictions of observed fluxes. These predictions are then used to infer the redshift from the observed fluxes. Maximum likelihood methods (e.g. Bolzonella, Miralles & Pelló 2000; Ilbert et al. 2006) find the best-fitting template by minimizing χ² to estimate the redshift. Bayesian methods, introduced by Benítez (2000), marginalize over all templates to produce a posterior redshift distribution. This correctly accounts for the uncertainty in the galaxy template that is ignored by maximum likelihood methods. Bayesian methods also include prior distributions that can reduce catastrophic outliers.

Empirical methods estimate redshifts by fitting for the mapping between flux and redshift from a set of training data, rather than specifying it a priori through a template set. This mapping is typically found using machine learning methods such as neural networks (e.g. Collister & Lahav 2004; Sadeh, Abdalla & Lahav 2016), Gaussian processes (e.g. Way & Srivastava 2006; Almosallam, Jarvis & Roberts 2016), and random forests (e.g. Carliles et al. 2010; Carrasco Kind & Brunner 2013). These methods are examples of supervised learning; they require large data sets of fluxes and associated spectroscopic redshifts that are representative of the sample they are applied to.

If representative data are available, empirical methods are typically more accurate than template-based methods (Hildebrandt et al. 2010). However, redshift estimates of galaxies not represented by the training data are much less reliable (Beck et al. 2017). The common case where spectroscopic training data are shallower than the photometry can lead to biases where the redshifts of high-redshift galaxies are underestimated (Rivera et al. 2018).

In practice, template-based and empirical methods are not so distinct. The priors of Bayesian methods typically include a set of parameters that are fitted using a set of training data (e.g. Benítez 2000; Schmidt & Thorman 2013, also see Section 2.6). In addition, recent applications of photometric redshifts have used hybrid methods that combine a template-based approach with machine learning methods (e.g. Speagle & Eisenstein 2017; Duncan et al. 2018; Leistedt et al. 2018).

In addition to these methods, clustering redshift methods (e.g. Newman 2008; Ménard et al. 2013; Schmidt et al. 2013) cross-correlate the angular positions of a photometric sample with a spectroscopic sample to estimate the redshift distribution. Clustering redshift methods do not model fluxes as a function of redshift, instead only using the spatial information of photometric data. As such, clustering redshifts are complementary to other photometric redshift methods; Cawthon et al. (2017) use clustering redshifts to calibrate biases from other photometric redshift methods, for example.

Ensuring that photometric redshifts are accurate and precise is necessary for obtaining unbiased constraints on cosmological parameters. Huterer et al. (2006) found that future tomographic surveys would require the mean of each redshift bin to be known to a precision of 0.003, though this requirement can be reduced by self-calibration (e.g. Huterer et al. 2006; Sun, Zhan & Tao 2015; Samuroff et al. 2017) and combining weak lensing data with other cosmological probes such as baryonic acoustic oscillations (LSST Science Collaboration 2009). Photometric redshifts are also important in the calibration of other systematics. Multiplicative biases in the measurement of shear can be detected and corrected for, provided that photometric redshifts of galaxies in the sample are unbiased (Hoekstra, Viola & Herbonnet 2017). Weak lensing shape measurement biases can themselves also be redshift dependent; without unbiased redshift estimates to make corrections, these can lead to biases of a few per cent in the cosmological parameters σ₈ and |$w$|₀ (Semboloni et al. 2009).

Another key part of precision cosmology is an accurate understanding of uncertainties in parameter constraints. To enable this, uncertainties arising from each step of the analysis should be accounted for and propagated onwards. In cosmological analyses, this is typically accomplished using a Bayesian framework (e.g. Hildebrandt et al. 2017; Troxel et al. 2017), allowing these uncertainties to be combined and marginalized over for the final constraints. It is therefore essential that photometric redshift methods provide not only point estimates of redshifts, but also a measure of their uncertainties.

The uncertainty associated with a redshift estimate can be represented by a single number, i.e. a point estimate with an error bar. However, doing this necessitates making an assumption about how the error is distributed. Uncertainties in photometric redshifts can be highly non-Gaussian, and so are poorly described by a single number such as the variance. Photometric redshift methods that instead characterize their results using a probability distribution function (PDF) can capture all of this information.

Photometric redshifts can also suffer from degeneracies that result in high-redshift galaxies having similar colours to those at low redshifts (e.g. Graham et al. 2018). As a result, several well-separated redshifts are plausible, and an accurate representation of the uncertainty should reflect this. While this can be easily described with a multimodal PDF, a single number can be misleading. Error bars that cover the full range of parameter space between the low- and high-redshift estimates do not show that redshifts between these are disfavoured, inflating uncertainties.

Several photometric redshift methods are able to produce PDFs as their result. Bayesian template-based methods (e.g. Benítez 2000) produce a posterior distribution, a PDF of the model parameters conditioned on the data and any model assumptions. In addition to the galaxy redshift, these model parameters can include other quantities of interest such as the galaxy template (Feldmann et al. 2006). A joint posterior over all of these parameters contains information about the uncertainty of each, including any correlation between them. Machine learning methods can also produce PDFs by utilizing ensemble techniques, where the predictions of several models are combined to produce a distribution. Examples of this technique include the combination of decision trees in a random forest (Carrasco Kind & Brunner 2013) and committees of neural networks constructed with different network architectures and initialized randomly (Firth, Lahav & Somerville 2003).

Future galaxy surveys such as the Large Synoptic Survey Telescope (LSST; Ivezic et al. 2008) will obtain extremely high-precision constraints on cosmological parameters. By utilizing deeper photometry, these surveys will probe greater volumes than previously, resulting in an increased number density of galaxies imaged. While this increased depth drives the high precision these surveys will achieve, it also increases the fraction of objects that overlap with others along the line of sight, known as blending (Chang et al. 2013).

Most existing deblending methods do not utilize the colour information from photometry, instead using the spatial information contained in an image from a single band. The commonly used SExtractor (Bertin & Arnouts 1996) searches for adjacent pixels on a flux-thresholded map that separate into disjoint regions as the threshold is increased. Doing this for many thresholds allows each pixel to be assigned to a single object, contributing the entirety of its flux to that object. The SDSS deblender (Lupton 2005) lifts this restriction, splitting the flux proportionally between objects based on object templates. These templates are constructed by finding peaks in the image and assuming symmetry around them, comparing pairs of pixels and setting them to be equal. Profile fitting methods (e.g. Peng et al. 2002; Robotham et al. 2017) forward model the image using physical profile models, deblending by directly fitting for the galaxy properties. In far-infrared astronomy, blending is common due to the reduced angular resolution of these instruments compared to optical telescopes. As a result, galaxies that are well resolved in optical observations may become blended in the far-infrared. Deblending methods designed for this case such as Hurley et al. (2017) can use the unblended observations to place strong priors on the number and position of sources. We refer to this mix of blended and unblended observations as partial blending throughout this paper.

The ability for most deblending methods to successfully identify blended galaxies depends on their angular separation. Galaxies with too small an angular separation are instead identified as single sources. Dawson & Schneider (2014) estimate that |$45\hbox{--}55{{\ \rm per\ cent}}$| of sources in LSST will be blended, with |$15\hbox{--}20{{\ \rm per\ cent}}$| of all sources being misidentified as a single source, referred to as ambiguously blended objects.

Blending of sources can have an impact that is significant for constraining cosmological parameters. Dawson et al. (2016) estimate that ambiguously blended objects in LSST will result in an increase in shear noise of |$14{{\ \rm per\ cent}}$| for the deepest photometry (i < 27) and |$7{{\ \rm per\ cent}}$| for the gold standard sample (i < 25.3). Since these ambiguous blends are difficult to separate due to their small angular separation, deblending methods that incorporate colour information could be beneficial. Recent deblending methods such as MuSCADeT (Joseph, Courbin & Starck 2016) and scarlet (Melchior et al. 2018) incorporate this colour information by using wavelet transforms, enforcing that the representation of components is sparse in this space.

Deblending methods that produce a set of component-separated maps are useful for later applying existing analysis techniques designed for individual components to. However, splitting the analysis in this way can lose information, such as the correlation between deblending parameters and the parameters in a subsequent analysis. An analysis method that jointly constrains parameters directly from blended data provides a self-consistent, principled way to characterize and propagate this information.

In this paper, we present a method that generalizes the Benítez (2000) Bayesian photometric redshift (BPZ) method to the case of blended observations. This is a template-based method where the task of determining the component redshifts is cast as a Bayesian parameter inference problem. The product of such an inference is a joint posterior distribution of the redshift and magnitude of each component in the blended source. This distribution characterizes the complete statistical uncertainty in the result in a way that can be propagated through the rest of the cosmological analysis. Determining the number of components in an observed source, i.e. whether or not it is blended, is treated as a model comparison problem. In this way, our method allows the identification of blended sources from aperture photometry alone.

Throughout, we use source to refer to the (possibly) blended system that is observed, and component to refer to the underlying physical objects that make up this source. For parameters defined for each component in a source, we index over component using Greek letters and indicate the collection of these using sets, i.e. {θ} ≡ {θ_α, θ_β, …θ_N}. Vector quantities defined for each filter band are in bold |$\boldsymbol {q}$|⁠, and observed quantities are denoted with a hat |$\hat{q}$|⁠. Where necessary, quantities defined for a specific number of components are distinguished by a subscript number in brackets, i.e. q⁽¹⁾ is the definition of q for a single component. A summary of our notation is provided in Table 1.

Table 1.

A summary of the notation used throughout this paper.

Symbol	Description
N	Number of components
T	Number of templates
B	Number of filter bands
|$z$|α	Redshift of component α
m0, α	Reference band magnitude of component α
tα	Template index of component α
{|$z$|}	Set of redshifts of each component
{m0}	Set of reference band magnitudes of each component
{t}	Set of template indices of each component
b	Index over filter bands
b0	Index of reference band filter
|$\hat{F}_0$|	Observed flux in reference band
|$\hat{\boldsymbol {F}}$|	Vector of observed fluxes, excluding the reference band
σ0	Error on the reference-band flux
σb	Error on the flux in band b
|$F^{(1)}_{ t, b } (z, m_0 )$|	Model flux for a single component in band b, at redshift |$z$|⁠, with reference band magnitude m0 and templates t
|$F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace m_0\rbrace )$|	Model flux for N-component blended source in band b, at redshifts {|$z$|}, with reference band magnitudes {m0} and templates {t}
χ	Set of cosmological parameters Ωm, |$\Omega _\Lambda$|⁠, and H0
|$\xi ^{(N)}_\chi (\lbrace z\rbrace )$|	Combination of up to N-point correlation functions describing the extra probability of N galaxies jointly sitting at redshifts {|$z$|} due to clustering

This paper is organized as follows. In Section 2, we describe our formalism for estimating redshifts as a parameter inference problem, describing its application to partially blended systems in Section 3. In Section 4, we discuss our inference methods, detailing how we use model comparison to identify blended objects in Section 4.1. In Section 5, we test our method on simulated observations. Section 6 describes a test of our method on the Galaxy And Mass Assembly survey (GAMA; Baldry et al. 2017) blended sources catalogue (Holwerda et al. 2015), for which spectroscopic redshifts are available. We conclude in Section 7.

2 BLENDED PHOTO-Z FORMALISM

2.1 Flux model

In the same way as other template-based photometric redshift methods, we assume that each observed component is well represented by one of a set of T templates. Each template t is defined by its rest-frame spectral flux density F_t(λ_em) as a function of the emitted wavelength λ_em. This template is redshifted and observed through a broad-band filter b, the response of which is denoted W_b(λ_obs) as a function of observed wavelength λ_obs.

The flux of template t, at redshift |$z$| and observed in band b, is then given by

\begin{eqnarray*} T_{t,b}(z) = \frac{1}{c g^{\rm {AB}} C_b} \int _0^\infty F_{t} \left(\frac{\lambda }{1+z} \right) W_b(\lambda) \lambda \textrm{d} \lambda , \end{eqnarray*}

(1)

where g^AB = 3631 Jy is the zero-point of the AB-magnitude system and the normalization |$C_b \equiv \int _0^\infty \frac{W_b(\lambda)}{\lambda } \textrm{d} \lambda$|⁠. By including g^AB, our fluxes are dimensionless throughout, and the conversion between magnitudes and fluxes defined in the way is given by F ≡ 10^−0.4m. This template is then scaled by a normalization a so that the flux of an object modelled with template t, at a redshift |$z$| and observed in band b, is given by

\begin{eqnarray*} F^{(1)}_{t, b} (z, a ) = a T_{t, b} (z ) . \end{eqnarray*}

(2)

We model the flux of blended sources as a linear combination of individual component fluxes. For a blend of N components, the flux observed in band b is given by

\begin{eqnarray*} F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace a\rbrace ) = \sum _{\alpha =1}^{N} a_\alpha T_{t_\alpha , b} (z_\alpha ) , \end{eqnarray*}

(3)

where a_α is the normalization for component α. For the reasons specified in Section 2.3, we sample m_{0, α}, the apparent magnitude of each component in the reference band b₀ rather than this normalization directly. The normalization a_α is then defined such that the model flux in the reference band is equal to m_{0, α}. Thus, the model flux is given by

\begin{eqnarray*} F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace m_0\rbrace ) = \sum _{\alpha =1}^{N} \frac{10^{-0.4 m_{0, \alpha }}}{ T_{t_\alpha , b_0} (z_\alpha )} T_{t_\alpha , b} (z_\alpha ) . \end{eqnarray*}

(4)

2.2 Fully blended posterior

For a fixed number of components, photometric redshift determination is a parameter inference problem; we wish to infer the joint posterior distribution of the redshifts and apparent magnitudes of each component given a data vector |$\hat{\boldsymbol {D}}$| of B broad-band fluxes. This data vector is split into two parts |$\hat{\boldsymbol {D}} = (\hat{\boldsymbol {F}} , \hat{F}_0)$|⁠, where |$\hat{F}_0$| is the flux of the reference band and |$\hat{\boldsymbol {F}}$| is the vector of the remaining B − 1 fluxes. This is done since the normalization of each component is defined in the reference band, and it is the flux of this band on which the priors are conditioned.

Following BPZ (Benítez 2000), we set the flux of non-detections to zero. Likewise, bands that are not observed are given a flux of zero, with the corresponding error set to an extremely large value. As discussed in Section 2.4, we assume that sources are selected using a magnitude limit on a single selection band. We therefore require that the source is detected in this band, by definition.

We start by writing our desired posterior as a marginalization over templates for each component. For N components, we marginalize over sets of N template indices {t}_i = {t_α, t_β…t_N}_i. Each template index can take a value 1 ≤ t ≤ T and components may share the same template, so there are T^N of these sets to marginalize over, giving

\begin{eqnarray*} P (\lbrace z\rbrace , \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N) = \sum _{i=1}^{T^{N}} P (\lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N) .\nonumber\\ \end{eqnarray*}

(5)

We have emphasized that our posterior is defined for a fixed number of components by conditioning on N. In the general case where this number is unknown a priori, it can be inferred from the data; this is discussed in Section 4.1. We have also made the dependence on cosmological parameters, which are required for converting between distance and redshift, explicit in the above expression. These parameters are denoted by |$\chi = \lbrace \Omega _{\textrm{m}}, \Omega _\Lambda , H_0\rbrace$| for brevity. Applying Bayes rule, the posterior becomes

\begin{eqnarray*} &&{P (\lbrace z\rbrace , \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N)} \nonumber \\ &&{\quad\propto\, \sum _{i=1}^{T^{N}} P (\hat{\boldsymbol {F}}, \hat{F}_0 | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace , N ) P (\lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace | \chi , N).} \end{eqnarray*}

(6)

Since only the prior is dependent on cosmological parameters, we have removed the conditioning on χ from the likelihood. We then factorize the likelihood so that it is split in the same way as the data vector, giving

\begin{eqnarray*} &&{P (\lbrace z\rbrace , \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N) \propto \sum _{i=1}^{T^{N}} \,\, P (\hat{\boldsymbol {F}} | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace , N)}\nonumber\\ &&{\quad\times P (\hat{F}_0 | \lbrace m_0\rbrace , N) P (\lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace | \chi , N ).} \end{eqnarray*}

(7)

Since the magnitude of each component in the reference band is a sampled parameter in the posterior, our model for the reference-band flux is simply the sum of these after converting from magnitudes to fluxes. As a result, the conditioning on {|$z$|} and {t}_i in the reference band likelihood is unnecessary and so has been removed. We assume that the error on the observed reference-band flux is normally distributed with variance |$\sigma _{0}^2$|⁠. Thus, the reference band likelihood is given by

\begin{eqnarray*} P (\hat{F}_0 | \lbrace m_0\rbrace , N ) = \frac{1}{\sqrt{2 \pi \sigma _0^2}} \textrm{exp} \left[ \frac{ - \left(\hat{F}_0 - \sum _{\alpha =1}^{N} 10^{-0.4 m_{0, \alpha }} \right)^2 }{ 2 \sigma _{0}^2}\right], \nonumber\\ \end{eqnarray*}

(8)

where m_{0, α} is the sampled reference band magnitude for component α. Similarly, we use an uncorrelated multivariate Gaussian likelihood for |$\hat{\boldsymbol {F}}$|⁠,

\begin{eqnarray*} &&{P (\hat{\boldsymbol {F}} | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace , N )} \nonumber\\ &&{\quad=\, \prod _{b=1}^{B} \frac{1}{\sqrt{2 \pi \sigma _b^2}} \textrm{exp} \left[ - \frac{ \left(\hat{F}_b - F^{(N)}_{ \lbrace t\rbrace _i, b } \big (\lbrace z\rbrace , \lbrace m_0\rbrace \big) \right)^2 }{ 2 \sigma _{b}^2} \right] ,} \end{eqnarray*}

(9)

where |$F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace m_0\rbrace )$| is the model flux specified in equation (4) and |$\sigma _{b}^2$| is the variance on the observed flux in band b.

2.3 Separating the joint prior

We now develop the prior so that it can be written in terms of individual components. We start by separating the joint prior into a product over priors on redshift, template, and magnitude. Removing unnecessary conditioning, the joint prior becomes

\begin{eqnarray*} P (\lbrace z\rbrace , \lbrace t\rbrace , \lbrace m_0\rbrace | \chi , N ) &=& P (\lbrace z\rbrace | \lbrace t\rbrace , \lbrace m_0\rbrace , \chi , N ) \nonumber\\ &&\times \, P (\lbrace t\rbrace | \lbrace m_0\rbrace , N ) P (\lbrace m_0\rbrace | N). \end{eqnarray*}

(10)

This splitting up of the joint prior is similar to the approach of Benítez (2000). There are two important differences, however. First, we include a prior on the apparent magnitude of each component. This differs from the approach of Benítez (2000), who considers the magnitude on which the redshift and template priors are conditioned to be exactly the observed reference band magnitude. The uncertainty in the scaling of the template is then represented by marginalizing over a normalization factor with an assumed flat prior. However, while this normalization is not defined as such, it is acting to set the apparent magnitude of the source in the reference band. This magnitude is a quantity about which prior information is known.

The prior information on the apparent magnitude of components is particularly important in the blended case, as we need to consider more than just the overall magnitude of the source. The individual magnitudes of each component are necessary for scaling the model fluxes when predicting the model flux |$F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace m_0\rbrace )$|⁠. In addition, motivated by existing galaxy observations and following Benítez (2000), our redshift and template priors for each component are magnitude-dependent. The individual component magnitudes are not directly observed in the blended case, and must therefore be considered as random variables in our model.

An alternative to sampling the magnitudes directly would be to make the fraction each component contributes to the total flux a model parameter. However, the combination of intrinsic magnitude distributions and survey-specific selection effects would give the distribution of this fraction a highly complicated shape. Instead, including a prior on the magnitude of each component allows these effects to be easily accounted for.

The other important difference in the blended case is that each term in equation (10) is a joint prior over all components in the source. The redshift, type, and magnitude properties of individual galaxies are much more well studied than those of blended sources. To make use of this information, we write these joint priors in terms of priors on the individual components.

First, we assume that the template priors for each component are independent, i.e. galaxy types are not correlated. This allows us to split the template prior as

\begin{eqnarray*} P (\lbrace t\rbrace | \lbrace m_0\rbrace , N ) = \prod _{\alpha =1}^{N} P (t_\alpha | m_{0, \alpha } ) . \end{eqnarray*}

(11)

We also make the assumption that the redshift of each component depends only on its own type, not the types of other components. The redshifts of each component cannot be assumed to be independent however, as galaxies are distributed in a correlated way. The additional probability of finding N galaxies within a separation r over a random Poisson process is described by galaxy correlation functions of up to order N (Peebles 2001). We denote the combination of correlation functions describing this extra correlation as |$\xi ^{(N)}_\chi (\lbrace z\rbrace )$|⁠, i.e. the excess probability for two galaxies is given by

\begin{eqnarray*} 1 + \xi ^{(2)}_\chi (z_\alpha , z_\beta) \equiv 1 + \xi (r_{\alpha \beta }) , \end{eqnarray*}

(12)

where the separation |$r_{\alpha \beta } \equiv |\vec{r}_\alpha - \vec{r}_\beta |$| is the comoving distance between components α and β. In the two-component case, only the two-point correlation function ξ(r) is necessary. However, for three galaxies, higher order correlation functions are needed, i.e.

\begin{eqnarray*} 1 + \xi ^{(3)}_\chi (z_\alpha , z_\beta , z_\gamma) &\equiv& 1 + \xi (r_{\alpha \beta }) + \xi (r_{\beta \gamma }) + \xi (r_{\alpha \gamma }) \nonumber\\ &&+\, \zeta (r_{\alpha \beta }, r_{\beta \gamma }, r_{\alpha \gamma }), \end{eqnarray*}

(13)

where ζ(r_αβ, r_βγ, r_αγ) is the connected three-point galaxy correlation function.

The excess probability term |$\xi ^{(N)}_\chi (\lbrace z\rbrace )$| is defined in the posterior as a function of the component redshifts {|$z$|}, though the galaxy correlation function ξ (and higher order correlations) is defined in terms of comoving separation r. We therefore need to convert between the redshifts of each component and the comoving distance separating them. The line-of-sight comoving distance as a function of redshift is given by (e.g. Hogg 1999)

\begin{eqnarray*} r(z) = \frac{c}{H_0} \int _{0}^{z} \frac{\textrm{d} z^{\prime }}{E(z^{\prime })} , \end{eqnarray*}

(14)

where, neglecting radiation density and rewriting Ω ≡ Ω_m + Ω_Λ,

\begin{eqnarray*} E(z) = \sqrt{\Omega _{\textrm{m}} (1+z)^3 + (1 - \Omega) (1+z)^2 + \Omega _\Lambda } . \end{eqnarray*}

(15)

We assume a flat Planck¹ (Planck Collaboration XIII 2016) cosmology throughout; Ω_m = 0.3065, |$\Omega _\Lambda = 0.6935$|⁠, and |$H_0 = 67.9 \rm{\,km} \rm{\,s}^{-1} \rm{\,Mpc}^{-1}$|⁠.

However, the comoving distance separating components will depend not only on their redshifts, but also on their angular separation on the sky. As a result, we derive an effective correlation function ξ_eff that takes this angular dependence into account.

Consider the case of a two-component blend, as shown in Fig. 1. The two components are at comoving distances r_α and r_β from the observer, with separation Δr ≡ r_β − r_α. From the definition of the correlation function, we can write the ratio of the expected number of galaxies in a region with clustering N^ξ and that without N⁰ as

\begin{eqnarray*} 1 + \xi _{\rm {eff}} = \frac{N^{\xi }}{N^{0}} . \end{eqnarray*}

(16)

Figure 1.

Diagram showing the set-up of the ξ_eff calculation. We assume that two galaxies, represented by grey circles, will be blended if their angular separation is within θ. Given that these two galaxies are blended, the galaxy at a comoving distance r_β will lie within the disc.

Given that these components are blended, there is some maximum angular separation θ between them; we assume this to be small. We therefore compare the expected number of galaxies in a disc of width dr and radius ρ_max = r_βθ. The expected number without clustering is given by

\begin{eqnarray*} N_{\rm {disc}}^0 = \bar{n} \pi r_\beta ^2 \theta ^2 . \end{eqnarray*}

(17)

To find the expected number with clustering, we integrate over the disc using the volume element of an annulus with radius ρ, i.e.

\begin{eqnarray*} N_{\rm {disc}}^\xi = \int _{\rho =0}^{r_\beta \theta } \bar{n} [1 + \xi (r)] 2 \pi \rho \mathrm{d} \rho \mathrm{d} r . \end{eqnarray*}

(18)

Thus, writing |$r = \sqrt{\Delta r^2 + \rho ^2}$|⁠, the ratio becomes

\begin{eqnarray*} 1 + \xi _{\rm {eff}} = \frac{2}{r_\beta ^2 \theta ^2} \int _{\rho =0}^{r_\beta \theta } \left[1 + \xi \left(\sqrt{\Delta r^2 + \rho ^2} \right) \right] \rho \mathrm{d} \rho . \end{eqnarray*}

(19)

As described below, the effect of clustering is small. As a result, we adopt a simple power law for the two-point correlation function,

\begin{eqnarray*} \xi (r) \propto \left(\frac{r}{r_0} \right)^{-\gamma } . \end{eqnarray*}

(20)

Inserting this into equation (19) and integrating, the effective correlation function is given by

\begin{eqnarray*} \xi _{\rm {eff}} (r_\alpha , r_\beta) = \frac{r_0^2}{\left(1 {-} \frac{\gamma }{2} \right) r_\beta ^2 \theta ^2}\!\! \left[ \left(\frac{\Delta r^2 {+} r_\beta ^2 \theta ^2}{r_0^2} \right)^{1 {-} \frac{\gamma }{2}} - \left(\frac{\Delta r^2}{r_0^2} \right)^{1 - \frac{\gamma }{2}}\right].\nonumber\\ \end{eqnarray*}

(21)

The effect of the strength of clustering evolving with redshift can be included in this formalism by allowing the parameters r₀ and γ to vary with redshift (e.g. Sołtan 2016). We test the effect of this on the redshift inference by using a toy model where γ = 1.92 is kept constant, while r₀ linearly varies between |$r_0 = 5 \rm{\,Mpc} \,h^{-1}$| at redshift |$z$| = 2 and |$r_0 = 6 \rm{\,Mpc} \, h^{-1}$| at redshift |$z$| = 0.5, with a linear extrapolation outside of this range. Since the value of ξ_eff is non-negligible only when |$z$|_α ≈ |$z$|_β, this interpolation of r₀ is evaluated using |$z$|_α only.

We then simulated two-component blends from a prior with ξ_eff included as described in Section 5. Results assuming ξ_eff = 0 showed negligible differences from those where the effect was included. At the population level, the RMS scatter defined in equation (45) changed by |$0.205{{\ \rm per\ cent}}$| between results including and excluding the correlation function. There were also negligible changes to the results at the individual source level. A comparison of maximum a posteriori results in each case is shown in Fig. 2. The vast majority of sources show negligible differences, and visually inspecting the posteriors with larger changes shows these are highly multimodal, with modes of comparable heights. In these cases, small differences in the posteriors result in larger differences in point estimates as the maximum a posteriori value moves between modes. This is a limitation of point estimates, and can be mitigated by using the full information content of the posterior distributions, which do not vary strongly.

Figure 2.

Comparison of the maximum a posteriori point estimates including the effective correlation function and neglecting it, for sources simulated from a prior that includes it. The lower redshift components |$z$|_α are plotted with closed blue markers, and |$z$|_β are plotted with open green markers. Most sources show negligible differences, while sources that show large differences are multimodal. In these sources, small differences in the posterior result in point estimates moving between modes of slightly different heights, illustrating a limitation of point estimates.

Due to the small effect, our results throughout include a simple non-evolving correlation function with |$r_0=5 \rm{\,Mpc} \, h^{-1}$| and γ = 1.77 (Peebles 2001). A plot of this is given in Fig. 3.

Figure 3.

Plot of the effective correlation function ξ_eff versus Δ|$z$| ≡ |$z$|_β − |$z$|_α for various |$z$|_α used for the results throughout.

Inserting the correlation function allows us to write the joint redshift prior separated by component as

\begin{eqnarray*} P (\lbrace z\rbrace | \lbrace t\rbrace , \lbrace m_0\rbrace , \chi , N ) = \left[ 1 + \xi ^{(N)}_\chi (\lbrace z\rbrace ) \right] \prod _{\alpha =1}^{N} P (z_\alpha | t_\alpha , m_{0, \alpha}).\nonumber\\ \end{eqnarray*}

(22)

We separate the joint magnitude prior by assuming that the only correlation between the component magnitudes is from the effect of a selection function S({m₀}) applied to the total magnitude, as discussed in Section 2.4. The magnitude prior can then be written as

\begin{eqnarray*} P (\lbrace m_0\rbrace | N ) = S (\lbrace m_0\rbrace ) \prod _{\alpha =1}^{N} P (m_{0, \alpha } ) , \end{eqnarray*}

(23)

Finally, we impose a sorting condition. Without this, the components would be exchangeable, i.e. swapping the component labels α, β… would have no effect on the prediction of the model. As a result, the marginalized posterior for the redshift of a single component would contain contributions from every component in the source, as demonstrated in Fig. 4.

Figure 4.

By not imposing a sorting condition, the components in a source are exchangeable. This is demonstrated here for a simple two-component blend with redshifts |$z$|_α = 0.31, |$z$|_β = 1.19 as indicated by the orange lines. As a result of the exchangeability, the 2D marginal redshift distribution is symmetric about the dashed black line, and each 1D posterior contains a distinct peak for each component.

Imposing a sorting condition on either the magnitudes or the redshifts would have the same effect of breaking the exchangeability of the components. In our tests, sorting by redshift produced posteriors that recovered the true redshift more successfully. However, in high-redshift samples, there is an intrinsic colour degeneracy that can occasionally cause problems with a redshift sorting condition.

The Lyman break and Balmer break are absorption features occurring at 912 and 3650 Å, respectively. If photometry over a sufficiently wide wavelength range is not available, a Lyman break at high redshift can be confused with a Balmer break at low redshift (e.g. Graham et al. 2018). If the sample is deep enough that these high-redshift solutions are not unlikely a priori, this can cause bimodal posteriors and contribute to catastrophic outliers (Brimioulle et al. 2008).

Consider the case of a two-component blend where the redshift of one component is well constrained but the other has a bimodal posterior. If the well-constrained redshift happens to lie between these two modes, it will appear in the 1D marginal distributions of each component redshift, as whether it is the lower or higher redshift object depends on which of the two degenerate peaks is being sampled. In this case, sorting by magnitudes would result in a posterior more representative of the underlying system, where the redshift of one component is well constrained while the other has two well-separated modes. We did not find this to be a problem in our tests however, and so apply redshift sorting throughout.

The sorting condition Λ_α is imposed by introducing Heaviside step functions Θ into the product over components, and is defined as

\begin{eqnarray*} \Lambda _\alpha &=& 1 \quad \textrm{for} \quad \alpha = 1 \nonumber\\ &=& \Theta (q_{\alpha -1} - q_{\alpha }) \quad \textrm{otherwise,} \end{eqnarray*}

(24)

where q is either |$z$| or m₀ depending on whether redshift or magnitude sorting is used. In summary, the posterior for the fully blended case is given by

\begin{eqnarray*} &&{P (\lbrace z\rbrace , \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N) \propto \sum _{i=1}^{T^{N}} \,\, P (\hat{\boldsymbol {F}} | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace , N )}\nonumber\\ &&{\quad\times \, P (\hat{F}_0 | \lbrace m_0\rbrace ) \left[ 1 + \xi ^{(N)}_\chi (\lbrace z\rbrace ) \right] S (\lbrace m_0\rbrace )}\nonumber \\ &&{\qquad\times \,\prod _{\alpha =1}^{N} \Lambda _\alpha P (z_\alpha | t_\alpha , m_{0, \alpha } ) P (t_\alpha | m_{0, \alpha } ) P (m_{0, \alpha}).} \end{eqnarray*}

(25)

2.4 Accounting for selection effects

When considering the total apparent magnitude of a source, we must account for the selection effect of the survey observing it. Galaxy surveys typically select sources by imposing cuts on the apparent magnitude they observe m < m_lim since they cannot observe arbitrarily faint sources. As we are sampling intrinsic magnitudes rather than observed magnitudes, these selection effects do not impose a hard cut in our magnitude prior.

Consider a source with an intrinsic apparent magnitude exactly equal to the survey magnitude limit. Assuming a normal distribution for the observational error, the probability of observing this source is 1/2, since its observed apparent magnitude is equally likely to have been scattered above and below the magnitude cut. However, since objects in the sample have been detected by definition, we know the source must have been scattered brighter, effectively breaking the symmetry of the error distribution. As a result, intrinsic apparent magnitudes around the magnitude limit are less probable and should be downweighted.

To account for this, we follow the approach described in Leistedt, Mortlock & Peiris (2016) for including a selection effect. A discrete variable D representing the fact that an object was detected is introduced, and each term in the posterior is conditioned on it. We assume that our selection effect is imposed on a single selection band. Without loss of generality, we derive the effect by assuming that the selection band is the reference band b₀ and so only the reference band likelihood is affected. Conditioning on D, the likelihood can be written using Bayes rule as

\begin{eqnarray*} &&{P (\hat{F}_0 | \lbrace m_0\rbrace , N, D )}\nonumber\\ &&{\quad = \frac{ P (D | \hat{F}_0, \lbrace m_0\rbrace , N ) P (\hat{F}_0 | \lbrace m_0\rbrace , N ) }{ \int _0^\infty P (D | \hat{F}_0, \lbrace m_0\rbrace , N ) P (\hat{F}_0 | \lbrace m_0\rbrace , N ) \textrm{d} \hat{F}_0}.} \end{eqnarray*}

(26)

The numerator of equation (26) is equal to the likelihood defined in equation (8) since the probability of detection for an object that we know has been observed is |$P (D | \hat{F}_0, \lbrace m_0\rbrace , N, ) = 1$|⁠. After integrating over |$\hat{F}_0$|⁠, the denominator depends only on {m₀} and represents the effect of the magnitude selection. We therefore choose to write this term as part of the joint magnitude prior, defining the selection effect

\begin{eqnarray*} S (\lbrace m_0\rbrace ) = \int _0^\infty P (D | \hat{F}_0, \lbrace m_0\rbrace , N, ) P (\hat{F}_0 | \lbrace m_0\rbrace , N ) \textrm{d} \hat{F}_0 \end{eqnarray*}

(27)

that appears in the posterior in equation (25). The selection is a hard cut based on the observed flux, and so

\begin{eqnarray*} P (D | \hat{F}_0, \lbrace m_0\rbrace , N, ) &=& 1 \quad \textrm{for} \quad \hat{F}_0 \gt 10^{-0.4 m_{\textrm{lim}}} \nonumber\\ &=& 0 \quad \textrm{otherwise.} \end{eqnarray*}

(28)

Thus, the integral becomes

\begin{eqnarray*} S (\lbrace m_0\rbrace ) = \int _{ 10^{-0.4 m_{\textrm{lim}} }}^\infty P (\hat{F}_0 | \lbrace m_0\rbrace , N ) \textrm{d} \hat{F}_0 . \end{eqnarray*}

(29)

Since the reference band likelihood is assumed Gaussian, this can be written in terms of the normal cumulative distribution function as |$S\left(\lbrace m_0\rbrace \right) = 1 - \Phi (\hat{F}_0)$|⁠, where Φ is defined for a Gaussian distribution with mean μ and standard deviation σ to be

\begin{eqnarray*} \Phi (x) = \frac{1}{2} \left[ 1 + \textrm{erf} \left(\frac{x - \mu }{\sigma \sqrt{2}} \right) \right] . \end{eqnarray*}

(30)

Inserting this into equation (29), the effect of the magnitude selection can be written as

\begin{eqnarray*} S (\lbrace m_0\rbrace ) = \frac{1}{2} - \frac{1}{2} \textrm{erf} \left(\frac{ 10^{-0.4 m_{\rm {lim}}} - \sum _{\alpha =1}^{N} 10^{-0.4{m_{0, \alpha }}}}{\sigma _{0} \sqrt{2}} \right) . \end{eqnarray*}

(31)

By replacing the reference-band flux |$\sum _{\alpha =1}^{N} 10^{-0.4{m_{0, \alpha }}}$| with the model flux |$F^{(N)}_{ \lbrace t\rbrace , b } (\lbrace z\rbrace , \lbrace m_0\rbrace )$|⁠, this selection can be performed on any band. The selection function would then also be dependent on the redshifts and templates, i.e. |$S (\lbrace z\rbrace , \lbrace t\rbrace , \lbrace m_0\rbrace )$|⁠. This choice of selection band is included in the implementation described in Section 4.3.

A plot of this selection function for a galaxy from the GAMA blended sources catalogue, described in Section 6, is shown in Fig. 5.

Figure 5.

Plot of the selection function for a typical source from the GAMA blended sources catalogue used in Section 6. The dashed line shows the magnitude limit for this source m_lim < 19.

2.5 Specifying the priors

Like all Bayesian methods, the choice of priors should be problem dependent. For ease of comparison, we use the parametric forms given by Benítez (2000) with an additional magnitude prior. However, we stress that this choice is not a necessary one for our method and any joint |$P (z, t, m_0)$| prior may be used.

The Benítez (2000) template and redshift priors are given by

\begin{eqnarray*} P (t | m_{0} ) = f_t \textrm{e}^{- k_t (m_0 - m_{\textrm{min}})} \end{eqnarray*}

(32)

and

\begin{eqnarray*} P (z | t, m_{0} ) \propto z^{\alpha _t} \textrm{exp} \left\lbrace - \left[ \frac{z}{ z_{0, t} + k_{\textrm{m},t}(m_0 - m_{\textrm{min}}) } \right]^{\alpha _t} \right\rbrace \,, \end{eqnarray*}

(33)

respectively, where m_min is the bright-end magnitude cut as described below. The parameters α_t, |$z$|_{0, t}, k_{m, t}, f_t, and k_t are set separately for early, late, and irregular template types. Their values are found using the procedure discussed in Section 2.6 and are listed in Table 2.

Table 2.

The maximum a posteriori values of the prior parameters for the GAMA blended sources catalogue after calibrating using 26 782 unblended sources.

Parameters	Early	Late	Irregular	Type-independent
αt	1.59	1.53	1.30	–
|$z$|0, t	0.016	0.019	0.066	–
km, t	0.048	0.048	0.022	–
kt	0.044	0.024	–	-
ft	0.45	0.51	–	–
ϕ	–	–	–	0.71

We use a magnitude prior given by

\begin{eqnarray*} P (m_0 ) \propto 10^{ \phi m_0} . \end{eqnarray*}

(34)

The value ϕ = 0.6 gives the expression for the expected galaxy number counts in a homogeneous, Euclidean universe (Yasuda et al. 2001), though we leave ϕ free to also be found using the procedure discussed in Section 2.6. This was found to be ϕ = 0.705, though the difference in results compared to fixing ϕ = 0.6 was negligible.

Since the selection effect applies to the total source flux, individual components may be fainter than the survey magnitude limit, and so unobservable outside of a blended source. As a result, an analytic magnitude prior is required to describe the distribution of component magnitudes so that it can be used at faint magnitudes, where observations of individual components are unavailable.

For the reasons discussed in Section 4.2, we also apply a hard minimum and maximum cut to each component redshift |$z$|_α and component magnitude m_{0, α}. This cut has little effect on the redshift priors, which already go towards zero at large redshift; the same is true of the magnitude prior at bright magnitudes. The faint end of the magnitude prior of the brightest component is also already forced towards zero by the selection function. This is because in a N-component blend, the flux of the brightest component must be at least 1/N that of the total source flux, by definition. The magnitudes of the other components are not constrained in this way however, and so this cut represents a sharp boundary in the prior.

In our tests, the results of the redshift estimation were not strongly dependent on the position of this faint-end magnitude cut m_max. However, the evidence calculation described in Section 4.1 is dependent on its position, as changing the position of the cut alters the prior volume integrated over in equation (42). As a result, the position of this cut must be decided; it defines the limit where a galaxy is considered to contribute to a blend, and is therefore problem-dependent.

In principle, one could consider a galaxy to be blended if another arbitrarily faint galaxy lies along the same line of sight. In practice however, observations have limited precision, and the flux of an extremely dim galaxy cannot be detected. In other words, a sufficiently dim galaxy should no longer be considered a blended component, but rather a contribution to the noise.

In practice, a simple method to set this cut is to fix it for the entire sample. However, the argument above suggests that this cut should be dependent on the noise of the observation, i.e. that m_max should be set to the faintest magnitude that would have an observable effect. Fixing m_max is effectively an assumption that the sample has sufficiently homogeneous noise properties that the change in this faintest magnitude is negligible. For a sample where this is not the case, the magnitude cut can be set as a nσ₀ flux deviation, i.e.

\begin{eqnarray*} m_{\textrm{max}} = - 2.5 \log _{10} (n \sigma _0) , \end{eqnarray*}

(35)

where σ₀ is the error on the reference-band flux that varies for each source. In the tests in Section 6, we test both of these methods of setting m_max.

2.6 Calibrating the priors using spectroscopic information

The joint prior is conditioned on a set of parameters θ, i.e. |$P (z, t, m_0 | \theta )$|⁠, the posterior distribution of which we wish to infer. We can use spectroscopic information of a sample of galaxies from the population of interest to calibrate the above priors as suggested by Benítez (2000).

We assume here that this calibration is done with unblended galaxies, though this procedure can be extended to include blended galaxies too, provided that the number of components N is known a priori. In that case, the reference band magnitudes of each component would need to be included as a parameter in this model, and either sampled along with θ or marginalized out of the posterior analytically.

We consider a sample of G galaxies with photometry and spectroscopic redshifts |$\hat{z}_{s}$|⁠. These redshifts are assumed to be exact, i.e. we neglect the error on |$\hat{z}_{s}$|⁠. The set notation here now runs over each independently observed galaxy, not the blended components as before.

We start by writing this posterior as a marginalization over the photometric redshift model parameters for each galaxy and applying Bayes rule. Since the likelihood is independent of the prior parameters, we condition on θ in the prior only, giving

\begin{eqnarray*} &&{P (\theta | \lbrace \hat{z}_s\rbrace , \lbrace \hat{\boldsymbol {F}}\rbrace , \lbrace \hat{F}_0\rbrace ) \propto \int \textrm{d}^{G} \lbrace z\rbrace \int \textrm{d}^{G} \lbrace m_0\rbrace}\nonumber\\ &&{\quad\times \, \sum _{i=1}^{T^{G}} P (\lbrace \hat{z}_s\rbrace , \lbrace \hat{\boldsymbol {F}} \rbrace , \lbrace \hat{F}_0\rbrace | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace ) P (\theta , \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace ).}\nonumber\\ \end{eqnarray*}

(36)

We apply product rule to separate the joint prior and remove other unnecessary conditioning. We also assume that the galaxies in the sample are independent, and so all terms not shared across the population [i.e. P(θ)] can be written as a product over galaxies. The posterior then becomes

\begin{eqnarray*} &&{P (\theta | \lbrace \hat{z}_s\rbrace , \lbrace \hat{\boldsymbol {F}}\rbrace , \lbrace \hat{F}_0\rbrace ) \propto P (\theta ) \prod _{g=1}^{G} \int \textrm{d} z_g \int \textrm{d} m_{0, g}}\nonumber \\ &&{\quad \times \,\sum _{i_g =1}^{T} P (\hat{F}_{0, g} | m_{0, g} ) P (\hat{\boldsymbol {F}}_g | z_g, t_{i_g}, m_{0, g} )}\nonumber\\ &&{\quad\times \,P (\hat{z}_{s, g} | z_g ) P (z_g, t_{i, g}, m_{0, g} | \theta).} \end{eqnarray*}

(37)

By assuming that the spectroscopic redshifts are exact, the redshift likelihood can be written as a delta function, i.e. |$P (\hat{z}_{s, g} | z_g ) = \delta (z_{g} - \hat{z}_{s, g} )$|⁠. We also assume that the error on the reference band magnitude is negligible, allowing us to write |$P (\hat{F}_{0, g} | m_{0, g} ) = \delta (m_{0, g} - \hat{m}_{0, g} )$|⁠, where |$\hat{m}_{0, g} = -2.5 \log _{10} \left(\hat{F}_{0, g} \right)$| is the reference-band flux of galaxy g, converted to magnitudes. Replacing these likelihoods with delta functions, the marginalization can be done analytically using the sifting property of the delta function to give

\begin{eqnarray*} &&{P (\theta | \lbrace \hat{z}_s\rbrace , \lbrace \hat{\boldsymbol {F}}\rbrace , \lbrace \hat{F}_0\rbrace ) \propto P (\theta )}\nonumber \\ &&{\quad\times \,\prod _{g=1}^{G} \sum _{i_g =1}^{T} P (\hat{\boldsymbol {F}}_g | \hat{z}_{s, g}, t_{i_g}, \hat{m}_{0, g} ) P (\hat{z}_{s, g}, t_{i_g}, \hat{m}_{0, g} | \theta ).} \end{eqnarray*}

(38)

To find the prior parameters θ that maximize this posterior, we use L-BFGS-B (Byrd et al. 1995), a local optimization algorithm that approximates the Hessian of the objective function and optimizes the parameters subject to simple box constraints; we use these constraints to ensure our parameters are positive. This method requires first-order derivatives, which we approximate through a finite-difference method.

The result of this procedure is an estimate of the maximum a posteriori values of the prior parameters. Throughout this paper, we use these values in the priors directly. In principle, these parameters could form part of a hierarchical model and be marginalized out as nuisance parameters. However, this would significantly increase the dimensionality of the parameter space to be sampled and, thus, the computation time required for each source. Table 2 lists the values of these prior parameters GAMA test described in Section 6. A plot of samples drawn from the resulting prior is shown in Fig. 6.

Figure 6.

Plot of the prior found for the test on the GAMA blended sources catalogue after calibrating using 26 782 unblended sources. The dashed line in the bottom panel shows a magnitude limit of r < 19.8.

3 PARTIALLY BLENDED SOURCES

We can modify the formalism above for the case of sources for which every component does not contribute to every observation. We refer to these as partially blended sources. This can be the case when combining photometry from a wide range of wavelengths, e.g. optical and far-infrared observations. This partial blending may also occur for some sources observed in both ground-based and space-based surveys, as the latter does not suffer from atmospheric seeing and so can achieve a higher spatial resolution. An example of a pair of such surveys is LSST (Ivezic et al. 2008) and Euclid (Laureijs et al. 2011). Utilizing resolved photometry from Euclid could improve the precision of photometric redshifts of sources that are blended in the higher signal-to-noise observations of LSST. This possibility is explored using simulated observations in Section 5.2.

To generalize the method for this case, we introduce the measurement-component mapping δ_{α, m}, a N × N_m matrix, where N_m is the number of measurements, a generalization of the number of bands in the fully blended case. This measurement-component mapping acts as an indicator variable, consisting only of zeros and ones indicating whether a particular component is present in a particular measurement.

An example of such a matrix is given below. Consider data containing N_m = 6 photometric measurements of N = 2 components. The first four measurements are of individually resolved components, while the final two measurements are blended. In a typical use case, we might expect the resolved measurements of each component to share filter bands, though the model does not require this. In this example, the measurement-component mapping is given by

\begin{eqnarray*} \delta = {\left[\begin{array}{c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c}1 & 0 & 1 & 0 & 1 & 1 \\ 0 & 1 & 0 & 1 & 1 & 1 \\ \end{array}\right]}. \end{eqnarray*}

(39)

We can then write the blended flux of N components at a redshift |$z$| in measurement m as

\begin{eqnarray*} F^{(N)}_{ \lbrace t\rbrace , m, \delta } (\lbrace z\rbrace , \lbrace m_0\rbrace ) = \sum _{\alpha =1}^{N} \delta _{\alpha , m} \frac{10^{-0.4 m_{0, \alpha }}}{ T_{t_\alpha , b_0} (z_\alpha )} T_{t_\alpha , m} (z_\alpha ) . \end{eqnarray*}

(40)

The only modification to the posterior of the fully blended case needed to accommodate the partial blending is to the sorting condition. As described in Section 2.3, the purpose of this condition is to prevent the exchangeability of components. However, this is not necessary in the partially blended case. Here, the components are intrinsically different as they appear individually in separate measurements and so are not exchangeable. As a result, we drop the sorting condition for the partially blended case, i.e. Λ_α = 1 over the entire parameter space. The posterior for the partially blended case is then given by

\begin{eqnarray*} &&{P (\lbrace z\rbrace , \lbrace m_0\rbrace | \hat{\boldsymbol {F}}, \hat{F}_0, \chi , N, \delta ) \propto \sum _{i=1}^{T^{N}} P (\hat{\boldsymbol {F}} | \lbrace z\rbrace , \lbrace t\rbrace _i, \lbrace m_0\rbrace , N, \delta )}\nonumber \\ &&{\quad\times \, P (\hat{F}_0 | \lbrace m_0\rbrace ) \left[ 1 + \xi ^{(N)}_\chi (\lbrace z\rbrace ) \right] S (\lbrace m_0\rbrace )}\nonumber\\ &&{\quad\times \,\prod _{\alpha =1}^{N} P (z_\alpha | t_\alpha , m_{0, \alpha } ) P (t_\alpha | m_{0, \alpha } ) P (m_{0, \alpha } ).} \end{eqnarray*}

(41)

4 INFERENCE USING NESTED SAMPLING

4.1 Determining the number of components with model comparison

The posteriors in equations (25) and (41) are defined for a specific number of components N. In general however, this number of components is not known a priori. We therefore need a method to determine how many components are present in a source. Since our model is defined for a fixed number of components, we treat finding the number of components in a source as a model comparison problem.

Bayesian model comparison involves the calculation of the evidence |$\mathcal {Z}$|⁠, an integral over the product of the prior and the likelihood (e.g. Trotta 2008). Given a data vector |$\boldsymbol {d}$|⁠, a model m, and a set of model parameters {θ}, the evidence is defined as

\begin{eqnarray*} \mathcal {Z} \equiv P (\boldsymbol {d} | m ) = \int P (\boldsymbol {d} | \lbrace \theta \rbrace , m ) P (\lbrace \theta \rbrace | m ) \textrm{d} \lbrace \theta \rbrace . \end{eqnarray*}

(42)

This evidence term plays the role of the normalization of the posterior and so is typically ignored in parameter inference problems where this normalization is irrelevant. However, the evidence is the quantity of interest for model comparison problems. The ratio of the posterior probabilities of two models is proportional to the ratio of their evidences, a quantity known as the Bayes factor. By considering the number of components in a source as the model, we can write the relative probability of the source containing n components compared to m components as

\begin{eqnarray*} \mathcal {P}_{n,m} = \frac{P (N = n | \hat{\boldsymbol {F}}, \hat{F}_0 )}{P (N = m | \hat{\boldsymbol {F}}, \hat{F}_0 )} = \frac{P (\hat{\boldsymbol {F}}, \hat{F}_0 | N = n )}{P (\hat{\boldsymbol {F}}, \hat{F}_0 | N = m)} \frac{P (N = n )}{P (N = m )} . \nonumber\\ \end{eqnarray*}

(43)

Considering the cases of either isolated galaxies or blends of two components, the model prior ratio |$P (N = 2 ) / P (N = 1 )$| represents the probability that a galaxy will be blended. Dawson & Schneider (2014) estimate the number of sources observed by LSST that will be blended by convolving Hubble Space Telescope images with a Gaussian point spread function (PSF) like that of LSST. They found this number to be |$45\text{--}55{{\ \rm per\ cent}}$| of the total sources observed, with |$15\text{--}20{{\ \rm per\ cent}}$| of observed sources classified as catastrophic blends that would be identified as single sources by fitting a profile template to a galaxy image. Chang et al. (2013) estimate that the rejection of blended sources will reduce the number density of LSST sources by |$16{{\ \rm per\ cent}}$|⁠, though this estimate does not include the catastrophic blends of above. Studies such as these using existing high-resolution data or simulated observations can inform the blending prior ratio. Throughout this paper, we present results where this prior ratio is |$P (N = 2 ) / P (N = 1 ) = 1$|⁠, i.e. we do not prefer either of the blended or single-component models a priori, though this information can be trivially included.

4.2 Nested sampling using multinest

Calculating the evidence directly through numerical integration presents a difficult technical problem, particularly as the number of dimensions increases. To avoid this, we use Nested sampling (Skilling 2006), a Monte Carlo method for estimating the evidence while also sampling the posterior for parameter inference. Nested sampling reduces the problem of estimating the evidence to sampling a series of increasing likelihood thresholds, i.e. progressively smaller prior volumes nested within one another. Equation (42) can then be calculated using a one-dimensional quadrature integration method over this prior volume.

The computationally difficult part of the nested sampling algorithm is sampling a new point from within the potentially complicated boundary defined by the likelihood threshold. The multinest sampler (Feroz, Hobson & Bridges 2009) does this efficiently by sampling from a collection of ellipses approximating this boundary rather than the prior itself. This collection of ellipses is formed by performing a clustering analysis on a fixed-sized set of the previous samples, known as the live points. A new sample is drawn from these ellipses, replacing the lowest likelihood point that is removed and stored as a posterior sample. Samples are rejected until the likelihood boundary is respected, though this occurs less frequently than when naively rejection sampling the prior.

The use of multiple ellipses when sampling has another distinct advantage in that it naturally enables efficient sampling of multimodal posteriors, since each mode is assigned a separate ellipse while low-probability regions between these modes are avoided. Multimodality is a feature that can cause difficulties for Markov Chain Monte Carlo (MCMC) samplers, as moving from one mode to another requires a move across the low-probability region separating them. As a result, these samplers can fail to explore the full posterior distribution, instead sampling only a single mode. We expect our problem to exhibit this multimodal behaviour due to the degeneracies described in Section 2.3, and so require a sampling method suited to this case.

The need for nested sampling methods to sample from the prior imposes some constraints on our choice of prior. multinest natively samples from a unit side-length hypercube and these samples are transformed into samples of the prior using a prior transform function. However, due to the discrete marginalization over template, we cannot separate the posterior to define a prior transform function. As a result, we take the approach suggested by Feroz et al. (2009) of defining a uniform prior to sample from, and defining the ‘likelihood’ for multinest as our marginalized posterior.

This has two main effects. First, the sampling is likely to be less efficient, as the prior sampling step is not guided by the true prior, and so low-prior regions may be sampled frequently. Secondly, sampling from a uniform prior necessitates imposing a hard cut on the prior range of each parameter. Since the location of these cuts affects the value of the evidence |$\mathcal {Z}$|⁠, they should not be imposed thoughtlessly. At high redshift and bright magnitudes, the priors tend to zero, meaning that the exact positions of these cuts have negligible effect on the evidence. However, this is not the case for the faint end of the magnitude priors; setting this cut is discussed is Section 2.5.

4.3 blendz package

We have written a Python package blendz to perform the redshift inference of blended sources described in Sections 2 and 3, and the identification of the number of components using model comparison described in Section 4.1. The package supports analysis of blends with an arbitrary number of components using either the included or user-supplied template sets. The output of such an analysis is a set of samples from the joint posterior for each number of components considered, and an estimate of the Bayes factor for model comparison. The model comparison can then easily include a model prior through multiplication of the Bayes factor.

The package is also written in an object-orientated way, allowing the user to easily redefine the priors. While the supplied prior is used in this work with galaxies of either early, late, or irregular types, it is written to be calibrated and used with any number of possible types. For blended sources of more than two components, the excess probability term |$\xi _{\chi }^{(N)}$| is defined recursively to use the correct combination of two-point terms and assumes higher order correlations are negligible.

Documentation and instructions for installation can be found at http://blendz.readthedocs.io. The package can also be immediately installed from the official Python Package Index² by using the pip install blendz command. Finally, the source is available in a git repository hosted at https://github.com/danmichaeljones/blendz.

5 RESULTS FROM MOCK OBSERVATIONS

5.1 Fully blended sources

As an initial test of the method, we used a Monte Carlo simulation to create a set of mock photometric observations to test our method against. These mock observations simulate an optical survey using the six LSST optical filters u, g, r, i, |$z$|⁠, Y (LSST Science Collaboration 2009), with an r-band magnitude selection of m_lim = 24. We also applied hard cuts to the component magnitudes of m_min = 19 and m_max = 26. We then generated 1000 sources, each of which is a blend of two components in all bands. This was done by sampling a prior describing this distribution of objects using the MCMC sampler emcee (Foreman-Mackey et al. 2013) to generate the true parameters {|$z$|}, {t}, {m₀} for each simulated source. A plot of this prior distribution, plotted using corner.py (Foreman-Mackey 2016), is shown in Fig. 7.

Figure 7.

Corner plot of the prior sampled to create the mock catalogue. As described in the text, the bimodal shape of the marginal magnitude distributions is a result of both the selection effect and sorting components by redshift. The redshift sorting condition can be seen as a hard diagonal cut in the joint redshift distribution.

The effect of the selection function and the faint-end magnitude cut can be seen clearly in the two-peaked shape of the marginal distributions of m_{0, α} and m_{0, β}. The brighter-magnitude peak is a result of the selection function. In the single-component case, this would cause the prior to tend to zero at faint magnitudes. In the two-component case however, the magnitude priors of each component extend beyond m_lim as the selection effect is applied to the combined magnitude of both components. The brighter component in a two-component blend must, by definition, contribute at least half of the total flux. As a result, the selection effect prevents the magnitude of this component from being too faint. Since we impose the sorting condition on redshifts, and the brightest component in a source is not exclusively the lower redshift one, this action of the sorting condition causes the brighter peak in the marginal distributions of both m_{0, α} and m_{0, β}. If we instead impose the sorting condition on the magnitudes, these distributions become unimodal. Fig. 7 also shows the effect of the redshift sorting condition in the (⁠|$z$|_α, |$z$|_β) marginal distribution as a hard diagonal cut.

The model fluxes for these sampled parameters were then generated using the template responses defined in Section 2.1. We use the template set of Coe et al. (2006), containing one early-type, two late-type, and one irregular-type templates from Coleman, Wu & Weedman (1980), two starburst templates from Kinney et al. (1996), and two starburst templates from Coe et al. (2006). This same template set is then used during the inference. This allows a test of the method without the effect of unrepresentative templates, a source of error that is not unique to the case of blended sources.

Finally, we add an observational error to each observation. The flux error in band b is randomly drawn from an uncorrelated, zero-centred normal distribution |$\sigma _{b} \sim \mathcal {N}\left(\sigma _{b} | 0, \Sigma \right)$|⁠. The noise is set for all sources to be the final 1σ depth expected from LSST (Ivezic et al. 2008). We use these noisy observations to draw samples from the one- and two-component posteriors to test both the redshift determination and model comparison performance, setting the prior to the true distribution the photometry was sampled from.

Fig. 8 shows two examples of the 4D posterior that is the output from our method for each sample. For plotting purposes, the number of live points used for sampling is larger than that used for the inference and model comparison results throughout this paper. However, the change in the results is negligible. The left-hand panel shows an example of a well-constrained source with a unimodal posterior. This posterior shows correlations between the component parameters; this is expected, since the total flux of each band that is well constrained by the observations is split between the components. Reducing the model flux in a band of one component will result in a compensation in the other component, correlating their parameters.

Figure 8.

The 4D posterior distribution output from our method for two example sources. The true parameter values are shown in orange. The left-hand panel shows a well-constrained source with some correlations between components, though the true redshift is well recovered. The right-hand panel shows an example of a bimodal posterior that can arise in photometric redshift problems.

The right-hand panel of Fig. 8 shows a particularly prominent example of the curved degeneracies that can arise in the blended posteriors. This is due to the total magnitude of the source being well constrained by noisy observations, while component magnitudes are not themselves observable. This leads to a degeneracy that is curved due to the non-linearity of adding magnitudes. A result of this curved degeneracy is bimodality in the marginalized posterior of |$z$|_β. However, there is still significant probability density around the true redshift, highlighting the importance of not compressing the information content of a full posterior distribution into only a small set of numbers.

Fig. 9 shows a comparison of the photometric estimation of the component redshifts against their true simulated values. Point estimates of the redshift |$z$|_MAP are obtained by taking the maximum a posteriori (MAP) value of each component redshift posterior, marginalizing over the other three parameters. The method recovers the true redshift of each component from simulated photometry well. The performance of photometric redshift methods is often summarized by the RMS scatter σ_RMS. We first define the normalized error for galaxy g as

\begin{eqnarray*} \tilde{\delta }z_g = \frac{\hat{z}_{s, g} - z_{\textrm{MAP}, g}}{1+ \hat{z}_{s, g}} , \end{eqnarray*}

(44)

where each galaxy g is a single component of a blended source. Writing the total number of galaxies in our test catalogue as N_g, we then define the RMS scatter as

\begin{eqnarray*} \sigma _{\textrm{RMS}} = \sqrt{\frac{1}{N_{g}} \sum _{g} \left(\tilde{\delta }z_g \right)^2 } . \end{eqnarray*}

(45)

Figure 9.

Scatter plot comparing the maximum a posteriori point estimates from the photometric redshift estimation with the true redshifts for the mock observations. The left-hand panels distinguish the components, with |$z$|_α plotted with closed blue markers and |$z$|_β plotted with open green markers. The centre panels show the blend identification, with sources identified as blends plotted with closed purple markers, and those misidentified as single sources plotted with open red markers. The right-hand panels show a 2D histogram of the combined sample. The panels in the top row show the results for the full mock catalogue, while the bottom row only includes sources where the standard deviation of samples from each redshift marginal posterior is sufficiently small, |$\sigma _\alpha \le 0.2 \forall \alpha$|⁠. The dashed lines in each panel show an error of 0.15(1 + |$z$|⁠).

Computing this quantity for our mock blended observations, we find an RMS scatter of σ_RMS = 0.163. This compares to a scatter of σ_RMS = 0.0267 for photometric redshifts of mock observations of single sources, demonstrating the added difficulty of blending.

This scatter can be improved by excluding sources with photometric redshifts that, using the uncertainty information of the posterior distribution, are identified as untrustworthy. This is done by comparing a summary statistic against a threshold that controls the stringency of the test; we use the standard deviation of redshift marginal-posterior samples σ_α separately for each component, though a variety of summary statistics are available. Keeping only sources with |$\sigma _\alpha \le 0.2 \forall \alpha$|⁠, the RMS scatter is reduced to σ_RMS = 0.064, with |$37{{\ \rm per\ cent}}$| of sources removed. The effect of this is shown in the bottom row of Fig. 9.

The percentage of outliers can also be quantified. Outliers are defined as sources where either component has an error |$|z_{\textrm{MAP}} - \hat{z}_s| \ge 0.15 (1 + \hat{z}_s)$|⁠. For the full set of mock observations, this percentage was found to be |$18.6{{\ \rm per\ cent}}$|⁠. By keeping only sources with |$\sigma _\alpha \le 0.2 \forall \alpha$| as described above, the percentage of outliers falls to only |$6.0{{\ \rm per\ cent}}$|⁠.

The results of the detection of blends are also shown in the centre panels of Fig. 9. By using equation (43), we calculate |$\mathcal {P}_{2, 1}$|⁠, the relative probability that a source is a two-component blend compared to a single source. The interpretation of this probability is problem-dependent; a probability of |$\ln (\mathcal {P}_{2, 1} ) \gt 0$| indicates a preference towards the source being blended, while a threshold to |$\ln (\mathcal {P}_{2, 1} ) \gt 5$| indicates strong evidence (Kass & Raftery 1995). Likewise, probabilities of |$\ln (\mathcal {P}_{2, 1} ) \lt 0$| and |$\ln (\mathcal {P}_{2, 1} ) \lt -5$| indicate a preference and strong evidence for the single source case, respectively. As the blended and single thresholds are pushed more positive and negative, respectively, there are sources with values of |$\ln \left(\mathcal {P}_{2, 1} \right)$| that fall between these thresholds. In these cases, the source is assigned neither label.

As described in Section 4.1, we assume the relative prior probability of a blend to be P(N = 2)/P(N = 1) = 1, i.e. we give no preference to either model. Under this assumption, the method identifies |$92.7{{\ \rm per\ cent}}$| of sources as blends and |$7.3{{\ \rm per\ cent}}$| as single sources. Increasing the threshold to strong evidence, we find that |$89.9{{\ \rm per\ cent}}$| of sources are identified as blends and |$0.2{{\ \rm per\ cent}}$| as single sources; the remaining |$9.9{{\ \rm per\ cent}}$| fall between these thresholds. The distribution of the relative probability of blending and the effect on blend identification of changing the threshold are shown in Fig. 10.

Figure 10.

The left-hand panel shows the distribution of the relative blend-to-single probability for the mock catalogue, with the inset showing the same distribution, zoomed around lower relative probabilities and binned more finely. The right-hand panel shows the percentage of sources assigned as either blended, single sources, or not assigned to either as the threshold for deciding between each label is changed.

These results show that the method can both recover the redshifts from broad-band observations of blended objects, and detect the blending of a large fraction of these objects from their photometry alone. In addition, the output from these tests is not just point estimates of redshifts, but the full four-dimensional posterior distributions that capture the correlations between components that can be lost by working with component-separated maps. These are the results of simulated observations, however; real data has the complication that the flux model is no longer exact, i.e. the templates are not perfectly representative of all galaxies observed. As such, we test the method on real data in Section 6.

5.2 Partially blended sources

To test the effect of adding resolved data, we created a set of mock photometric observations of two-component partially blended systems. These observations simulate the same six-band optical survey as described in Section 5.1, combined with a four-band optical and infrared space-based survey using the Euclid filters |$vis, Y, J, H$| (Racca et al. 2016). This latter survey is assumed to have made resolved measurements of each component, while the former is fully blended as before. Thus, our partially blended data vector contains 14 fluxes for each source.

For a comparison with the partially blended results, we repeat the inference several times. First, we compare against the fully blended LSST-like case described above. Next, we compare to an inference using the resolved Euclid bands only, testing the effect of removing the difficulty of blending but using lower signal-to-noise data. Finally, we test against the case of using both the LSST- and Euclid-like data, but assuming that sources are blended in all bands. This allows us to separate the improvement as a result of adding resolved data from that of simply having more bands available.

For the fully blended bands, we reuse the simulated fluxes described in Section 5.1. For the resolved bands, we generate observed fluxes using the same randomly sampled source parameters. The observed fluxes are then generated using the flux model described in Section 2.1 with added observational errors drawn randomly from an uncorrelated, zero-centred normal distribution. The noise in these resolved bands is set to the final 1σ depths expected from Euclid observations (Laureijs et al. 2011).

We use the same prior as described in Section 5.1 for both the simulation and inference steps. The reference band over which the prior is defined is set to be the rband of the blended observations. This band is not present during the inference step using only the resolved data. As a result, we use the flux model from Section 2.1 to convert between r- and Y-band magnitudes before evaluating the prior.

Fig. 11 shows a comparison of the photometric redshift point estimates with the true simulated values for the four sets of inferences. As before, these point estimates are the MAP values of the marginal redshift posteriors.

Figure 11.

Scatter plot comparing the maximum a posteriori point estimates for the fully blended, resolved, and partially blended cases. The closed blue markers represent the redshift of the closer component, |$z$|_α, while the open green markers represent the redshift of the more distant component, |$z$|_β.

The left-hand panel shows the fully blended case, the same results as Section 5.1. Again, the redshifts of many components are well recovered, though there is a significant fraction of outliers. The RMS scatter in this case was found to be σ_RMS = 0.163. The percentage of outliers, sources where either component has an error |$|z_{\textrm{MAP}} - \hat{z}_s| \ge 0.15 (1 + \hat{z}_s)$|⁠, was found to be |$18.6{{\ \rm per\ cent}}$|⁠.

The centre-left panel of Fig. 11 shows the results for the resolved observations. Though finding the photometric redshift of resolved sources is an easier inference problem, this is counteracted by the significant reduction in the signal-to-noise of this data. As a result, we find an RMS scatter of σ_RMS = 0.212 with |$55.0{{\ \rm per\ cent}}$| of sources marked as outliers.

The centre-right panel of Fig. 11 shows the results for combination of LSST- and Euclid-like data in the fully blended case. We find that the addition of the four Euclid bands significantly improves the precision of the redshift inference, which has an RMS scatter of σ_RMS = 0.073. The fraction of outliers has also improved significantly, with only |$6.6{{\ \rm per\ cent}}$| of sources marked as outliers.

Finally, the right-hand panel of Fig. 11 shows the results for the partially blended case, combining the high-precision blended observations with the resolved data. Here, we find that the RMS scatter has reduced to σ_RMS = 0.065, a factor of 2.5 improvement over the blended LSST-like data alone, and a factor of 1.12 over the combined LSST- and Euclid-like blended data. The percentage of outliers has also been reduced. Here, only |$3.4{{\ \rm per\ cent}}$| of sources are found to be outliers, a factor of 5 improvement over the fully blended case, and a factor of 1.9 improvement over the combined LSST- and Euclid-like blended data.

While the most significant improvement was obtained through the increase in the number of bands, these results show that the quality of photometric redshifts of blended sources can be improved through the inclusion of resolved data. This is particularly apparent in the reduction of outliers. One advantage conferred by the addition of resolved data is a constraint on the relative magnitudes of blended components. In the fully blended case, the reference-band magnitude of each of the components must be inferred from the combined magnitude of the blended source only. This can lead to the degenerate distributions shown in Fig. 8. Adding resolved photometry can help to break this degeneracy by providing information about each component individually. The precision of the photometric redshift inferences is therefore improved.

An example of this phenomenon is shown in Fig. 12. The left-hand panel shows a corner plot of the posterior distribution for a fully blended source. The marginal distributions for each component redshift are highly multimodal, with well-separated redshifts occurring at distinct magnitudes. Though there is significant posterior density around the true redshifts, the MAP point estimate of |$z$|_β would show a significant error, as the incorrect mode has a higher posterior. The right-hand panel of Fig. 12 shows the same source analysed in the partially blended case after the addition of the resolved photometry. Here, the width of the posterior has been significantly reduced by the removal of the incorrect mode. The posterior now shows that the redshift of the source has been well constrained, and the redshift point estimates would no longer have a large error.

Figure 12.

The 4D posterior distributions for a two-component blended source in the fully blended and partially blended cases. The left plot shows the result of inference using blended data only. While there is significant posterior density around the true parameter values shown by the orange line, this posterior is highly bimodal, with two distinct solutions that cannot be distinguished. The plot on the right shows the result of the partially blended case that includes both blended and resolved observations. The addition of information about the magnitude of each component separately has removed the incorrect mode, resulting in a posterior that recovers the true solution well.

6 GAMA BLENDED SOURCES CATALOGUE

The GAMA survey (Baldry et al. 2017) is a spectroscopic galaxy survey that observed 286 deg² of sky over several regions to a magnitude limit of between r < 19 and r < 19.8. In doing so, it obtained precise redshifts of |$\gt 150\,000$| sources. The observed regions were chosen to overlap with existing imaging surveys such as Sloan Digital Sky Survey (SDSS) (Stoughton et al. 2002) and VISTA Kilo-degree Infrared Galaxy (VIKING) Survey (Edge et al. 2013). As a result, the spectroscopic data are accompanied by a set of aperture-matched photometry covering nine filter bands u, g, r, i, |$z$|⁠, Y, J, H, K from optical to infrared wavelengths (Hill et al. 2011).

The GAMA blended sources catalogue (Holwerda et al. 2015) contains 280 sources from the GAMA survey that have been spectroscopically identified as blended objects. These were selected using an automated template-based spectrum fitting method (Baldry et al. 2014) that cross-correlates galaxy templates with the observed spectra to determine the galaxy redshift. Sources where two different redshifts showed strong cross-correlations were visually inspected, resulting in a selection of blended galaxies. The motivation of Holwerda et al. (2015) was the identification of strong lens candidates. However, a catalogue of spectroscopically identified blended galaxies with accompanying nine-band photometry gives us a useful test case for the blended photometric redshift estimation method on non-simulated photometry with secure redshifts available for both components.

We first calibrate the prior using the procedure described in Section 2.6. To do this, we used 26 782 unblended, well-observed galaxies. These were selected by enforcing every band to be free from SExtractor (Bertin & Arnouts 1996) error flags and excluding all galaxies in the blended source catalogue. The resulting prior from the calibration procedure is shown in Fig. 6. As discussed in Section 2.5, we test two methods of setting the faint-end magnitude cut m_max, first as a 5σ₀ flux deviation using equation (35), and secondly, fixing m_max = 20.8. Throughout, we refer to these as the sigma-m_max case and fixed-m_max case, respectively. We then proceed with the inference using the same template set³ described in Section 5.

The resulting redshift point estimates are shown in Fig. 13. While noisier than the simulated case, the method still recovers reasonable estimates; using equation (45), we find an RMS scatter of σ_RMS = 0.156 in both the sigma-m_max and fixed-m_max cases. This compares to the scatter for a set of unblended GAMA sources of σ_RMS = 0.116. Obtaining this value even without the added complication of blending suggests a mismatch between the sources and the template set.

Figure 13.

Scatter plot comparing the maximum a posteriori point estimates from the photometric redshift estimation with the spectroscopic redshifts for sources from the GAMA blended sources catalogue. The left-hand panels distinguish the components, with |$z$|_α plotted with closed blue markers and |$z$|_β plotted with open green markers. The centre panels show the blend identification, with sources identified as blends plotted with closed purple markers, and those misidentified as single sources plotted with open red markers. The right-hand panels show a 2D histogram of the combined sample. The panels in the top row show the results for the sigma-m_max case, while those in the bottom row show the fixed-m_max case. The dashed lines in each panel show an error of 0.15(1 + |$z$|⁠).

We also compute the inferred blend probability |$\mathcal {P}_{2, 1}$| for these galaxies. The distribution of these probabilities is shown in Fig. 14. As described in Section 5, |$\ln (\mathcal {P}_{2, 1}) \gt 0$| and |$\ln (\mathcal {P}_{2, 1}) \gt 5$| show a preference and strong evidence for a blended source, respectively, while |$\ln (\mathcal {P}_{2, 1}) \lt 0$| and |$\ln (\mathcal {P}_{2, 1}) \lt -5$| show the same for the single source case. The distribution of the blend probability and the effect of the evidence threshold on blend identification are shown in Fig. 14.

Figure 14.

Plots showing the differences in the model comparison results between the two methods tested of setting the faint-end magnitude cut m_max, labelled the sigma-m_max and fixed-m_max cases. The left-hand panel shows the distribution of the relative blend-to-single probability for the mock catalogue, with the inset showing the same distribution, zoomed around lower relative probabilities and binned more finely. The solid line shows the sigma-m_max case, and the dashed line shows the fixed-m_max case. The two right-hand panels show the percentage of sources assigned as either blended, single sources, or not assigned to either as the threshold for deciding between each label is changed.

In our tests of the sigma-m_max case, |$71.6{{\ \rm per\ cent}}$| of sources showed a preference for being blended, with |$28.4{{\ \rm per\ cent}}$| preferring a single source. Increasing the threshold to strong evidence, these percentages fall to |$61.8$| and |$18.2{{\ \rm per\ cent}}$|⁠, respectively. Finally, the incorrectly identified single sources can be excluded entirely by increasing the threshold to |$\left|\ln \left(\mathcal {P}_{2, 1}\right)\right| \lt 12.5$|⁠, with |$50.7{{\ \rm per\ cent}}$| of sources identified as blends at this level.

The identification of blends was very similar in the fixed-m_max case. We found that |$71.1{{\ \rm per\ cent}}$| of sources showed a preference for being blended, and |$28.9{{\ \rm per\ cent}}$| preferred a single source. At the strong evidence threshold, |$60.4{{\ \rm per\ cent}}$| of sources are correctly identified as blends, with |$18.2{{\ \rm per\ cent}}$| misidentified as single sources. The threshold to exclude misidentified sources completely in the fixed-m_max case is |$\left|\ln \left(\mathcal {P}_{2, 1}\right)\right| \lt 13.9$|⁠, slightly higher than the sigma-m_max case. At this level, |$48.0{{\ \rm per\ cent}}$| of sources are still correctly identified as blends.

These results show that photometric redshift estimates can be obtained for blended sources, and that the method can identify many blended sources from just their broad-band photometry. By adjusting the threshold of the probability |$\mathcal {P}_{2, 1}$|⁠, blended sources can be selected in a way that trades off completeness and purity.

Several techniques for improving the scatter of photometric redshifts have been proposed, such as rest-frame template error functions (Brammer, van Dokkum & Coppi 2008), iterative methods to modify templates to be more representative (Feldmann et al. 2006), using clustering-based redshift estimation to calibrate systematic biases using galaxies (Gatti et al. 2018) and intensity mapping observations (Alonso et al. 2017), and constructing priors in terms of physical galaxy properties (Tanaka 2015). While an investigation of these methods is beyond the scope of this paper, they could also be applied while using this method. This could help to reduce the scatter of the blended photometric redshift estimates to a level necessary for future surveys, while retaining the full information of the posterior for accurate error propagation.

7 CONCLUSIONS

Blended sources will become far more common in future galaxy surveys than are found currently due to increases in the depth of photometry and as a result, the number density of galaxies. We present a BPZ method that generalizes the existing BPZ (Benítez 2000) method to the case of blended observations. We derive a posterior for the redshift and magnitude of each component that we sample to obtain estimates of the redshift. We also use this posterior in a model comparison procedure to infer the number of components in a source.

By doing this, the method is able to infer both the redshift of each component within a blended source, and identify that a source is blended from its broad-band photometry alone. The joint posterior distribution of the redshifts of all components in a blend provides a complete accounting of the correlations in the final result, information that can be lost when separating components and estimating redshifts for each component separately. This uncertainty information is essential for obtaining accurate uncertainties on cosmological parameters that rely on the photometric redshift estimates. A Python implementation of the method, blendz, is available to download.

By inferring the redshifts of components directly from their blended photometry, the method presented here is directly applicable to ambiguously blended objects that cannot otherwise be deblended. The partial-blending formalism described in Section 3 also enables the catalogue-level joint analysis of sources in space- and ground-based surveys such as Euclid and LSST. The complementarity of these surveys will allow cosmological parameters to be constrained more precisely than either survey could individually, and analysis of blended sources from their aperture photometry will be simpler than a joint pixel-level analysis (Rhodes et al. 2017).

The method presented here could also be combined with existing deblending methods that utilize the spatial information of images directly. These methods are complementary; image-based deblending methods are effective provided that components are sufficiently well separated. If this is not the case, there is too little spatial information to be able to separate components, and colour information is necessary. Combining these methods could allow future surveys to identify a greater proportion of blended sources, reducing their effects on cosmological constraints. Deblending methods that also incorporate colour information would need to be combined with this method more carefully however, as the colour information would be used twice and thus the blending probabilities would not be independent. This method could instead be extended to incorporate imaging data by constructing a forward model of the galaxy in each band and constraining both morphology and redshift simultaneously.

ACKNOWLEDGEMENTS

We thank Andrew Jaffe, Daniel Mortlock, and Josh Greenslade for helpful discussions, and the referee Michael Troxel for helpful comments. DMJ acknowledges funding from STFC through training grant ST/N504336/1. GAMA is a joint European–Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the SDSS and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT, and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 179.A-2004.

Footnotes

1

We use the TT + lowP + lensing  + ext values from Table 4.

2

https://pypi.org/

3

The templates used to fit the spectroscopic redshifts as described in Holwerda et al. (2015) do not cover the full wavelength range of the photometry. As a result, we do not use them for the photometric redshift inference.