diff options
Diffstat (limited to 'utils')
| -rwxr-xr-x | utils/chi2 | 16 | ||||
| -rwxr-xr-x | utils/dm-search | 16 |
2 files changed, 30 insertions, 2 deletions
@@ -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 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(): """ |