BayeSED-GALAXIES I. Performance test for simultaneous photometric redshift and stellar population parameter estimation of galaxies in the CSST wide-field multiband imaging survey

The first method to generate mock photometric catalog is built by randomly draw samples from the parameter space of a particular SED model while under the constraints of some empirical statistical properties of galaxies. The sampling is performed with the same nested sampler MultiNest as in the Bayesian SED analysis mode of BayeSED. A sample from this catalog will be used in §4 to test the performance of redshift and stellar population parameter estimation in the case where the SED modeling is perfect, since exactly the same SED modeling method will be used in the Bayesian SED analysis of it.

The second method to generate mock photometric catalog is based on an SED library which is built by the post-processing of galaxies from a hydrodynamical simulation. This catalog will be used in §5 to test the performance of redshift and stellar population parameter estimation in the case where the SED modeling is imperfect, since the SED modeling method employed in the Bayesian SED analysis will be very different from the one used to built it.

We start from the rest-frame spectra of galaxies which are produced using the light-cone from the cosmological hydrodynamical simulation Horizon-AGN (Dubois et al., 2014). The computation of these spectra, which accounts for the complex star formation history and metal enrichment of Horizon-AGN galaxies, and consistently includes dust attenuation, is described in details by Laigle et al. (2019) and Davidzon et al. (2019). The dust attenuation of galaxies in Horizon-AGN simulation is modelled for each stellar particle, assumed to be a SSP, by using the gas metal mass distribution as an approximation of the dust mass distribution, assuming a constant dust-to-metal mass ratio (Laigle et al., 2019). Besides, to obtain the amount of extinction at a given wavelength, the Weingartner & Draine (2001) model of Milky Way dust grain with $R_{\rm V}=3.1$ and the prominent $2175{\rm\AA}$-graphite bump is employed for post-processing the simulated galaxies. As mentioned in Laigle et al. (2019), the overall attenuation curve becomes less steep and the bump tends to be reduced when summing up the contribution of all the SSPs to obtain the resulting galaxy spectrum. They also noticed that the averaged attenuation curve in Horizon-AGN simulation can not be well reproduced by either the model of Calzetti et al. (2000) or Arnouts et al. (2013). The more flexible form of DAL as given by Equation 9 is more likely to reproduce the attenuation curves of galaxies in the Horizon-AGN simulation. In order to isolate the possible differences in the convolution with filter response function, observational error modeling, and consideration of IGM absorption, we choose to convert their rest-frame spectra of mock galaxies to corresponding mock photometries with BayeSED, instead of using their virtual photometries directly¹¹1We have compared the two set of photometries careful and found that they are actually very similar, only with some differences at the very faint end.. The consideration for the effects of IGM absorption is the same as in §2.1.1.Therefore, the difference between the empirical statistics-based (§2.1) and hydrodynamical simulation-based photometric mock catalog are only driven by the different SFH, CEH, DAL of mock galaxies and their different distribution of redshift and physical parameters.

The modeling of realistic errors on the flux is crucial for a meaningful performance test of redshift and stellar population parameter estimation. Here, we introduce the method we have employed to compute flux errors of mock galaxies and perturb their fluxes accordingly. The flux error for a wavelength band $i$ with N$\sigma$ AB magnitude limit $m_{{\rm lim},i}=-2.5*{\rm log}(F_{\rm lim,i})+23.9$ is given by:

where the flux limit $F_{\rm lim,i}$ and the systematic flux error $\sigma_{F,i,{\rm sys}}$ are in unit of $\mu{\rm Jy}$. From the relation that magnitude error $\sigma_{m}\approx 1.08574/SNR$ and signal-to-noise ratio $SNR=F/\sigma_{F}$, we can obtain:

As in Cao et al. (2018), we assume a systematic magnitude error $\sigma_{m,{\rm sys}}=0.02$. The final mock flux is obtained by the original flux perturbed by a Gaussian noise $\epsilon\sim\mathcal{N}\left(0,\sigma_{F,i}^{2}\right)$. In practice, the magnitude limit may have a dispersion $\sigma_{m,{\rm lim},i}$ for galaxies with different sizes. So, the actually used magnitude limit is drawn from the Gaussian distribution $\mathcal{N}\left(m_{{\rm lim},i},\sigma_{m,{\rm lim},i}^{2}\right)$. In this work, we set $\sigma_{m,{\rm lim},i}=0.1$. We have generated three sets of mock catalog for CSST-like, Euclid-like, and COSMOS-like surveys, respectively. A summary of the adopted depths in all bands of the three surveys is shown in Table 2. The response functions and modeled relation between magnitude and magnitude error are shown in the panels of Figure 2 for the 7 CSST bands, 3 Euclid bands and 26 COSMOS bands, respectively. To separate the effects of observational errors on the accuracy of parameter estimation, we also generated another two sets of mock catalog without adding observational errors. In this case (the no noise case in Figure 2), the magnitude errors are all fixed to be $0.01$, but the photometries have not been perturbed accordingly.

