CIGALEMC: GALAXY PARAMETER ESTIMATION USING A MARKOV CHAIN MONTE CARLO APPROACH WITH CIGALE

The spectral energy distribution (SED) of galaxies depends on many physical processes related to the emission from different stellar populations, absorption, and re-emission from dust and gas as well as the possible presence of active galactic nuclei (AGNs). Each process has been studied by many authors; libraries of stellar population models (Fioc & Rocca-Volmerange 1997; Bruzual & Charlot 2003; Maraston 2005), fitting curves for dust emission (Calzetti et al. 1994, 2000; Witt & Gordon 2000), and studies of emission of dust grains (Chary & Elbaz 2001; Dale & Helou 2002; Lagache et al. 2003, 2004; Siebenmorgen & Krügel 2007; Silva et al. 1998; Dopita et al. 2005; da Cunha et al. 2008) are the basis of sophisticated fitting codes which derive physical parameters such as stellar mass, star formation rate (SFR), dust luminosity, and so on.

Many parameters are usually necessary to describe these processes and model theoretical SEDs of galaxies. A grid of theoretical SED models is usually built and fitted to the data and statistical properties are derived for the parameters of interest. A big drawback of any grid-based method is that, for any fitting process, the time to build models grows linearly with the number of models and then about exponentially with the number of parameters involved: such approaches are difficult to implement for complex models involving a sufficiently large number of parameters or when a fine sampling of the parameter space is necessary in order to retrieve statistically robust results. In the past few years, Markov Chain Monte Carlo (MCMC) techniques have started being widely used in physics. In cosmology, parameter estimation from cosmic microwave background data with MCMC methods has been introduced in Christensen et al. (2001) and has been implemented in the publicly available code COSMOMC (Cosmological Monte Carlo; Lewis & Bridle 2002)⁴; in astrophysics, an MCMC approach to the stellar population synthesis modeling has been introduced in Conroy et al. (2009).

Here we use COSMOMC as a generic sampler and interface it to the publicly available code Code Investigating GALaxy Emission⁵ (CIGALE; Noll et al. 2009) in order to allow a fast and accurate evaluation of the multidimensional parameter space probed by this code.⁶ The main advantage of this method is that the computing time to fit the data scales approximately linearly with the number of parameters involved, allowing the user to consider complex models with many parameters for only small additional computational time. MCMC techniques also allow us to probe the shape of the probability distribution, giving far more information than just best fit and marginalized values for the parameters.

The paper is organized as follows; in the next section we briefly describe CIGALE, introducing the main parameters used in the subsequent sections. We then explain the MCMC technique implemented in the modified version of CIGALE, which we call CIGALEMC. We test our code using a mock sample of 62 galaxies already used in Giovannoli et at. (2011) and apply it to a real galaxy sample with data from the Spitzer Infrared Nearby Galaxy Survey (SINGS; see Kennicutt et al. 2003). We always consider a flat cosmological model with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. Finally, we give our results and conclusions.

CIGALE calculates a grid of theoretical SEDs and fits to observational input data constituted by photometric filter fluxes ranging from UV to IR. For a detailed description of the code and its application to real data, we refer the interested reader to a number of other papers (Burgarella et al. 2005; Noll et al. 2009; Giovannoli et al. 2011; Buat et al. 2010, 2011b). In the following, we briefly summarize its main characteristics and the basic parameters used in the next sections. Our notation follows the one introduced in Giovannoli et al. (2011).

2.1. Stellar Populations and Star Formation Rate

CIGALE combines both old and young stellar populations using single stellar populations of Maraston (2005) or PEGASE (Fioc & Rocca-Volmerange 1997). In this paper we will only use Maraston models; we assume star formation histories (SFHs) with either exponentially decreasing SFR in function of time t ("τ models") as

$$\begin{eqnarray} {\rm SFR}_{{\rm old}}(t)&=&{\rm SFR}_{0,{\rm old}}\cdot \,e^{-\frac{(t-t_1)}{\tau _1}}\,\nonumber\\ {\rm SFR}_{{\rm young}}(t)&=&{\rm SFR}_{0,{\rm young}}\cdot \,e^{-\frac{(t-t_2)}{\tau _2}} \end{eqnarray} \tag{ 1 }$$

or "box models" characterized by constant SFR over a limited period of time; in this case the instantaneous SFR at look-back time t' = 0 is given by the galaxy mass divided by the age t of the population, i.e., M_gal/t. Labels 1 and 2 refer to the old and young stellar populations while τ_{1, 2} and t_{1, 2} (both is units of Gyr) are their e-folding time and age, respectively. The two stellar populations are linked and weighted through their mass fraction; the parameter f_ySP represents the fraction of the young stellar mass over the total mass, so that the total instantaneous SFR (output parameter of CIGALE) is expressed as

