add estimate_errors to chi2 analysis

This commit updates the estimate_errors() function so that it works
without a list of constraints and uses arrays of low and high bounds
passed in instead of hardcoded constraints.

I can now call this function from the chi2 analysis to get the stepsizes
before running the MCMC.

1 files changed, 18 insertions, 1 deletions

diff --git a/utils/chi2 b/utils/chi2
index 9970bfd..b1c3c17 100755
--- a/utils/chi2
+++ b/utils/chi2
@@ -31,6 +31,18 @@ from scipy.stats import iqr, norm, beta
 from scipy.special import spence
 from itertools import izip_longest
 from sddm.stats import *
+import contextlib
+from sddm.dc import estimate_errors
+
+# from https://stackoverflow.com/questions/2891790/how-to-pretty-print-a-numpy-array-without-scientific-notation-and-with-given-pre
+@contextlib.contextmanager
+def printoptions(*args, **kwargs):
+    original = np.get_printoptions()
+    np.set_printoptions(*args, **kwargs)
+    try:
+        yield
+    finally: 
+        np.set_printoptions(**original)
 
 # Uncertainty on the energy scale
 # FIXME: These are just placeholders! Should get real number from stopping
@@ -364,9 +376,14 @@ def do_fit(data,muon,data_mc,bins,steps):
     nll_xopt = nll(xopt)
     print("nll(xopt) = ", nll(xopt))
 
+    stepsizes = estimate_errors(nll,xopt,low,high,constraints)
+
+    with printoptions(precision=3, suppress=True):
+        print("Errors: ", stepsizes)
+
     pos = np.empty((20, len(x0)),dtype=np.double)
     for i in range(pos.shape[0]):
-        pos[i] = xopt + np.random.randn(len(x0))*xopt*0.1
+        pos[i] = xopt + np.random.randn(len(x0))*stepsizes
         pos[i,:] = np.clip(pos[i,:],low,high)
 
     nwalkers, ndim = pos.shape