added plotting function of prob inconsist consist

validphys2/src/validphys/closuretest/multiclosure_nsigma.py

@@ -183,7 +183,8 @@ def def_of_nsigma_alpha(
             fit_idxs = np.where(nsigma_values < z_alpha)[0]
         else:
             fit_idxs = np.where(nsigma_values > z_alpha)[0]
-        set1_alpha[alpha] = df.columns[fit_idxs].tolist()
+        # save it as set to allow for easy intersection with other sets
+        set1_alpha[alpha] = set(df.columns[fit_idxs])
 
     return NsigmaAlpha(alpha_dict=set1_alpha, is_weighted=multiclosurefits_nsigma.is_weighted)
 
@@ -326,9 +327,7 @@ def def_set_2(
         else:
             columns_bools = np.any((df_weight - df_ref).values > z_alpha, axis=0)
 
-        columns = df_weight.columns[columns_bools].to_list()
-
-        set2_alpha[alpha] = columns
+        set2_alpha[alpha] = set(df_weight.columns[columns_bools])
 
     return set2_alpha
 

validphys2/src/validphys/closuretest/multiclosure_nsigma_output.py

@@ -5,6 +5,62 @@
 import matplotlib.pyplot as plt
 
 
+@figure
+def plot_probability_inconsistent(set_2_alpha, set_1_alpha, set_3_alpha, comp_set_1_alpha, n_fits):
+    """
+    P(inconsistent) = P(set 2 | set 1) + P(set 3) + P(set 2 | ~ set 1)
+    """
+
+    fig, ax = plt.subplots()
+
+    tagged_rates = []
+    for alpha in set_2_alpha.keys():
+        set_2_inters_1 = set_2_alpha[alpha].intersection(set_1_alpha[alpha])
+        set_3 = set_3_alpha[alpha]
+        set_2_inters_comp_1 = set_2_alpha[alpha].intersection(comp_set_1_alpha[alpha])
+
+        set_tagged_fits = set_2_inters_1.union(set_3).union(set_2_inters_comp_1)
+        tagged_rates.append(len(set_tagged_fits) / n_fits)
+
+    ax.plot(set_2_alpha.keys(), tagged_rates, label="P(inconsistent)")
+    ax.legend()
+    return fig
+
+@figure
+def plot_probability_consistent(set_2_alpha, set_1_alpha, set_3_alpha, comp_set_1_alpha, comp_set_2_alpha, comp_set_3_alpha, n_fits):
+    """
+    P(consistent) = 1 - P(inconsistent) ?= P(~set 2 | ~set 1) + P(~set 2 and ~set 3| set 1)
+    """
+    fig, ax = plt.subplots()
+
+    untagged_rates = []
+    untagged_rates_method2 = []
+    for alpha in set_2_alpha.keys():
+        set_2_inters_1 = set_2_alpha[alpha].intersection(set_1_alpha[alpha])
+        set_3 = set_3_alpha[alpha]
+        set_2_inters_comp_1 = set_2_alpha[alpha].intersection(comp_set_1_alpha[alpha])
+
+        comp_set_2_inters_comp_set1 = comp_set_2_alpha[alpha].intersection(comp_set_1_alpha[alpha])
+        comp_set_2_inters_comp_set3 = comp_set_2_alpha[alpha].intersection(comp_set_3_alpha[alpha])
+        comp_set_2_inters_comp_set3_inter_set_1 = comp_set_2_inters_comp_set3.intersection(set_1_alpha[alpha])
+
+        set_untagged_fits = comp_set_2_inters_comp_set1.union(comp_set_2_inters_comp_set3_inter_set_1)
+
+
+        set_tagged_fits = set_2_inters_1.union(set_3).union(set_2_inters_comp_1)
+        untagged_rates.append(1 - len(set_tagged_fits) / n_fits)
+
+        untagged_rates_method2.append(len(set_untagged_fits) / n_fits)
+
+
+    ax.plot(set_2_alpha.keys(), untagged_rates, label="1 - P(inconsistent)")
+    ax.plot(set_2_alpha.keys(), untagged_rates_method2, label="P(~set 2 | ~set 1) + P(~set 2 and ~set 3| set 1)")
+    ax.legend()
+
+    return fig
+
+
+
 @figure
 def plot_set1_vs_set3_alpha(set_1_alpha, set_3_alpha, weighted_dataset, n_fits):
     """