$$\begin{equation} {{\rm SFR}=(1-f_{{\rm ySP}}){\rm SFR_{old}}(t)+f_{{\rm ySP}}{\rm SFR_{young}}(t)}. \end{equation} \tag{ 2 }$$

2.2. Absorption and Emission by Dust and Gas

In CIGALE, the absorption of star light by dust is described by a Calzetti attenuation curve (Calzetti et al. 1994, 2000); possible modifications of the curve include both the addition of a UV bump at about 2175 Å (Fitzpatrick & Massa 1990, 2007; Noll et al. 2009) and the change of the slope through the multiplication by a power law (λ/λ_V)^δ, where λ_V = 5500 Å is the reference wavelength for the V filter; both the amplitude of the bump and the slope δ are free parameters of the code. The attenuation correction is applied to both stellar populations individually using the visual attenuation parameter of the young stellar populations (A_ySP, in units of magnitudes) and a reduction factor of the attenuation for the old model (f_V) as free parameters.⁷ The code also calculates three additional parameters, derived from the final model SEDs; A_FUV and A_V are defined as the effective attenuation factors in magnitudes at 1500 ± 100 Å and 5500 ± 100 Å while the age t_D4000 is derived from the D4000 break (see Balogh et al. 1999) of the unreddened SED for a single stellar population.⁸

Dust emission in the IR is taken into account using 64 templates of Dale & Helou (2002). These models are parameterized by α, the power-law slope of the dust mass over heating intensity, defined as follows:

$$\begin{equation} {\it dM}_d(U)=U^{-\alpha }{\it dU}, \end{equation} \tag{ 3 }$$

where M_d(U) is the dust mass heated by a radiation field of intensity U.

Bolometric and dust luminosities (L_bol and L_dust, respectively) are derived from all the basic parameters (see Noll et al. 2009) and dust emission due to non-thermal sources such as AGNs can also be added; the fraction f_AGN of dust luminosity L_dust (in L_☉) due to an AGN is estimated using AGN templates from Siebenmorgen et al. (2004a, 2004b).

The spectral line correction due to interstellar gas is performed as in Noll et al. (2009): for the optical band, empirical line templates are taken from the Kinney et al. (1996) starburst spectra while for the UV we use templates derived from SEDs presented in Noll et al. (2004). A correction for the redshift-dependent absorption of the intergalactic medium shortward of the Lyα line is also included using the algorithm of Meiksin (2006).

2.3. Comparison with Data

A grid of theoretical photometric fluxes is calculated at the redshift of the objects considered and a Bayesian analysis is performed through the calculation of the χ² of each model:

$$\begin{equation} \chi ^2({M_{\rm gal}}) = \sum _{i=1}^{k} \frac{(M_{\rm gal} f_{{\rm mod},i} - f_{{\rm obs},i})^2}{\sigma _{{\rm obs},i}^2}\equiv -2 \ln (L), \end{equation} \tag{ 4 }$$

here the galaxy mass M_gal (in M_☉) is treated as a free parameter, f_{mod, i} and f_{obs, i} are the theoretical and experimental fluxes, respectively, the statistical photometry errors are considered in the term σ_{obs, i}, and L is the normalized likelihood function.

In Bayesian inference, the posterior probability of the parameters (${\rm \vec{\theta }}$) of a model in the light of the observed data (${\vec{d}}$) is given by

$$\begin{equation} P(\vec{\theta }|\vec{d})=\frac{P(\vec{d}|\vec{\theta })P(\vec{\theta })}{P(\vec{d})}; \end{equation} \tag{ 5 }$$

here $P{(\vec{d}|\vec{\theta })}\equiv \,L(\vec{\theta })$ is the likelihood of the data given the model, $P{\rm (\vec{\theta })}$ is the prior on the parameters, which quantifies our a priori knowledge of the parameters, and $P{(\vec{d})}$ (called evidence) is a normalization factor. In our case, ${\vec{d}}$ represents the SED of each galaxy while ${\rm \vec{\theta }}$ represents the astrophysical parameters of CIGALE, as {θ}_i ≡ {τ₁, t₂, f_ySP, ...}. An MCMC sampler provides an efficient way to explore the posterior distribution and ensures that the number density of samples is asymptotically proportional to the probability density.

3.1. Metropolis–Hastings Algorithm

