add sddm/stats.py

This commit adds the new file sddm/stats.py to and adds a function to
correctly sample a Monte Carlo histogram when computing p-values. In
particular, I now take into account the uncertainty on the total number
of expected events by drawing from a gamma distribution, which is the
posterior of the Poisson likelihood function with a prior of 1/lambda.

1 files changed, 13 insertions, 9 deletions

diff --git a/utils/chi2 b/utils/chi2
index 0460993..72ad89b 100755
--- a/utils/chi2
+++ b/utils/chi2
@@ -32,6 +32,7 @@ from matplotlib.lines import Line2D
 from scipy.stats import iqr, norm, beta
 from scipy.special import spence
 from itertools import izip_longest
+from sddm.stats import *
 
 # Uncertainty on the energy scale
 # FIXME: Should get real number from stopping muons
@@ -124,16 +125,11 @@ def chi2(samples,expected):
     return np.sum((samples-expected)**2/expected,axis=-1)
 
 def get_mc_hist(data,x,bins):
-    if len(data):
-        return np.histogram(data,bins=bins)[0]*x[0]/100.0
-    else:
-        return np.zeros(len(bins)-1,dtype=np.int)
+    hist = np.histogram(data,bins=bins)[0]
+    return sample_mc_hist(hist,norm=x[0]/100.0)
 
 def get_data_hist(data,x,bins):
-    if len(data):
-        return np.histogram(data*x[1],bins=bins)[0]
-    else:
-        return np.zeros(len(bins)-1,dtype=np.int)
+    return np.histogram(data*x[1],bins=bins)[0]
 
 def get_multinomial_prob(data, data_mc, x_samples, bins, size=10000):
     """
@@ -190,7 +186,15 @@ if __name__ == '__main__':
 
     import matplotlib.pyplot as plt
 
-    ev = get_events(args.filenames,merge_fits=True,nhit_thresh=args.nhit_thresh)
+    # Loop over runs to prevent using too much memory
+    evs = []
+    rhdr = pd.concat([read_hdf(filename, "rhdr").assign(filename=filename) for filename in args.filenames],ignore_index=True)
+    for run, df in rhdr.groupby('run'):
+        evs.append(get_events(df.filename.values, merge_fits=True, nhit_thresh=args.nhit_thresh))
+    ev = pd.concat(evs)
+
+    ev_mc = get_events(args.mc, merge_fits=True)
+
     ev_mc = get_events(args.mc,merge_fits=True,nhit_thresh=args.nhit_thresh)
 
     ev = ev.reset_index()