# f(a) and f(b) must have different signs

41 views

### Timothy Anderson

Apr 26, 2024, 6:39:20â€¯AMApr 26
to lmfit-py
Hi @mailinglist,

I have started using lmfit in order to get some condensed information along with my fits (especially confidence intervals) but I am running into an error on some of my fits.

To make it short, I have 10 data points (with y_err) per fit, and 12 independent experiments. Each experiment are described by the same mathematical equation, so I loop over the 12 experiments and create a new lmfit model instance each time then fit the model instance with the 10 experimental data points associated to the experiment.

In 9 out of 12 cases, I can useÂ ci = result.conf_interval()Â  to get my confidence intervals. However, and for no obvious reason (as they all look kind of similar), 3 out of 12 give me aÂ ValueError: f(a) and f(b) must have different signs, but f(4.09...e-00) = 3.151...e-01, f(4.12...e-00) = 3.136...e-01 Â

associated toÂ File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/confidence.py:218, in ConfidenceInterval.calc_all_ci(self)

I was wondering if you had an idea as to how I could get this confidence interval computation to work, as for now I just have it in a try... except block so that the error is not fatal.

Note that depending on how the input data was preprocessed, some of the experiments used to not generate that error...

I read from the documentation and from online question boards that this may happen if the parameter being evaluated is close to its boundary but none of my 3 parameters have boundaries and their initial guess is somewhat close to the expected fit value.

Below a reduced sample of the code I'm using, with the data from a pathological case:

from lmfit import Model, printfuncs
from numpy import array

def invrecov(Ti, a, b, T1):
Â  Â  if T1 * 1e3 < 50: # Force T1 in valid range
Â  Â  Â  Â  return 1e7
Â  Â  return np.abs(a + b * np.exp(-Ti/T1))

fullMeans = array(
Â  Â  [2965.61052096, 1743.81963565, 1067.84200776,
Â  Â  Â 482.80343508, 59.56659829, 1540.88095421,
Â  Â  Â 4019.87643396, 4134.87512142, 4137.39869382,
Â  Â  Â 4145.13411215]
Â  Â  )
fullSEMs = array(
Â  Â  [0.13002311890348395, 0.09963669736078981, 0.07788739802154705,
Â  Â  Â 0.07584347244366009, 0.031136842850181763, 0.09370153130305435,
Â  Â  Â 0.15141882065887413, 0.15357232888620656, 0.1536191959417515,
Â  Â  Â 0.1537628769325479]
Â  Â  )

mod = Model(invrecov)
params = mod.make_params(a=4, b=-8, T1=0.1)
result = mod.fit(
np.array(fullMeans),
params,
Ti=np.array(time)/1e3,
weights=1.0/np.array(fullSEMs)
)
print(result.fit_report())
try:
Â  Â  ci = result.conf_interval()
Â  Â  print("[[Confidence Intervals]]")
Â  Â  printfuncs.report_ci(ci)
except ValueError:
pass

### Timothy Anderson

Apr 26, 2024, 8:53:50â€¯AMApr 26
to lmfit-py

Apologies for the lack of full traceback and fit report.
Here the full traceback:
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[133], line 21 18 result = mod.fit(np.array(fullMeans), params, Ti=np.array(time)/1e3, weights=1.0/np.array(fullSEMs)) print(result.fit_report()) ---> 21 ci = result.conf_interval() print("[[Confidence Intervals]]") 23 printfuncs.report_ci(ci) File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/model.py:1656, in ModelResult.conf_interval(self, **kwargs) def conf_interval(self, **kwargs): """Calculate the confidence intervals for the variable parameters. Confidence intervals are calculated using the (...) """ -> 1656 self.ci_out = conf_interval(self, self, **kwargs) return self.ci_out File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/confidence.py:143, in conf_interval(minimizer, result, p_names, sigmas, trace, maxiter, verbose, prob_func) 139 sigmas = [1, 2, 3] 141 ci = ConfidenceInterval(minimizer, result, p_names, prob_func, sigmas, 142 trace, verbose, maxiter) --> 143 output = ci.calc_all_ci() if trace: return output, ci.trace_dict File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/confidence.py:218, in ConfidenceInterval.calc_all_ci(self) 215 out = {} for p in self.p_names: --> 218 out[p] = (self.calc_ci(p, -1)[::-1] + 219 [(0., self.params[p].value)] + self.calc_ci(p, 1)) if self.trace: self.trace_dict = map_trace_to_names(self.trace_dict, self.params) File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/confidence.py:258, in ConfidenceInterval.calc_ci(self, para, direction) 255 ret.append((prob, direction*np.inf)) continue --> 258 sol = root_scalar(calc_prob, method='toms748', bracket=sorted([limit, a_limit]), rtol=.5e-4, args=(prob,)) if sol.converged: 260 val = sol.root File ~/anaconda3/envs/math/lib/python3.12/site-packages/scipy/optimize/_root_scalar.py:279, in root_scalar(f, args, method, bracket, fprime, fprime2, x0, x1, xtol, rtol, maxiter, options) 277 a, b = bracket[:2] try: --> 279 r, sol = methodc(f, a, b, args=args, **kwargs) except ValueError as e: # gh-17622 fixed some bugs in low-level solvers by raising an error # (rather than returning incorrect results) when the callable # returns a NaN. It did so by wrapping the callable rather than # modifying compiled code, so the iteration count is not available. if hasattr(e, "_x"): File ~/anaconda3/envs/math/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1394, in toms748(f, a, b, args, k, xtol, rtol, maxiter, full_output, disp) 1392 f = _wrap_nan_raise(f) 1393 solver = TOMS748Solver() -> 1394 result = solver.solve(f, a, b, args=args, k=k, xtol=xtol, rtol=rtol, 1395 maxiter=maxiter, disp=disp) 1396 x, function_calls, iterations, flag = result return _results_select(full_output, (x, function_calls, iterations, flag)) File ~/anaconda3/envs/math/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1240, in TOMS748Solver.solve(self, f, a, b, args, xtol, rtol, k, maxiter, disp) r"""Solve f(x) = 0 given an interval containing a root.""" self.configure(xtol=xtol, rtol=rtol, maxiter=maxiter, disp=disp, k=k) -> 1240 status, xn = self.start(f, a, b, args) if status == _ECONVERGED: return self.get_result(xn) File ~/anaconda3/envs/math/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1140, in TOMS748Solver.start(self, f, a, b, args) return _ECONVERGED, b if np.sign(fb) * np.sign(fa) > 0: -> 1140 raise ValueError("f(a) and f(b) must have different signs, but " "f({:e})={:e}, f({:e})={:e} ".format(a, fa, b, fb)) self.fab[:] = [fa, fb] return _EINPROGRESS, sum(self.ab) / 2.0 ValueError: f(a) and f(b) must have different signs, but f(-8.483104e+00)=1.044658e-03, f(-8.478857e+00)=9.178023e-04