The code COSMOMC uses the Metropolis–Hastings algorithm to generate samples; each chain moves according to a transition probability $T(\vec{\theta }_i,\vec{\theta }_{i+1})$ which is determined so that the Markov Chain has a stationary asymptotic distribution equal to the posterior distribution $P(\vec{\theta })$ that we want to sample from. Given an arbitrary proposal density distribution ${q(\vec{\theta }_i,\vec{\theta }_{i+1})}$ to propose a new point $\vec{\theta }_{i+1}$ when the chain is at the point ${\vec{\theta }_{i}}$, the probability of transition is given by β:

$$\begin{equation} \beta (\vec{\theta }_i,\vec{\theta }_{i+1})\,=\,\min \left \lbrace 1,\frac{P(\vec{\theta }_{i+1})q(\vec{\theta }_{i+1},\vec{\theta }_i)}{P(\vec{\theta }_{i})q(\vec{\theta }_i,\vec{\theta }_{i+1})}\right \rbrace \end{equation} \tag{ 6 }$$

so that

$$\begin{equation} T(\vec{\theta }_i,\vec{\theta }_{i+1})=\beta (\vec{\theta }_i,\vec{\theta }_{i+1})q(\vec{\theta }_i,\vec{\theta }_{i+1}). \end{equation} \tag{ 7 }$$

This ensures that the detailed balance holds

$$\begin{equation} P(\vec{\theta }_{i+1})T(\vec{\theta }_{i+1},\vec{\theta }_i)=P(\vec{\theta }_i)T(\vec{\theta }_i,\vec{\theta }_{i+1}) \end{equation} \tag{ 8 }$$

and that the distribution converges to $P{\rm (\vec{\theta })}$. In practice, a random number x ∈ [0: 1] is generated in the process of moving from $\vec{\theta }_i$ to $\vec{\theta }_{i+1}$ so that the new point $\vec{\theta }_{i+1}$ is accepted if β ⩾ x. This ensures that each point of the chain depends only on its predecessor; in this sense the chain is a Monte Carlo Markov process.

3.2. Comparison with Grid-based Methods, Burn in, and Convergence Diagnostics

As an illustration of the sampling mechanism, in Figure 1 we plot samples from the posterior distribution for an MCMC run with a test galaxy taken from a mock sample at redshift z ∼ 0.7 (see the following section for details); the number density of samples in the plane is proportional to the probability density of these two parameters. The dust luminosity L_dust strongly depends on the SFR and the two parameters are degenerate, as shown by the colors in the figure. This plot clearly shows the efficiency of this MCMC method. In the grid-based approach, the parameter space is sampled in the same "blind" way for high and low values of the posterior: this can be an issue for both reliability of results and computation time, as also pointed out in Noll et al. (2009) and Acquaviva et al. (2011). In fact, local minima and degeneracies between parameters can be easily missed or undersampled without a good a priori knowledge of the parameter space; the oversampling of an ill-constrained parameter can also lead to a slight degradation of the estimates of well-constrained parameters and many points can be generated in a region where the posterior is low, resulting in a waste of CPU time. This is not the case when MCMC chains are used because each chain "learns" where to move in the parameter space through the Metropolis–Hastings algorithm so that the density of samples is proportional to the posterior distribution. Degeneracies between parameters are more easily found, especially if many chains, starting from different regions in the parameter space, are used. In other words, the CIGALEMC user needs to specify the prior parameter space (number of parameters and their limits) but not the density of points for each parameter. In the following section, we will provide a comparison of CPU time between the original CIGALE and CIGALEMC when evaluating physical properties of a mock sample. The code also calculates the covariance between various parameters so that an initial run can be made and the covariance matrix obtained can be used to improve the efficiency of sampling for subsequent runs.

Since each MCMC chain starts at a random position in the parameter space, it will take a little time before the chain equilibrates and starts sampling the posterior distribution. This period of initial convergence is called burn in period and the first burn in points of each chains will be discarded when doing any statistical analysis. In order to obtain uncorrelated samples of the posterior each chain is also "thinned" by using only occasional points of it; the thinning factor varies according to the number of parameters involved and it is typically in the range 25–50. The code allows to choose the burn in fraction of the chain we want to discard and automatically thins out the chains.

In our analyses, we would not use any initial covariance matrix for the parameters and, in order to be conservative, we compute statistical quantities using only the second-half of each chain. In Figure 2, we plot the points of an MCMC chain for a galaxy sample in the plane L_dust versus f_AGN; the chain reaches the sensitive region of the parameter space after only a few "burn in" points characterized by very low values of L_dust. Having a set of samples from the full posterior distribution, it is possible to calculate statistical quantities for the parameters of interest. Since the number density of sampling points is proportional to the posterior density, it is easy to calculate mean values and marginalized one-dimensional distribution for each parameter θⁱ by simply counting the number N of samples within binned ranges of parameter values:

$$\begin{equation} \langle \,\theta ^i\rangle \sim \frac{1}{{N}}\sum _{j=1}^N\theta ^i_j, \end{equation} \tag{ 9 }$$

while this is much more difficult in the context of the numerical grid integration because the calculation time grows exponentially with the number of dimensions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to be sure that MCMC chains are efficiently sampling the posterior distribution (and then obtain robust statistics for each parameter), it is important to check their convergence. The code COSMOMC provides two convergence criteria for runs with one single chain (Raftery & Lewis 1992) and with multiple chains (Gelman & Rubin 1992). In the following analysis, we will run multiple chains using the Gelman & Rubin diagnostic, which is characterized by the "variance of chain means"/"mean of chain variances" parameter R; |R − 1| ⩽ 0.03 is usually enough to reach convergence and stop the chains.

In this work we use COSMOMC as a generic MCMC sampler and link it to CIGALE in order to allow a faster exploration of the astrophysical parameter space. Our modified CIGALEMC code will be publicly available very soon⁹ and, since it is based on COSMOMC for the sampling options, convergence criteria and statistical quantities provided, we refer the reader to the Web site¹⁰ and to Lewis & Bridle (2002) and references therein for a detailed explanation of the code and MCMC methods in general.

We test our CIGALEMC code with a mock sample already used in Giovannoli et al. (2011). We consider 62 artificial galaxy SEDs corresponding to Luminous InfraRed Galaxies at redshift z ∼ 0.7 and obtained by varying the input parameters of CIGALE in the following ranges: [0.1 Gyr ⩽ τ₁ ⩽ 10 Gyr, 0.025 Gyr ⩽ t₂ ⩽ 0.7 Gyr, 0 ⩽ f_ySP ⩽ 1, 0.6 mag ⩽ A_ySP ⩽ 2.1 mag, 0 ⩽ f_V ⩽ 1, 1.0 ⩽ α ⩽ 2.5, 0 ⩽ f_AGN ⩽ 0.3]; such a wide multidimensional parameter space allows us to model very different SEDs of galaxies characterized by a wide range of possible SFHs, absorption, and emission by gas and dust and AGN contamination. All galaxies are based on a Salpeter initial mass function; the age of the old stellar population is fixed at 7 Gyr; we consider a constant SFR for the young population model (τ₂ = 20 Gyr) and do not add any modification to the original Calzetti attenuation curve (no UV bump and δ = 0).

Theoretical fluxes are calculated in the following 17 bands from UV to IR: 0.231 μm for the Galaxy Evolution Explorer (GALEX), [0.35–0.36–0.46–0.54–0.65–0.87–0.90–1.2–1.6–2.1] μm (corresponding to MUSYC bands), [3.6–4.5–5.8–8.0] μm (corresponding to IRAC photometry), and [24–70] μm (corresponding to MIPS photometry). We add a Gaussian distributed uncertainty σ to each theoretical flux; its value is 10% of the corresponding flux, a reasonable choice for the measurements considered above. For each galaxy we run eight chains with initial positions randomly chosen in the parameter space determined by the following set of input parameters:

$$\begin{equation} \lbrace {\tau _1,t_2,f_{{\rm ySP}},A_{{\rm ySP}},f_{V},\alpha,f_{{\rm AGN}},M_{{\rm gal}}}\rbrace, \end{equation} \tag{ 10 }$$

from which statistical quantities of interest can be calculated for the following set of derived parameters:

$$\begin{equation} \lbrace {{\rm SFR},t_{{\rm D4000}},L_{{\rm bol}},L_{{\rm dust}},A_{{\rm FUV}},A_{V},M_{{\rm star}}}\rbrace. \end{equation} \tag{ 11 }$$

We assume flat priors in the following parameters space: [0.1 Gyr ⩽ τ₁ ⩽ 10 Gyr, 0.025 Gyr ⩽ t₂ ⩽ 2 Gyr, 0 ⩽ f_ySP ⩽ 1, 0 mag ⩽ A_ySP ⩽ 3 mag, 0 ⩽ f_V ⩽ 1, 0.5 ⩽ α ⩽ 3, 0 ⩽ f_AGN ⩽ 1, 8 M_☉ < M_gal < 14 M_☉]; chains are stopped when the Gelman & Rubin R-1 parameter is |R − 1| ∼ 0.03.

First of all, we want to check that our results are statistically in agreement with the input values for the mock sample. As a tool derived from COSMOMC, CIGALEMC allows to calculate and plot the mean likelihood and marginalized distribution for each parameter.

The marginalized distribution in a given direction of the parameter space $\vec{d}={\bf h}(\vec{\theta })$ (where ${\bf h}(\vec{\theta })$ is the projector operator in one of the parameters considered, as ${\bf h}(\vec{\theta })=\theta _i$) is proportional to the number of samples at $\vec{d}$ and it can be expressed as

