diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py index f5479d00b6..fc9cf0822c 100644 --- a/validphys2/src/validphys/dataplots.py +++ b/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): diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index 3a9c3b3fb0..a49dc1a6fa 100644 --- a/validphys2/src/validphys/results.py +++ b/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)