Here the fit report:
[[Model]] Model(invrecov) [[Fit Statistics]] # fitting method = leastsq # function evals = 38 # data points = 10 # variables = 3 chi-square = 0.76333854 reduced chi-square = 0.10904836 Akaike info crit = -19.7263875 Bayesian info crit = -18.8186322 R-squared = 0.96833338 [[Variables]] a: 4.12311957 +/- 0.02568629 (0.62%) (init = 2) b: -8.38047280 +/- 0.05131575 (0.61%) (init = -4) T1: 0.20805633 +/- 0.00223634 (1.07%) (init = 1) [[Correlations]] (unreported correlations are < 0.100) C(b, T1) = +0.5472 C(a, T1) = +0.5385 C(a, b) = -0.3639 Calculating CI for a P(a=4.09743327967623) = 0.6489231252421723, max. prob=0.9973002039367398 P(a=4.071746985836196) = 0.994494268195663, max. prob=0.9973002039367398 P(a=4.0460606919961615) = 0.9976314999048017, max. prob=0.9973002039367398 Calculating CI for a P(a=4.148805867356299) = 0.6486322505079891, max. prob=0.9973002039367398 P(a=4.1744921611963335) = 0.9138781895005766, max. prob=0.9973002039367398 P(a=4.200178455036368) = 0.9798615529287257, max. prob=0.9973002039367398 P(a=4.225864748876402) = 0.9947399783060421, max. prob=0.9973002039367398 P(a=4.251551042716437) = 0.9984083715135426, max. prob=0.9973002039367398 Calculating CI for b P(b=-8.431788545236826) = 0.6483453541220551, max. prob=0.9973002039367398 P(b=-8.48310429428959) = 0.998344862193551, max. prob=0.9973002039367398

Moreover, the sample code above has a *1e-3 factor missing during the definition of fullMeans. It should read:
fullMeans = np.array([2965.61052096, 1743.81963565, 1067.84200776, Â 482.80343508,
Â  Â  Â  Â  Â 59.56659829, 1540.88095421, 4019.87643396, 4134.87512142,
Â  Â  Â  Â 4137.39869382, 4145.13411215])*1e-3

### Matt Newville

Apr 26, 2024, 12:41:47â€¯PMApr 26
I think that is signaling that `conf_interval()` call to `scipy.optimize.root_scalar` cannot find theÂ appropriate confidence level.Â

I was not able to run your example without guessing what the "times" value should be.Â  But.Â I'll guess that (as is often the case with fitting exponentials) that the uncertainties can be very large, and parameter values highly correlated.Â  In those cases, it can be hard to find confidence_levels by hand.

IÂ  am not sure we can always fix, that.Â  But, we should probablyÂ fail more gracefully here, and give a better message.

Separately, I would suggest removing your:
```
Â  Â  if T1 * 1e3 < 50: # Force T1 in valid range
Â  Â  Â  Â  return 1e7
```

and replacing it with a minimumÂ value for `T1`, perhaps something like:

```Â
Â  Â  params = mod.make_params(a=4, b=-8, T1=0.1)
Â  Â  params['T1'].min = 5.e-8 # note: AllowÂ the fit to explore smaller values.
```

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.