$$\begin{equation} P(\vec{v}) = M(P,\vec{v}) \equiv \int d\vec{\theta } \, P(\vec{\theta }) \delta ({\bf h}(\vec{\theta }) - \vec{v}), \end{equation} \tag{ 12 }$$

where $P(\vec{\theta })$ is the posterior distribution. Assuming flat priors on $\vec{\theta }$, the mean likelihood of samples with $\vec{d}={\bf h}(\vec{\theta })$ can be expressed as

$$\begin{equation} \langle \,P(\vec{v}:{\bf h}(\vec{\theta })=\vec{v})\rangle \equiv \frac{\int d\vec{\theta }\, P(\vec{\theta })^2 \delta ({\bf h}(\vec{\theta }) - \vec{v})}{\int d\vec{\theta } \,P(\vec{\theta }) \delta ({\bf h}(\vec{\theta }) - \vec{v})} = \frac{M(P^2,\vec{v})}{M(P,\vec{v})}. \end{equation} \tag{ 13 }$$

If $P(\vec{\theta })$ is a multivariate Gaussian distribution, it is possible to demonstrate that both mean likelihood and marginalized distribution are Gaussian and proportional so that they look the same: differences in these distributions will be then a signal of non-Gaussianity, which can be either intrinsic or due to parameters not very well constrained. In Figure 3 we show, for an artificial galaxy of the mock sample considered, both marginalized distributions (black solid lines) and mean likelihood (dotted lines) for some parameters of interest: we can see that τ₁, t₂, f_ySPf_V, and α are not very well constrained by the code. Similar results have been found in the analyses by Noll et al. (2009), Giovannoli et al. (2011), and Buat et al. (2011b).

In order to study the goodness of the fit in a quantitative way for the whole sample of galaxies, we introduce the quantity

$$\begin{equation} {X^{i}} \equiv \,\frac{1}{{N}} \sum _{j}\frac{{O^i_{j}-I^i_{j}}}{\sigma ^i_{j}}; \end{equation} \tag{ 14 }$$

here i runs for all the parameters considered by the code (so that {X}ⁱ = {τ₁, t₂, f_ySP, ...}), j runs for the N = 62 objects of the mock, Oⁱ is the set of the best-fit values obtained as output from CIGALEMC, Iⁱ is the vector of input parameters, and σⁱ is the vector of the 68% confidence level marginalized uncertainties. As we can see from Figure 4, all X values are compatible with X = 0, which means that the code is able to find the best-fit values of the parameters with high confidence. However, some distributions are slightly skewed: this can be due to chains being stuck in local minima of the likelihood function so that, in some cases, the best fit found does not correspond to the input value. This situation is typical when parameters are not well constrained (as for t₂) or in the presence of strong degeneracies between parameters. As an example, Figure 5 shows the 68% and 95% confidence levels in the plane f_ySP versus t₂ for one sample galaxy of the mock: as we can see, limits on these parameters are not very strong; we also note a partial degeneracy for high values of f_ySP, which is due to the fact that these two parameters affect the galactic SFR in the same way since we consider constant SFR for the young model:

$$\begin{equation} {\rm SFR}\sim {f_{{\rm ySP}}\frac{M_{{\rm gal}}}{t_2}}. \end{equation} \tag{ 15 }$$

In order to check the reliability of our MCMC algorithm, we tested the parameter uncertainty estimation. We created and fit the SEDs of 100 artificial galaxies in 21 bands from IR to UV, built with the same input parameter set and with a 10% scatter in flux, corresponding to the photometric error considered. We checked that in general we are able to find the true values within the region allowed at 68% (95%) confidence 68% (95%) of the time within the Poisson fluctuation error for different choices for both the input parameters and the prior ranges allowed in the fit.

Finally, we checked the consistency of our results by calculating the Pearson correlation coefficient r between the input values of the parameters used to generate the mock and the best-fit values found with CIGALEMC. The Pearson correlation coefficient quantifies the amount of correlation between two variables X and Y and it assumes values in the range r ∈ [ − 1: 1] (see Cohen 1988). For samples X_i and Y_i, it can be written as

$$\begin{equation} r=\frac{\sum _{i=1}^{N}(X_i-\bar{X})(Y_i-\bar{Y})}{(\sum _{i=1}^{N}(X_i-\bar{X})^2)^{{1}/{2}}(\sum _{i=1}^{N}(Y_i-\bar{Y})^2)^{{1}/{2}}}, \end{equation} \tag{ 16 }$$

