Adding action to plot chi2 dist for aggregate of datasets

validphys2/src/validphys/dataplots.py

@@ -29,24 +29,55 @@
 log = logging.getLogger(__name__)
 
 @figure
-def plot_chi2dist(results, dataset, abs_chi2_data, chi2_stats, pdf):
+def plot_chi2dist_experiments(total_experiments_chi2data, experiments_chi2_stats, pdf):
+    """Plot the distribution of chi²s of the members of the pdfset."""
+    fig, ax = _chi2_distribution_plots(total_experiments_chi2data, experiments_chi2_stats, pdf, "hist")
+    ax.set_title(r"Experiments $\chi^2$ distribution")
+    return fig
+
+
+@figure
+def kde_chi2dist_experiments(total_experiments_chi2data, experiments_chi2_stats, pdf):
+    """KDE plot for experiments chi2."""
+    fig, ax = _chi2_distribution_plots(total_experiments_chi2data, experiments_chi2_stats, pdf, "kde")
+    ax.set_ylabel(r"Density")
+    ax.set_title(r"Experiments $\chi^2 KDE plot$")
+    return fig
+
+
+@figure
+def plot_chi2dist(dataset, abs_chi2_data, chi2_stats, pdf):
     """Plot the distribution of chi²s of the members of the pdfset."""
     setlabel = dataset.name
+    fig, ax = _chi2_distribution_plots(abs_chi2_data, chi2_stats, pdf, "hist")
+    ax.set_title(r"$\chi^2$ distribution for %s" % setlabel)
+    return fig
+
+
+def _chi2_distribution_plots(chi2_data, stats, pdf, plot_type):
     fig, ax = plt.subplots()
     label = pdf.name
-    alldata, central, npoints = abs_chi2_data
+    alldata, central, npoints = chi2_data
     if not isinstance(alldata, MCStats):
         ax.set_facecolor("#ffcccc")
         log.warning("Chi² distribution plots have a "
                 "different meaning for non MC sets.")
         label += " (%s!)" % pdf.ErrorType
-    label += '\n'+ '\n'.join(str(chi2_stat_labels[k])+(' %.2f' % v) for (k,v) in chi2_stats.items())
-    ax.set_title(r"$\chi^2$ distribution for %s" % setlabel)
+    label += '\n'+ '\n'.join(str(chi2_stat_labels[k])+(' %.2f' % v) for (k,v) in stats.items())
+    ax.set_xlabel(r"Replica $\chi^2$")
+
+    if plot_type == "hist":
+        ax.hist(alldata.data, label=label, zorder=100)
+    elif plot_type == "kde":
+        # We need the squeeze here to change shape from (x, 1) to (x,)
+        ax = plotutils.kde_plot(alldata.data.squeeze(), label=label)
+    else:
+        raise ValueError(f"plot_type must either be hist or kde, not {plot_type}")
 
-    ax.hist(alldata.data, label=label, zorder=100)
     l = ax.legend()
     l.set_zorder(1000)
-    return fig
+    return fig, ax
+
 
 @figure
 def plot_phi(experiments, experiments_phi):

validphys2/src/validphys/results.py

@@ -819,8 +819,33 @@ def count_negative_points(possets_predictions):
     'chi2_per_data': r'$\frac{\chi^2}{N_{data}}$'
 }
 
+def experiments_chi2_stats(total_experiments_chi2data):
+    """Compute several estimators from the chi² for an
+    aggregate of experiments:
+
+     - central_mean
+
+     - npoints
+
+     - perreplica_mean
+
+     - perreplica_std
+
+     - chi2_per_data
+    """
+    rep_data, central_result, npoints = total_experiments_chi2data
+    m = central_result.mean()
+    rep_mean = rep_data.central_value().mean()
+    return OrderedDict([
+            ('central_mean'        ,  m),
+            ('npoints'             , npoints),
+            ('chi2_per_data', m/npoints),
+            ('perreplica_mean', rep_mean),
+            ('perreplica_std',  rep_data.std_error().mean()),
+           ])
+
 def chi2_stats(abs_chi2_data):
-    """Compute severa estimators from the chi²:
+    """Compute several estimators from the chi²:
 
      - central_mean
 
@@ -1071,11 +1096,10 @@ def total_experiments_chi2data(pdf: PDF, experiments_chi2):
     ndata = 0
     central_chi2 = 0
     nmembers = len(experiments_chi2[0].replica_result.error_members())  # ugh
-    member_chi2  = np.zeros(nmembers)
+    member_chi2  = np.zeros((nmembers, 1))
 
     for cd in experiments_chi2:
-        # not sure why the transpose or [0] are needed here
-        member_chi2  += cd.replica_result.error_members().T[0]
+        member_chi2  += cd.replica_result.error_members()
         central_chi2 += cd.central_result
         ndata += cd.ndata
     return Chi2Data(pdf.stats_class(member_chi2), central_chi2, ndata)