The formula from Sawicki et al. (2012) used to compute the χ² in the presence...

The formula from Sawicki et al. (2012) used to compute the χ² in the presence of upper limits was not correct. This led the χ² to depend directly on the absolute value of the upper limit. The formula has been rederived and corrected.

CHANGELOG.md

@@ -30,6 +30,7 @@
 - A more explicit message is now given when the flux table cannot be read properly. (Médéric Boquien)
 - Make sure that we do not try to fit data that have an error bar of 0 mJy. (Médéric Boquien)
 - An erroneous warning was displayed when using the `restframe\_parameters` module. (Médéric Boquien)
+- The formula from Sawicki et al. (2012) used to compute the χ² in the presence of upper limits was not correct. This led the χ² to depend directly on the absolute value of the upper limit. The formula has been rederived and corrected. (Médéric Boquien & Denis Burgarella)
 
 ### Optimised
 - The cache architecture has been simplified, making it somewhat faster. It speeds up the model generation by ~1%. (Médéric Boquien)

pcigale/analysis_modules/pdf_analysis/utils.py

@@ -201,9 +201,10 @@ def compute_chi2(model_fluxes, obs_fluxes, obs_errors, lim_flag):
     if limits == True:
         for i, obs_error in enumerate(obs_errors):
             if obs_error < 0.:
-                chi2 += -2. * np.log(np.sqrt(np.pi/2.) * (-obs_errors[i]) * (
-                        1.+erf((obs_fluxes[i] - model_fluxes[i, :]*scaling) /
-                            (np.sqrt(2)*(-obs_errors[i])))))
+                chi2 -= 2. * np.log(.5 *
+                                    (1. + erf(((obs_fluxes[i] -
+                                                model_fluxes[i, :] * scaling) /
+                                               (-np.sqrt(2.)*obs_errors[i])))))
 
     return chi2, scaling