here $\bar{X}$ and $\bar{Y}$ denote the mean values of X_i and Y_i. The interpretation of the strength of the correlation can depend on the context; however, a widely used standard introduced by Cohen (1988) considers Pearson values |r| > 0.5 as a signal of large correlation between variables, while |r| < 0.3 denotes poor correlation. Again, as we can see from Table 1, some parameters (τ₁, t₂, f_ySP) have small values of r; this is expected since these parameters are mostly unconstrained; however, even if the exact values of the input parameters are not found, they always fall inside the statistically significant region of the parameter space. We also explicitly checked that results do not significantly change when increasing the prior space for some parameters not very well constrained as τ₁, t₂, and α. Finally, we tested, for a galaxy of the mock, how results change when considering a Gaussian prior, with different central values and width σ = 0.2, on one of the most unconstrained parameters, t₂; we found that results on the other parameters are always statistically consistent with respect to the choice of a flat prior on t₂.

In general, it is a good idea to choose the largest possible prior space for unconstrained parameters in order to avoid possible biases due to the prior choice.

It is useful to compare the performance of CIGALE and CIGALEMC in terms of computing time, especially since computation can become prohibitive for any grid-based method if the number of parameters involved is sufficiently high. The CPU time required to obtain convergence of the chains for each galaxy mainly depends on both the quality of the data and the number of parameters considered. Running eight chains in parallel (each one on a 2.40 GHz Intel Xeon E5530), for each galaxy of the mock we typically reach a good convergence with ∼20,000 points for each chain, which means a total of ∼ 160,000 points. The grid built in Giovannoli et al. (2011) to analyze the same sample with the same number of free parameters and bounds contained ∼3.5 × 10⁶ points; this means a gain of order ∼20 in efficiency in the estimation of the parameters but a more dramatic efficiency can be easily reached when we need to use either more parameters or a fine sampling in a given direction of the parameter space or both. The average CPU time to reach a good convergence for a galaxy of this mock is about 35 s, which translates in 280 s of total CPU time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now use CIGALEMC to infer physical properties of the well-known SINGS (Kennicutt et al. 2003) sample. In order to make a comparison with results obtained using the grid-based CIGALE, we use the same 39 galaxies already considered in Noll et al. (2009) with the same spectral coverage: GALEX far-UV (∼1500 Å) and near-UV (∼2300 Å) filters (Gil de Paz et al. 2007), Two Micron All Sky Survey (2MASS) data for J, H, and K_s (Jarrett et al. 2003), IRAC, and MIPS filters for [3.6, 4.5, 5.8, 8.0, 24, 70, 160] μm (Dale et al. 2005, 2007, 2008) optical data for B, V, R, and I bands corrected as in Muñoz-Mateos et al. (2009) and fluxes from u', g', r', i', and z' filters of the Sloan Digital Sky Survey (SDSS; Stoughton et al. 2002); Dale et al. (2007, 2008) optical data are only used for 14 galaxies for which SDSS data are not available. The mean photometric relative uncertainties for the bands considered are shown in Table 2; in particular, very small uncertainties are associated with the 2MASS bands. We performed a preliminary run, realizing that the hard constraints coming from this filter set did not allow the code to properly fit for the other filters. Since these uncertainties do not take into account for possible calibration errors and since a systematic offset can also affect the CIGALE theoretical models, we decided to be conservative and, following Noll et al. (2009), we add 5% uncertainties in quadrature for each filter set.

In our analysis, we assume the following range of variation for a set of nine astrophysical parameters: [0.1 Gyr ⩽ τ₁ ⩽ 10 Gyr, 0.025 Gyr ⩽ t₂ ⩽ 2 Gyr, 0 ⩽ f_ySP ⩽ 1, −0.5 ⩽ δ ⩽ 0.5, 0 mag ⩽ A_ySP ⩽ 5 mag, 0 ⩽ f_V ⩽ 1, 0.5 ⩽ α ⩽ 3, 0 ⩽ f_AGN ⩽ 1, 8 M_☉, M_gal < 13 M_☉]. We keep fixed both the age of the old stellar population (t₁ = 10 Gyr) and the e-folding time for the young stellar population (τ₂ = 20 Gyr): these parameters are not well constrained by the data so that fixing them does not alter the fit. Finally, we only consider models with a Salpeter initial mass function and solar metallicity; metallicity measurements are quite uncertain due to many limiting factors (Noll et al. 2009; Moustakas & Kennicutt 2006, and references therein); Noll et al. (2009) checked the influence of different assumptions for the metallicity, concluding that deviations in the properties of the galaxies are within the uncertainties, which tend to increase by 0%–20% when half of or double the solar metallicity are considered. The only exception is the absolute value of t_D4000 because of the well-known age–metallicity degeneracy (e.g., Kodama & Arimoto 1997).

