From d0261d21c76a7d2fe740840d188775504195813e Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 31 Aug 2020 09:58:48 -0500 Subject: 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. --- utils/chi2 | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) (limited to 'utils/chi2') 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 -- cgit