Estimate the fluxes in individual bands through a Bayesian analysis.

CHANGELOG.md

@@ -6,6 +6,7 @@
 - It is now possible to fit any physical property indicated by the code (e.g. equivalent width, dust luminosity, etc.). For this the physical property needs to be given in the input file and the properties to be fitted must be given in the properties filed in pcigale.ini. (Héctor Salas & Médéric Boquien)
 - It is now possible to fit emission lines. For this the line has to be indicated in the same way as any other band both in the input flux file (in units of W/m²) and in the list of bands in `pcigale.ini`. Lines are prefixed with `line.` followed by the name of the line, for instance `line.H-alpha` for Hɑ. The following lines are supported at the moment: `Ly-alpha`, `CII-133.5`, `SiIV-139.7`, `CIV-154.9`, `HeII-164.0`, `OIII-166.5`, `CIII-190.9`, `CII-232.6`, `MgII-279.8`, `OII-372.7`, `H-10`, `H-9`, `NeIII-386.9` `HeI-388.9`, `H-epsilon`, `SII-407.0`, `H-delta`, `H-gamma`, `H-beta`, `OIII-495.9`, `OIII-500.7`, `OI-630.0`, `NII-654.8`, `H-alpha`, `NII-658.4`, `SII-671.6`, `SII-673.1`. (Médéric Boquien)
 - Two new dust attenuation modules have been added: `dustatt\_modified\_CF00` and `dustatt\_modified\_starburst`. The former implements a modified 2-component Charlot & Fall (2000) model whereas the latter implements a modified starburst law with the continuum attenuated with a Calzetti (2000) curve and the lines extincted with a Milky Way or a Magellanic Cloud law. The previous models `dustatt\_powerlaw`, `dustatt\_2powerlaws`, and `dustatt\_calzleit` are still available but are deprecated. (Médéric Boquien & David Corre)
+- In addition to the physical properties, the fluxes are now also estimated through a Bayesian analysis. (Médéric Boquien)
 
 ### Changed
 - The `sfhdelayed` module has been extended to optionally include an exponential burst to model the latest episode of star formation. (Médéric Boquien & Barbara Lo Faro)

pcigale/analysis_modules/pdf_analysis/workers.py

@@ -192,6 +192,15 @@ def analysis(idx, obs):
             if gbl_models.conf['analysis_params']['save_chi2'] is True:
                 save_chi2(obs, prop, gbl_models, chi2, values)
 
+        for band in gbl_results.bayes.fluxmean:
+            values = gbl_models.flux[band][wz]
+            mean, std = weighted_param(values[wlikely] * scaling_l,
+                                       likelihood)
+            gbl_results.bayes.fluxmean[band][idx] = mean
+            gbl_results.bayes.fluxerror[band][idx] = std
+            if gbl_models.conf['analysis_params']['save_chi2'] is True:
+                save_chi2(obs, band, gbl_models, chi2, values * scaling)
+
         best_idx_z = np.nanargmin(chi2)
         gbl_results.best.chi2[idx] = chi2[best_idx_z]
         gbl_results.best.scaling[idx] = scaling[best_idx_z]

pcigale/managers/results.py

@@ -47,6 +47,8 @@ class BayesResultsManager(object):
         self.interror = {prop: SharedArray(nobs) for prop in intpropnames}
         self.extmean = {prop: SharedArray(nobs) for prop in extpropnames}
         self.exterror = {prop: SharedArray(nobs) for prop in extpropnames}
+        self.fluxmean = {band: SharedArray(nobs) for band in models.flux}
+        self.fluxerror = {band: SharedArray(nobs) for band in models.flux}
         self.weight = SharedArray(nobs)
 
     @property
@@ -102,6 +104,12 @@ class BayesResultsManager(object):
         exterror = {prop: np.array([result.exterror[prop]
                                     for result in results])
                     for prop in merged.exterror}
+        fluxmean = {band: np.array([result.fluxmean[band]
+                                   for result in results])
+                   for band in merged.fluxmean}
+        fluxerror = {band: np.array([result.fluxerror[band]
+                                    for result in results])
+                    for band in merged.fluxerror}
         weight = np.array([result.weight for result in results])
 
         totweight = np.sum(weight, axis=0)
@@ -139,6 +147,19 @@ class BayesResultsManager(object):
                     np.maximum(0.05 * merged.extmean[prop],
                                merged.exterror[prop])
 
+        for band in merged.fluxmean:
+            merged.fluxmean[band][:] = np.sum(
+                fluxmean[band] * weight, axis=0) / totweight
+
+            # We compute the merged standard deviation by combining the
+            # standard deviations for each block. See
+            # http://stats.stackexchange.com/a/10445 where the number of
+            # datapoints has been substituted with the weights. In short we
+            # exploit the fact that Var(X) = E(Var(X)) + Var(E(X)).
+            merged.fluxerror[band][:] = np.sqrt(np.sum(
+                weight * (fluxerror[band]**2. + (fluxmean[band]-merged.fluxmean[band])**2), axis=0) / totweight)
+
+
         return merged
 
 
@@ -356,6 +377,11 @@ class ResultsManager(object):
                                     name="bayes."+prop))
             table.add_column(Column(self.bayes.exterror[prop],
                                     name="bayes."+prop+"_err"))
+        for band in sorted(self.bayes.fluxmean):
+            table.add_column(Column(self.bayes.fluxmean[band],
+                                    name="bayes."+band))
+            table.add_column(Column(self.bayes.fluxerror[band],
+                                    name="bayes."+band+"_err"))
 
         table.add_column(Column(self.best.chi2, name="best.chi_square"))
         obs = [self.obs.table[obs].data for obs in self.obs.tofit]