Our reference AGN model has L = 10¹² L_☉, R = 125 pc, and A_V = 32 for the luminosity of the non-thermal source, the outer radius of a spherical dust cloud covering the AGN, and the amount of attenuation in the optical caused by the cloud, respectively.

Our findings can be summarized as follows.

• 1.  
  Very good constraints are derived for the AGN fraction of all sources (f_AGN ⩽ 0.10) except for NGC 0584 and NGC 1404 for which f_AGN ⩽ 0.6 at the 68% confidence level. We note that the flux at 160 μm for NGC 0584 can be contaminated by some foreground/background emission and the most recent Herschel data conclusively confirm that a background source contaminates both fluxes at 70 and 160 μm for NGC 1404 (D. A. Dale 2011, private communication); the low quality of these data is responsible for the big uncertainties obtained for other parameters such as SFR, L_bol, L_dust, and A_FUV; in Figure 6 we plot the two-dimensional marginalized distribution for SFR versus f_AGN for NGC 1404: the double-peaked likelihood function is most probably an artifact due to the low quality of data for this source. In general, we note that that NGC 0584 and NGC 1404 are elliptical galaxies which tend to have very weak but warm dust emission; it is not surprising that they show an apparently high AGN fraction. We decided not to use these sources in the rest of our analysis.