To test the performance of BayeSED, we selected two sets of samples of galaxies with $K_{s}<24.7$ (Laigle et al., 2019) and $i+<25$ (Cao et al., 2018)²²2In Cao et al. (2018) and Zhou et al. (2022a, b), only the high-quality sources with the $SNR\geq 10$ in g or i band were selected to test the performance of photometric redshift estimation. This selection is not used in this work to obtain a more full picture for the performance of the CSST wide-field multiband imaging survey. from the empirical statistics-based mock catalog (§2.1) and the hydrodynamical simulation-based mock catalog (§2.2), respectively. The first set of samples are obtained directly with BayeSED combined with SED models with different complexity by using the method present in §2.1.2. The second sample is selected from the Horizon-AGN hydrodynamical simulation-based photometric catalogs for COSMOS-like configuration³³3https://www.horizon-simulation.org/PHOTOCAT/HorizonAGN_LAIGLE-DAVIDZON+2019_COSMOS_v1.6.fits which contains $789,354$ galaxies. We find that a sample with $10,000$ galaxies is large enough to obtain stable results for the performance tests as presented in §4 and §5. The redshift and magnitude distributions of the two samples are presented in Figure 3. When employing different SED models with different complexity and the same set of empirical statistics, the empirical statistics-based samples show some differences, especially the redshift distribution. This is likely due to the different mapping relations from physical parameters to photometries and from free parameters to derived parameters for different SED models. Generally, the hydrodynamical simulation-based sample is consistent with the empirical statistics-based samples. We attribute the difference between the two set of samples to the different modeling of the SFH, CEH and DAL of galaxies and their different distribution of redshift and physical parameters.

For the Bayesian analysis of the mock data generated in the last section, we employ the same SED modeling procedure and setting of priors for free parameters as in §2.1.1, while the commonly used Gaussian form of likelihood function is employed. The performance of this Bayesian analysis, including its speed and quality, is crucial for the analysis of large sample of galaxies in the big data era. We need some metrics to quantify the performance of parameter estimation which is the main subject of this work.

While the speed of parameter estimation can be easily quantified by the running time, some metrics for the quality of parameter estimation are required. Similar to Acquaviva et al. (2015), we use three metrics to quantify the quality of parameter estimation. Bias, which characterizes the median separation between the predicted and the true values, is defined as:

while the precision, which describes the scatter between predicted and the true values, is defined as:

where $\Delta x=(x_{\mathrm{phot}}-x_{\mathrm{true}})|/(1+x_{\mathrm{true}})$ for redshift (Ilbert et al., 2009; Dahlen et al., 2013; Salvato et al., 2018), and $\Delta x=(x_{\mathrm{phot}}-x_{\mathrm{true}})/(x_{\mathrm{true}}^{\mathrm{max}}-x_{\mathrm{true}}^{\mathrm{min}})$ for other parameters. The median-base definition makes them to be less sensitive to outliers (sources with unexpectedly large errors). The fraction of outliers is defined as:

An important advantage of nested sampling-based algorithm, such as MultiNest, over MCMC-based method is the ability to carry out a simultaneous parameter estimation and model selection. While the main subject of this work is parameter estimation, it is also interesting to explore the effects of sampling parameters (namingly, $nlive$ and $efr$) of MultiNest on the computation of Bayesian evidence, the quantity which is crucial for Bayesian model selection.

In Han & Han (2019), we presented a mathematical framework to discriminate the different assumptions about SSP, DAL and SFH in the SED modeling of galaxies based on the Bayesian evidence for a sample of galaxies. In this work, since the SSP model employed in the generation of mock data is the same as that employed in their Bayesian analysis, we do not need to consider the different choices of SSP. So, the problem is significantly simplified. In this work, we focus on the computation of the Bayesian evidence for the SED modeling of a sample of galaxies with SSP, SFH, and DAL all being assumed to be universal (i.e. ${{\bm{M}}({ssp}_{0},{sfh}_{0},{dal}_{0})}$-like model (See Section 5.1 of Han & Han, 2019).) The sample Bayesian evidence in this case (as Equation 33 of Han & Han (2019)) is:

Although the detailed SFH and DAL of different galaxies can vary significantly, the sample Bayesian evidence computed in this manner remains valuable for identifying the most efficient combination of SFH and DAL for analyzing a vast sample of galaxies, such as the one provided by the CSST wide-field imaging survey.

In practice, the natural logarithmic of Bayesian evidence is commonly used for Bayesian model selection. Therefore, Equation 30 can be rewritten as:

where $\ln(p(\bm{d}_{g}|{{\bm{M}}({ssp}_{0},{sfh}_{0},{dal}_{0})},\bm{I}))$, the Bayesian evidence for an individual galaxy, can be directly obtained in BayeSED with MultiNest. However, the individual Bayesian evidences estimated with MultiNest contain errors. A more strict Bayesian model selection should consider the effects of error propagation. In our case, the error of the sample Bayesian evidence $\ln(BE)$ is simply the sum of errors for individual galaxies which is provided by MultiNest as well.

The minimum $\chi^{2}$ method is also widely used for model selection. For the case with Gaussian observational errors, there is only a constant difference between the minimum $\chi^{2}$ and the natural logarithmic of maximum likelihood. The sample maximum likelihood (as Equation 32 of Han & Han (2019)) is:

Then, the natural logarithmic of sample maximum likelihood is:

where $\ln({\mathop{\rm max}\limits_{\bm{\theta}_{g}}}[p(\bm{d}_{g}|{\bm{\theta}_{g}},{{\bm{M}}({ssp}_{0},{sfh}_{0},{dal}_{0})},\bm{I})])$, the natural logarithmic of maximum likelihood for an individual galaxy, can be directly obtained in BayeSED with MultiNest. Similar to the model selection with Bayesian evidence, only the difference of $\ln(ML)$ between different models is useful for the model selection. Therefore, the model selection with $\ln(ML)$ is equivalent to that with minimum $\chi^{2}$. In §6, we will discuss the difference between the two model selection methods.

As the Bayesian inference engine of BayeSED, MultiNest has some runtime parameters. The values of these runtime parameters have very important effects on the performance of BayeSED for redshift and stellar population parameter estimation of galaxies. Here, we briefly introduce the meaning of these runtime parameters of MultiNest algorithm.

Nested sampling (NS) (Skilling, 2004, 2006), as a Monte Carlo (MC) method primarily designed for the efficient computation of the Bayesian evidence, allows posterior inference as a by-product at the same time. So, it provides a way to carry out simultaneous Bayesian parameter estimation and model selection. As an algorithm built on the NS framework, MultiNest (Feroz & Hobson, 2008; Feroz et al., 2009) is special for its efficiency in sampling from posteriors that may contain several modes and/or degeneracies. It has been improved further by the implementation of importance nested sampling (INS) (Cameron & Pettitt, 2014; Feroz et al., 2019) to increase the efficiency for evidence computation. In the latest version of BayeSED, the V3.12 version of MultiNest, which includes the implementation of INS, is employed.

Similar to most nested sampling algorithms, MultiNest explores the posterior distribution by maintaining a fixed number (See also Higson et al., 2019; Speagle, 2020, for new methods using variable number) of samples drawn from the prior distribution, called live points, and iteratively replaces the point with the lowest likelihood value (the dead point) with another point drawn from the prior but has a higher value of likelihood. While there are many runtime parameters of MultiNest which can be set in BayeSED, only two of them are of particular importance. They largely determined the accuracy and computational cost for the running of MultiNest algorithm and therefore BayeSED. The first one is the total number of live points ($nlive$), which determines the effective sampling resolution. The second one is the target sampling efficiency ($efr$), which determines the ratio of points accepted to those sampled. Generally, the larger $nlive$ and lower $efr$ lead to more accurate posteriors and evidence values but higher computational cost. The optimal value of $nlive$ and $efr$ should be problem-dependent, although $efr$ equals to $0.8$ and $0.3$ are recommended by the authors of MultiNest for parameter estimation and evidence evalutaion, respectively.

In this work, we will explore the effects of $nlive$ and $efr$ on the estimation of photometric redshift and stellar population parameters. The results are presented in §4.1 and §4.2, respectively.

