From 67dd8aa0737bd362dab382500149ee9c3e8637ab Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 16 Nov 2020 08:29:32 -0600 Subject: use nlopt to find xopt This commit updates both the chi2 and dm-search scripts to run nlopt at the end of the MCMC to find the best fit point. The reason for this is that for the dm-search script, we rely on the best fit point being correct when doing the discovery threshold analysis. To really get the best fit point with the MCMC we would need to run a ridiculous number of steps, so it's better to run a fewer amount of steps to get close and then run the minimization from there. --- utils/chi2 | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) (limited to 'utils/chi2') diff --git a/utils/chi2 b/utils/chi2 index 12b602a..0326e1a 100755 --- a/utils/chi2 +++ b/utils/chi2 @@ -35,6 +35,7 @@ from sddm.dc import estimate_errors, EPSILON, truncnorm_scaled import emcee from sddm import printoptions from sddm.utils import fast_cdf, correct_energy_bias +import nlopt # Likelihood Fit Parameters # 0 - Atmospheric Neutrino Flux Scale @@ -390,7 +391,20 @@ def do_fit(data,muon,data_mc,atmo_scale_factor,muon_scale_factor,bins,steps,prin samples = sampler.get_chain(flat=True,thin=thin) - return sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()], samples + # Now, we use nlopt to find the best set of parameters. We start at the + # best starting point from the MCMC and then run the SBPLX routine. + x0 = sampler.get_chain(flat=True)[sampler.get_log_prob(flat=True).argmax()] + opt = nlopt.opt(nlopt.LN_SBPLX, len(x0)) + opt.set_min_objective(nll) + low = np.array(PRIORS_LOW) + high = np.array(PRIORS_HIGH) + opt.set_lower_bounds(low) + opt.set_upper_bounds(high) + opt.set_ftol_abs(1e-10) + opt.set_initial_step([0.01]*len(x0)) + xopt = opt.optimize(x0) + + return xopt, samples if __name__ == '__main__': import argparse -- cgit