• 2.  
  We are not able to put strong constraints on "phenomenological" parameters as τ₁, t₂, and δ; limits on these parameters depend on the assumed prior range. The poor determination of δ is essentially due to the low number of data in UV. In general, both δ and the possible UV bump allowed by the code are very difficult to constrain with few only broadband data (Buat et al. 2011a, 2011b). The fraction of the young stellar mass over the total mass is well constrained, with values f_ySP < 0.26 for all the galaxies, except for NGC 1705, NGC 2798, and NGC 4631, for which f_ySP is unconstrained. The mean value for the parameter α of the Dale & Helou (2002) templates is α = 2.44 (in agreement with Noll et al. 2009; we note that, while low values of α are well constrained, no constraints are found on high values, with α = 3 compatible at the 68% confidence level for all galaxies expect NGC 2798, for which α < 1.77 at the 68% confidence level; weak constraints on high values for α are due to the degeneracy of the dust emission models for wavelengths short of the emission peak, so that, for α > 2.5, the flux ratio f_ν(60) μm/f_ν(100) μm is almost constant (see Noll et al. 2009).
• 3.  
  A frequency-dependent χ² analysis shows which bands mostly contribute to the total χ² for each galaxy. We introduce the averaged χ²:

  $$\begin{equation} {\langle \chi ^2(\lambda ^j)\rangle =\frac{1}{N_j}\sum _{i=1}^{N_j}\left (\frac{{\rm Th}_i^j-d_i^j}{\sigma ^j_i}\right)^2}; \end{equation} \tag{ 17 }$$

  here Th^j_i and d^j_i are, respectively, the theoretical best-fit SEDs and the data points for the ith galaxies at the jth frequency and the sum runs over our sample of galaxies. From Figure 7, we see that the code is able to find a good agreement with the data for all the bands considered. The frequency-dependent χ² is particularly low for bands with large relative errors as UV GALEX bands and optical B, V, R, and I bands. The large uncertainty found for the MIPS 70 μm filter is mainly due to the galaxies NGC 1715 and NGC 5866; the code is not able to find a good fit for the 70 μm filter and the reduced χ² is high (χ²/degrees of freedom (dof) = 7.8 for both galaxies, see Table 3). A previous analysis by Cannon et al. (2006), based on Spitzer observations, found similar results for NGC 1705, showing in particular that the models of Li & Draine (2001, 2002) give a better fit to IR data than the Dale & Helou (2002) models used in this paper (see Figure 3 of Cannon et al. 2006). Interestingly enough, NGC 1715 and NGC 5866 are, respectively, the only irregular and S0 galaxies in our sample.In general, we find a small trend of worse fitting for galaxies with small SFR values. Galaxies with the lowest values of the χ² do not show peculiar properties, with total stellar masses of about 10¹¹ M_☉ and SFRs between 0.1 and 10 M_☉ yr⁻¹, typical of nearby spiral galaxies; an exception is the dwarf galaxy NGC 4625 (reduced χ²/dof = 1.2), characterized by smaller values for both the stellar mass and luminosity. Three other dwarf galaxies (NGC 1705, NGC 2976, and NGC 5474) are clearly identified in our sample by looking at their smaller values for both the luminosity and the stellar masses with respect to the rest of the sample. Finally, in order to check for a possible correlation between the filter set used and the goodness of the fit, we also calculated the mean $\bar{\chi ^2}$ for the subset of galaxies for which SDSS filters are available with respect to the subset of galaxies with optical Dale et al. filters; we found no correlation, with $\bar{\chi }\sim 2.5$ in both cases.
• 4.  
  Good constraints can be derived for the mass-dependent parameters M_star, SFR, L_bol, and L_dust as shown in Table 3. The Pearson correlation coefficient with results from Noll et al. (2009) is r = 0.99, 0.98, 0.94, 1.0, 1.0, and 0.94 for ${M_{{\rm star}},\,{\rm SFR},\,{\rm t_{{\rm D4000}}},\,L_{{\rm bol}},\,L_{{\rm dust}}, {\rm and}\,A_{{\rm FUV}}}$, respectively. In Figure 8, we show a quantitative comparison with the analysis performed by Noll et al. (2009) with the original code CIGALE by plotting the ratio, $Q\equiv (\mathrm{\texttt {CIGALEMC}-\mathrm{\texttt{CIGALE}}})/\sigma _{\texttt {CIGALEMC}}$, where CIGALEMC and CIGALE refer to the mean value of the parameters quoted in Table 3 for this work and Noll et al. (2009), respectively. Possible systematic differences between the results of both methods can be studied by considering mean, standard deviation, and skewness for the parameters of interest. As we can see from Table 4, no significant difference between results in this paper and in Noll et al. (2009) is found for SFR, log (t_D4000), log (L_bol), and A_FUV, with each value compatible with 0 at less or about the 68% confidence level; however, there is some skewness associated with SFR, t_D4000, and L_bol. In the cases of M_star and L_dust, the mean values are lower from 0 at about the 95% confidence level, indicating the existence of some offset between results in this paper and in Noll et al. (2009). A further comparison of the uncertainties between our results and Noll et al. (2009) shows that we are able to put stronger constraints on the galaxy SFR, while our constraints on M_star, L_bol, L_dust, andA_FUV are weaker. Differences between our results and results in Noll et al. (2009) can be due to the different choice of the input parameter space in terms of both number and type of parameters used in the analysis; in particular, we consider a wider range of variation for α, δ, A_ySP, t₂, and f_AGN, while we decide to fix the parameter τ₂ (see Table 3 of Noll et al. 2009 for their choice of the parameter space); however, it is reassuring to see how different theoretical assumptions lead to compatible results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Values of the First Three Moments of the $Q\equiv (\texttt {CIGALEMC}-{\texttt{CIGALE}})/\sigma _{\texttt {CIGALEMC}}$ Distributions for Some Parameters of Interest

Parameters	Mean	$\sqrt{\mathrm Var}$	Skewness
log10 Mstar	−1.0	0.52	−0.1
log10SFR	0.29	0.89	1.79
logtD4000	1.2	2.8	2.6
log10Lbol	−0.45	0.36	1.17
log10Ldust	−0.58	0.25	0.01
AFUV	−0.69	0.62	−0.25

Download table as:  ASCIITypeset image

In this paper, we have introduced an MCMC sampling method for the astrophysical parameter estimation from SED fitting with CIGALE. We have shown the following advantages of our modified CIGALEMC code over the usual grid-based CIGALE:

• 1.  
  Its efficiency, in terms of CPU time, through the Metropolis–Hastings algorithm. Most of the sampling points are drawn in the region where the posterior probability is high, while in the grid-based approach all regions are sampled in the same way. Moreover, marginalized one-dimensional probability distributions for the parameters of interest are calculated by simply counting the number of samples within a binned range of parameter values, the density of sampling points being proportional to the posterior probability. It is hard to do the same with the usual grid approach, since the integration calculation scales exponentially with the number of dimensions. The analysis of a mock sample shows that CIGALEMC needs 20 times less points than CIGALE to reach convergence for a given galaxy but, in general, results depend on both the number of parameters and the prior used.
• 2.  
  Its accuracy; degeneracies between parameters are easily found, convergence criteria (already implemented in the code) ensure that statistical quantities for the parameters of interest are robustly determined and cross-checks through statistical analysis of mock catalogs are not necessary.
• 3.  
  Its "user friendly" characteristics. The user does not need to decide a priori the number density of samples for each region, trying to find a compromise between the accuracy of the results and the speed of the code: only the prior range must be chosen in advance for CIGALEMC; in fact the Metropolis–Hastings algorithm automatically samples adequately the posterior probability according to its values.

Our code will be available very soon at this Web site: http://www.oamp.fr/cigale/.

P.S. thanks Denny Dale for very useful discussions. We thank Antony Lewis for providing a publicly available version of COSMOMC.