diff options
Diffstat (limited to 'lmfit/confidence.py')
| -rw-r--r-- | lmfit/confidence.py | 75 |
1 files changed, 41 insertions, 34 deletions
diff --git a/lmfit/confidence.py b/lmfit/confidence.py index 61ddc78..14d26d5 100644 --- a/lmfit/confidence.py +++ b/lmfit/confidence.py @@ -14,24 +14,18 @@ CONF_ERR_STDERR = '%s without sensible uncertainty estimates' % CONF_ERR_GEN CONF_ERR_NVARS = '%s with < 2 variables' % CONF_ERR_GEN -def f_compare(ndata, nparas, new_chi, best_chi, nfix=1): - """Return the probability calculated using the F-test. +def f_compare(best_fit, new_fit): + """Return the probability calculated using the F-test. The null model (i.e., best-fit solution) is compared to an alternate model where one or more parameters are fixed. Parameters ---------- - ndata : int - Number of data points: :math:`N`. - nparas : int - Number of variables in the alternate model. - new_chi : float - Chi-square of the alternate model. - best_chi : float - Chi-square of the null model. - nfix : int - Number of fixed parameters (default is 1). + best_fit: MinimizerResult + The result from the best-fit. + new_fit: MinimizerResult + The result from fit with the fixed parameter(s). Returns ------- @@ -39,10 +33,9 @@ def f_compare(ndata, nparas, new_chi, best_chi, nfix=1): Value of the calculated probality. """ - nparas = nparas + nfix - nfree = ndata - nparas - nfix = 1.0*nfix - dchi = new_chi / best_chi - 1.0 + nfree = best_fit.nfree + nfix = best_fit.nvarys - new_fit.nvarys + dchi = new_fit.chisqr / best_fit.chisqr - 1.0 return f.cdf(dchi * nfree / nfix, nfix, nfree) @@ -161,7 +154,7 @@ def map_trace_to_names(trace, params): return out -class ConfidenceInterval(object): +class ConfidenceInterval: """Class used to calculate the confidence interval.""" def __init__(self, minimizer, result, p_names=None, prob_func=None, @@ -170,7 +163,7 @@ class ConfidenceInterval(object): self.verbose = verbose self.minimizer = minimizer self.result = result - self.params = result.params + self.params = result.params.copy() self.org = copy_vals(self.params) self.best_chi = result.chisqr @@ -182,20 +175,18 @@ class ConfidenceInterval(object): # check that there are at least 2 true variables! # check that all stderrs are sensible (including not None or NaN) - nvars = 0 + for par in self.fit_params: - if par.vary: - nvars += 1 - try: - if not (par.vary and par.stderr > 0.0): - raise MinimizerException(CONF_ERR_STDERR) - except TypeError: + if par.vary and (par.stderr is None or par.stderr is np.nan): raise MinimizerException(CONF_ERR_STDERR) + nvars = len([p for p in self.params.values() if p.vary]) if nvars < 2: raise MinimizerException(CONF_ERR_NVARS) - if prob_func is None or not hasattr(prob_func, '__call__'): + if prob_func is None: self.prob_func = f_compare + else: + self.prob_func = prob_func if trace: self.trace_dict = {i: [] for i in self.p_names} @@ -255,7 +246,6 @@ class ConfidenceInterval(object): try: val = brentq(calc_prob, a_limit, limit, rtol=.5e-4, args=prob) - except ValueError: self.reset_vals() try: @@ -293,16 +283,35 @@ class ConfidenceInterval(object): old_prob = 0 limit = start_val i = 0 + bound_reached = False + max_prob = max(self.probs) - while old_prob < max(self.probs): + while old_prob < max_prob: i = i + 1 limit += step * direction + if limit > para.max: + limit = para.max + bound_reached = True + elif limit < para.min: + limit = para.min + bound_reached = True new_prob = self.calc_prob(para, limit) - rel_change = (new_prob - old_prob) / max(new_prob, old_prob, 1.e-12) + rel_change = (new_prob - old_prob) / max(new_prob, old_prob, 1e-12) old_prob = new_prob + if self.verbose: + msg = "P({}={}) = {}, max. prob={}" + print(msg.format(para.name, limit, new_prob, max_prob)) # check for convergence + if bound_reached: + if new_prob < max(self.probs): + errmsg = ("Bound reached with " + "prob({}={}) = {} < max(sigmas)" + ).format(para.name, limit, new_prob) + warn(errmsg) + break + if i > self.maxiter: errmsg = "maxiter={} reached ".format(self.maxiter) errmsg += ("and prob({}={}) = {} < " @@ -331,8 +340,7 @@ class ConfidenceInterval(object): self.params[para.name] = para self.minimizer.prepare_fit(self.params) out = self.minimizer.leastsq() - prob = self.prob_func(out.ndata, out.ndata - out.nfree, - out.chisqr, self.best_chi) + prob = self.prob_func(self.result, out) if self.trace: x = [i.value for i in out.params.values()] @@ -391,7 +399,7 @@ def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10, best_chi = result.chisqr org = copy_vals(result.params) - if prob_func is None or not hasattr(prob_func, '__call__'): + if prob_func is None: prob_func = f_compare x = params[x_name] @@ -423,8 +431,7 @@ def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10, result.params[y.name] = y minimizer.prepare_fit(params=result.params) out = minimizer.leastsq() - prob = prob_func(out.ndata, out.ndata - out.nfree, out.chisqr, - best_chi, nfix=2.) + prob = prob_func(result, out) result.params[x.name] = save_x result.params[y.name] = save_y return prob |