The results of performance tests for photometric redshift estimation are shown in Figure 4. In Figure 4(a), we show the results for only the simplest SED model (SFH=$\tau$,-CEH,DAL=Cal+00) employed in this work. As shown in the top right panel of this figure, in the case without noise, there is a clear anti-correlation between the computation time (or the sampling resolution $nlive$) and the error $\sigma_{\rm NMAD}$ (defined in Equation 28) of photometric redshift estimation. Meanwhile, a smaller $efr$ makes the anti-correlation converge faster with the increasing of $nlive$. There is a clear lower limit for the value of $\sigma_{\rm NMAD}$, which is about $0.006$. As shown in the top left panel of Figure 4(a), in the case with noise, the error of photometric redshift estimation does not always decrease with the sampling resolution $nlive$. When we set $efr=0.1$, the lowest error ($\lesssim 0.056$) of redshift estimation is obtained when $nlive$ is about $25$.When $nlive>25$, the error of redshift estimation start to increase with $nlive$ and finally converge to $\sim 0.058$. This is most likely due to the overfitting to the noise added to the mock data.

The middle two panels of Figure 4(a) show the relation between the computation time (or $nlive$) and the bias (defined in Equation 27) of photometric redshift estimation. In the case with noise, the relation between the computation time (or $nlive$) and bias has almost the opposite profile of that of the error $\sigma_{\rm NMAD}$. However, the bias of photometric redshift estimation is generally very small, which is almost zero in the noise-free case.

The bottom two panels of Figure 4(a) show the relation between the computation time (or $nlive$) and the fraction of outliers OLF (defined in Equation 29) of photometric redshift estimation. Similar to that for $\sigma_{\rm NMAD}$, in the noise-free case, there is a clear anti-correlation between the computation time (or $nlive$) and OLF. In this case, the lower limit for the value of OLF is about $0.002$. In the case with noise, the relation between the computation time (or $nlive$) and OLF has the same profile as that of the error $\sigma_{\rm NMAD}$. When we set $efr=0.1$, the lowest OLF ($\lesssim 0.215$) of redshift estimation is also obtained when $nlive$ is about $25$. When $nlive>25$, the OLF of redshift estimation start to increase with $nlive$ and finally converge to $\sim 0.225$.

In Figure 4(b), we show the results for all of the six SED models with increasing complexity. Here, only the results with $efr=0.1$ are shown. In the case without noise, as shown in the top right panel of this figure, the error $\sigma_{\rm NMAD}$ of photometric redshift estimation tends to converge to a larger value when more complicated SED model is employed. This is not strange, since more complicated SED models have more free parameters and thus suffer from more severe parameter degeneracies. Besides, more complicated SED models apparently require longer running time. The bias of redshift estimation is always very small no matter which SED is employed. In general, when more complicated SED model is employed, the OLF of redshift estimation also increases significantly, and decreases much slower with the increasing of $nlive$.

In the case with noise, as shown in the left panels of Figure 4(b), the results are a little more complicated. For the first three simplest SED models, the error $\sigma_{\rm NMAD}$ of photometric redshift estimation apparently increases with the increasing of model complexity. However, when more complicated forms of SFH is considered, $\sigma_{\rm NMAD}$ start to decreases with the increasing of model complexity, although not very significantly. Meanwhile, the most complicated SED model (SFH=$\alpha$-$\beta$-$\tau$-r,+CEH,DAL=Sal+18) lead to the smallest absolute value of bias, although the bias is actually very small in all cases. The situation for OLF is somewhat similar to that of $\sigma_{\rm NMAD}$. No matter which SED model is employed, when $nlive\gtrsim 25$, both $\sigma_{\rm NMAD}$ and OLF of redshift estimation start to increase, and then slowly decrease to a stable value.

In this subsection, we show the results of performance tests for the photometric stellar population parameter estimation. While the estimates of many stellar population parameter are available, we only show the results for stellar mass and SFR, which are two of the most important physical parameters for the study of the formation and evolution of galaxies.

In this subsection, we present the results of performance test for the computation of Bayesian evidence, a quantity which is crucial for Bayesian model selection.

In Figure 7(a), we show the results for only the simplest SED model (SFH=$\tau$,-CEH,DAL=Cal+00) employed in this work. As shown in the top and middle panels of this figure, the Bayesian evidence computed with importance sampling is more stable than that without importance sampling, especially in the case with noise. So, hereafter and especially in §5 and 6, all mentioned Bayesian evidences are computed with importance sampling. The value of Bayesian evidence increases with the number of live points ($nlive$) which determines the effective sampling resolution. In all cases, it eventually converges to a stable value when $nlive$ is very large, while a smaller sampling efficiency ($efr$) leads to faster convergence rate. A good balance between the speed and quality of Bayesian evidence estimation can be achieved when the MultiNest runtime parameters $efr$ equals to $0.1$ and $nlive$ equals to $50$.

