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/dm-search | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) (limited to 'utils/dm-search') diff --git a/utils/dm-search b/utils/dm-search index 2a008d5..4e8772e 100755 --- a/utils/dm-search +++ b/utils/dm-search @@ -36,6 +36,7 @@ from sddm.utils import fast_cdf, correct_energy_bias from scipy.integrate import quad from sddm.dm import * from sddm import SNOMAN_MASS +import nlopt # Likelihood Fit Parameters # 0 - Atmospheric Neutrino Flux Scale @@ -312,7 +313,20 @@ def do_fit(dm_particle_id,dm_mass,dm_energy,data,muon,data_mc,atmo_scale_factor, 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 def sample_priors(): """ -- cgit