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

[Timothy Anderson]

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)

which is called by File ~/anaconda3/envs/math/lib/python3.12/site-packages/lmfit/confidence.py:143,

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]

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)) 19 print(result.fit_report()) ---> 21 ci = result.conf_interval() 22 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) 1646 def conf_interval(self, **kwargs): 1647 """Calculate the confidence intervals for the variable parameters. 1648 1649 Confidence intervals are calculated using the (...) 1654 1655 """ -> 1656 self.ci_out = conf_interval(self, self, **kwargs) 1657 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() 144 if trace: 145 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 = {} 217 for p in self.p_names: --> 218 out[p] = (self.calc_ci(p, -1)[::-1] + 219 [(0., self.params[p].value)] + 220 self.calc_ci(p, 1)) 221 if self.trace: 222 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)) 256 continue --> 258 sol = root_scalar(calc_prob, method='toms748', bracket=sorted([limit, a_limit]), rtol=.5e-4, args=(prob,)) 259 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] 278 try: --> 279 r, sol = methodc(f, a, b, args=args, **kwargs) 280 except ValueError as e: 281 # gh-17622 fixed some bugs in low-level solvers by raising an error 282 # (rather than returning incorrect results) when the callable 283 # returns a NaN. It did so by wrapping the callable rather than 284 # modifying compiled code, so the iteration count is not available. 285 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 1397 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) 1238 r"""Solve f(x) = 0 given an interval containing a root.""" 1239 self.configure(xtol=xtol, rtol=rtol, maxiter=maxiter, disp=disp, k=k) -> 1240 status, xn = self.start(f, a, b, args) 1241 if status == _ECONVERGED: 1242 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) 1137 return _ECONVERGED, b 1139 if np.sign(fb) * np.sign(fa) > 0: -> 1140 raise ValueError("f(a) and f(b) must have different signs, but " 1141 "f({:e})={:e}, f({:e})={:e} ".format(a, fa, b, fb)) 1142 self.fab[:] = [fa, fb] 1144 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]

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.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/39d885d0-3672-47fe-8306-3704f022625bn%40googlegroups.com.

--

--Matt Newville <newville at cars.uchicago.edu> 630-327-7411