diff --git a/pcigale/analysis_modules/pdf_analysis/__init__.py b/pcigale/analysis_modules/pdf_analysis/__init__.py
index 6fd74b7e954eac0240524875d67f0fb2588dd9b4..b7a46af63fe15f25b96c0f6fd2f20cf042ec0cce 100644
--- a/pcigale/analysis_modules/pdf_analysis/__init__.py
+++ b/pcigale/analysis_modules/pdf_analysis/__init__.py
@@ -195,6 +195,7 @@ class PdfAnalysis(AnalysisModule):
         obs.save('observations')
 
         results = self._compute(conf, obs, params)
+        print("\nSanity check of the analysis results...")
         results.best.analyse_chi2()
 
         print("\nSaving the analysis results...")
diff --git a/pcigale/analysis_modules/pdf_analysis/workers.py b/pcigale/analysis_modules/pdf_analysis/workers.py
index d8e5f474322056e78cf3f73350c506d1d7d2e07a..a47f897102a6033770f501939b33bd3d251f8b11 100644
--- a/pcigale/analysis_modules/pdf_analysis/workers.py
+++ b/pcigale/analysis_modules/pdf_analysis/workers.py
@@ -207,11 +207,6 @@ def analysis(idx, obs):
         gbl_results.best.scaling[idx] = scaling[best_idx_z]
         gbl_results.best.index[idx] = (wz.start + best_idx_z*wz.step +
                                        gbl_models.block.start)
-    else:
-        # It sometimes happens because models are older than the Universe's age
-        print(f"No suitable model found for the object {obs.id}. It may be that"
-              f" models are older than the Universe or that your chi² are very "
-              f"large.")
 
     gbl_counter.inc()
 
diff --git a/pcigale/managers/results.py b/pcigale/managers/results.py
index 9abd36bafb05436496d656eb7c0e70d5d097e98d..190f45dec0ccacb685bcc9137b072a616c6e954e 100644
--- a/pcigale/managers/results.py
+++ b/pcigale/managers/results.py
@@ -308,15 +308,23 @@ class BestResultsManager(object):
          objects seems to be overconstrainted.
 
         """
+        # If no best model has been found, it means none could be properly
+        # fitted. We warn the user in that case
+        bad = self.obs.table['id'][np.isnan(self.chi2)].tolist()
+        if len(bad) > 0:
+            print(f"No suitable model found for {', '.join(bad)}. It may be "
+                  f"that models are older than the universe or that your χ² are"
+                  f" very large.")
+
         obs = [self.obs.table[obs].data for obs in self.obs.tofit]
         nobs = np.count_nonzero(np.isfinite(obs), axis=0)
         chi2_red = self.chi2 / (nobs - 1)
         # If low values of reduced chi^2, it means that the data are overfitted
         # Errors might be under-estimated or not enough valid data.
-        print(f"\n{np.round((chi2_red < 1e-12).sum() / chi2_red.size, 1)}% of "
-              f"the objects have chi^2_red~0 and "
+        print(f"{np.round((chi2_red < 1e-12).sum() / chi2_red.size, 1)}% of "
+              f"the objects have χ²_red~0 and "
               f"{np.round((chi2_red < 0.5).sum() / chi2_red.size, 1)}% "
-              f"chi^2_red<0.5")
+              f"χ²_red<0.5")
 
 
 class ResultsManager(object):