On the other hand, as shown in the bottom panels of Figure 7(a), the error of Bayesian evidence decreases with $nlive$, while a larger sampling efficiency ($efr$) also leads to faster convergence rate. However, unlike the value of Bayesian evidence, the error of Bayesian evidence converges slower with the increasing of $nlive$ in all cases. As a result, if we set $efr=0.1$ and $nlive=50$, the error of Bayesian evidence would be overestimated. A much larger value of $nlive$ seems required to obtain a more reliable estimation for the error of Bayesian evidence with MultiNest, which would be very computationally expansive and not suitable for the analysis of massive photometric data. However, in practice, this may not be a serious issue, since an overestimated error of Bayesian evidence only leads to a more conservative conclusion about model comparison. We just need to keep this in mind.

In Figure 7(b), we show the results for all SED models with increasing complexity, where only the results with $efr=0.1$ are shown. Although the data used for the computation of Bayesian evidence is different for different SED models, the value of Bayesian evidence clearly decreases with the increasing of the complexity of SED model. This is reasonable. Since the same SED model is employed for the generation of mock data and their Bayesian analysis, the mock data can always be interpreted well. However, a more complicated SED model is penalized for being distributed over a larger space, of which only a smaller fraction is useful for the given mock data.

In Figure 8, we investigate the performance of BayeSED combined with the simplest SED model to estimate the photometric redshifts of hydrodynamical simulation-based mock galaxy sample for CSST-like imaging survey. Aa a reference, the panel a of this figure show the ideal case without observational noise and SED modeling errors. So, in this case, the effects of parameter degeneracies, the stochastic nature of MultiNest sampling algorithm and other potential errors in the BayeSED code are the sources of error. As shown clearly, the total error from all of these sources is very small. By comparing panels a and b of this figure, with only the observational noise added, the $\sigma_{\rm NMAD}$ of photometric redshift estimation increases by eight times and the OLF increases by more than forty times. The bias also increases, but is ignorable with respect to $\sigma_{\rm NMAD}$. By comparing panels a and c of this figure, with only the error from the imperfect SED modeling added, the $\sigma_{\rm NMAD}$ of photometric redshift estimation increases by more than three times and the OLF decreases slightly. Besides, there are some additional systematic patterns in the relation between true and estimated values of photometric redshift. The bias also increases and is comparable to $\sigma_{/rmNMAD}$. In general, the observational noise is more important source of error for the photometric redshift estimation of galaxies, although the other one is also very important.

Finally, as shown in panel d of Figure 8, when all sources of error are included, the $\sigma_{\rm NMAD}$ of photometric redshift estimation increases to $0.097$, the OLF increases to $0.264$, and the bias becomes $0.003$. The systematic patterns shown in panel c seems being hidden due to the added noise. The algorithm seems to be struggling to estimate the photometric redshifts correctly by only using the seven-band photometries from CSST-like imaging survey, especially for galaxies with redshift larger than one.

In Figure 9, we investigate the performance of BayeSED combined with the simplest SED model to estimate the stellar population parameters of hydrodynamical simulation-based mock galaxies for CSST-like imaging survey. Aa a reference, the panel a of this figure show the ideal case without observational noise and SED modeling errors. In this ideal case, the $\sigma_{\rm NMAD}$, bias and OLF of stellar mass estimation are $0.047$, $0.005$ and $0.052$, respectively. As shown in panel b, with only the observational noise added, the results become $0.115$, $0.02$ and $0.283$, respectively. As shown in panel c, with only the error from the imperfect SED modeling added, the results become $0.103$, $0.003$ and $0.141$, respectively. By comparing panels a and c, the performance of stellar mass estimation is severely affected by the simplified assumptions in the SED modeling. However, by comparing panels b and c, the observational noise is more important source of error for the photometric stellar mass estimation of galaxies, although the other one is also very important. Finally, as shown in panel d of this figure, when all sources of error are included, the $\sigma_{\rm NMAD}$ of photometric stellar mass estimation increases to $0.135$, the bias increases to $0.034$, and the OLF increases to $0.341$. The algorithm seems to be even more struggling to estimate the photometric stellar mass correctly by only using the seven-band photometries from CSST-like imaging survey.

Similarly, Figure 10 shows the performance of BayeSED combined with the simplest SED model to estimate the SFR of hydrodynamical simulation-based mock galaxy sample for CSST-like imaging survey. Comparing with the results for stellar mass estimation in Figure 9, the photometric SFR estimation is even more severely affected by the simplified assumptions in the SED modeling. Actually, by comparing panels b and c of this figure, it is clear that the error from the imperfect SED modeling is more important source of error for the photometric SFR estimation of galaxies, although the other one is also very important. Finally, it becomes even more struggling to estimate the photometric SFR correctly by only using the seven-band photometries from CSST-like imaging survey.