diff options
Diffstat (limited to 'examples/doc_fitting_emcee.py')
| -rw-r--r-- | examples/doc_fitting_emcee.py | 47 |
1 files changed, 38 insertions, 9 deletions
diff --git a/examples/doc_fitting_emcee.py b/examples/doc_fitting_emcee.py index 1d396d6..5862804 100644 --- a/examples/doc_fitting_emcee.py +++ b/examples/doc_fitting_emcee.py @@ -8,6 +8,7 @@ try: HASPYLAB = True except ImportError: HASPYLAB = False +HASPYLAB = False try: import corner @@ -46,27 +47,45 @@ if HASPYLAB: mi.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2)) res = lmfit.minimize(residual, method='emcee', nan_policy='omit', burn=300, - steps=1000, thin=20, params=mi.params, is_weighted=False) + steps=1000, thin=20, params=mi.params, is_weighted=False, + progress=False) if HASPYLAB and HASCORNER: emcee_corner = corner.corner(res.flatchain, labels=res.var_names, truths=list(res.params.valuesdict().values())) plt.show() +if HASPYLAB: + plt.plot(res.acceptance_fraction) + plt.xlabel('walker') + plt.ylabel('acceptance fraction') + plt.show() + +if hasattr(res, "acor"): + print("Autocorrelation time for the parameters:") + print("----------------------------------------") + for i, par in enumerate(p): + print(par, res.acor[i]) + print("\nmedian of posterior probability distribution") print('--------------------------------------------') lmfit.report_fit(res.params) + # find the maximum likelihood solution highest_prob = np.argmax(res.lnprob) hp_loc = np.unravel_index(highest_prob, res.lnprob.shape) mle_soln = res.chain[hp_loc] for i, par in enumerate(p): p[par].value = mle_soln[i] -print("\nMaximum likelihood Estimation") -print('-----------------------------') -for _, vals in p.items(): - print(vals) + +print('\nMaximum Likelihood Estimation from emcee ') +print('-------------------------------------------------') +print('Parameter MLE Value Median Value Uncertainty') +fmt = ' {:5s} {:11.5f} {:11.5f} {:11.5f}'.format +for name, param in p.items(): + print(fmt(name, param.value, res.params[name].value, + res.params[name].stderr)) if HASPYLAB: plt.figure() @@ -76,7 +95,17 @@ if HASPYLAB: plt.legend() plt.show() -quantiles = np.percentile(res.flatchain['t1'], [2.28, 15.9, 50, 84.2, 97.7]) -print("1 sigma spread", 0.5 * (quantiles[3] - quantiles[1])) -print("2 sigma spread", 0.5 * (quantiles[4] - quantiles[0])) -# <end of examples/doc_fitting_emcee.py> +print('\nError Estimates from emcee ') +print('------------------------------------------------------') +print('Parameter -2sigma -1sigma median +1sigma +2sigma ') + +for name in p.keys(): + quantiles = np.percentile(res.flatchain[name], + [2.275, 15.865, 50, 84.135, 97.275]) + median = quantiles[2] + err_m2 = quantiles[0] - median + err_m1 = quantiles[1] - median + err_p1 = quantiles[3] - median + err_p2 = quantiles[4] - median + fmt = ' {:5s} {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f}'.format + print(fmt(name, err_m2, err_m1, median, err_p1, err_p2)) |