41 views

Skip to first unread message

Apr 26, 2024, 6:39:20â€¯AMApr 26

to lmfit-py

Hi @mailinglist,

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.

mod = Model(invrecov)

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]

Â Â )

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

Apr 26, 2024, 8:53:50â€¯AMApr 26

to lmfit-py

Apologies for the lack of full traceback and fit report.

---------------------------------------------------------------------------
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

[[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

Apr 26, 2024, 12:41:47â€¯PMApr 26

to lmfi...@googlegroups.com

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

Â Â Â Â 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.

Reply all

Reply to author

Forward

0 new messages