conf_interval NoneType error

[Ben Kupernik]

Hi,

I have a bunch of ModleResult objects that I saved to files using the save_modleresult function a while back. I need to get error bars on the amplitudes of certain peaks in these models. My thought is to load the models using the load_modleresults function and then use the conf_interval function to get the error at one standard deviation. The issue I am running into is when I call the conf_interval function it returns TypeError NoneType object is not iterable. I tried result.ci_report() and it resulted in the same error.

The file I am loading is attached, any help would be appreciated.

Here is the code

import lmfit as lm

result = lm.model.load_modelresult('model_result')

# manually set the error per the Calculation of confidence intervals documentation. I don't think this is relevant but i tried it just in case, it didn't work.

# for p in result.params:

#     result.params[p].stderr = abs(result.params[p].value * 0.1)

lm.conf_interval(result, result) # per the discussion here https://groups.google.com/g/lmfit-py/c/WZmntp8eaqI

Here is the fit report

[[Model]]
    (Model(linear, prefix='b_') + Model(pvoigt, prefix='v_bkg_'))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 71
    # data points      = 38
    # variables        = 6
    chi-square         = 3.52416173
    reduced chi-square = 0.11013005
    Akaike info crit   = -78.3618552
    Bayesian info crit = -68.5363382
    R-squared          = 0.97230752
[[Variables]]
    b_slope:          32.7574758 +/- 2.67122678 (8.15%) (init = 75)
    b_intercept:      12.9899327 +/- 6.43488238 (49.54%) (init = 50)
    v_bkg_amplitude:  0.17590785 +/- 0.07327093 (41.65%) (init = 0.2085684)
    v_bkg_center:     2.53268932 +/- 0.00137492 (0.05%) (init = 2.53288)
    v_bkg_sigma:      0.02067294 +/- 0.00327488 (15.84%) (init = 0.02055412)
    v_bkg_fraction:   0.87772000 +/- 0.52212039 (59.49%) (init = 0.5)
    v_bkg_fwhm:       0.04134588 +/- 0.00654977 (15.84%) == '2.0000000*v_bkg_sigma'
    v_bkg_height:     2.86606585 +/- 0.27935328 (9.75%) == '(((1-v_bkg_fraction)*v_bkg_amplitude)/max(1e-15, (v_bkg_sigma*sqrt(pi/log(2))))+(v_bkg_fraction*v_bkg_amplitude)/max(1e-15, (pi*v_bkg_sigma)))'
[[Correlations]] (unreported correlations are < 0.100)
    C(b_slope, b_intercept)            = -0.9992
    C(v_bkg_amplitude, v_bkg_fraction) = +0.9359
    C(v_bkg_amplitude, v_bkg_sigma)    = +0.7385
    C(b_slope, v_bkg_amplitude)        = -0.7339
    C(b_intercept, v_bkg_amplitude)    = +0.7068
    C(b_slope, v_bkg_sigma)            = -0.6705
    C(b_intercept, v_bkg_sigma)        = +0.6543
    C(b_intercept, v_bkg_center)       = +0.6158
    C(b_slope, v_bkg_center)           = -0.6150
    C(b_slope, v_bkg_fraction)         = -0.5827
    C(b_intercept, v_bkg_fraction)     = +0.5543
    C(v_bkg_sigma, v_bkg_fraction)     = +0.4878
    C(v_bkg_amplitude, v_bkg_center)   = +0.4348
    C(v_bkg_center, v_bkg_sigma)       = +0.3934

    C(v_bkg_center, v_bkg_fraction)    = +0.3484

Here is the error report

{
"name": "TypeError",
"message": "'NoneType' object is not iterable",
"stack": "---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[32], line 1
----> 1 lm.conf_interval(result, result)

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/confidence.py:146, in conf_interval(minimizer, result, p_names, sigmas, trace, maxiter, verbose, prob_func, min_rel_change)
    142     sigmas = [1, 2, 3]
    144 ci = ConfidenceInterval(minimizer, result, p_names, prob_func, sigmas,
    145                         trace, verbose, maxiter, min_rel_change)
--> 146 output = ci.calc_all_ci()
    147 if trace:
    148     return output, ci.trace_dict

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/confidence.py:257, in ConfidenceInterval.calc_all_ci(self)
    254 out = {}
    256 for p in self.p_names:
--> 257     out[p] = (self.calc_ci(p, -1)[::-1] +
    258               [(0., self.params[p].value)] +
    259               self.calc_ci(p, 1))
    260 if self.trace:
    261     self.trace_dict = map_trace_to_names(self.trace_dict, self.params)

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/confidence.py:287, in ConfidenceInterval.calc_ci(self, para, direction)
    284     self.trace_dict[para.name].append(x + [0])
    286 para.vary = False
--> 287 limit, max_prob = self.find_limit(para, direction)
    288 a_limit = float(para.value)
    289 ret = []

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/confidence.py:345, in ConfidenceInterval.find_limit(self, para, direction)
    342     limit = para.min
    343     bound_reached = True
--> 345 new_prob = self.calc_prob(para, limit)
    346 rel_change = (new_prob - old_prob) / max(new_prob, old_prob, 1e-12)
    347 old_prob = new_prob

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/confidence.py:384, in ConfidenceInterval.calc_prob(self, para, val, offset, restore)
    382 self.params[para.name] = para
    383 self.minimizer.prepare_fit(self.params)
--> 384 out = self.minimizer.leastsq()
    385 prob = self.prob_func(self.result, out)
    387 if self.trace:

File /opt/anaconda3/envs/lmfit/lib/python3.10/site-packages/lmfit/minimizer.py:1635, in Minimizer.leastsq(self, params, max_nfev, **kws)
   1631     warnings.warn(maxeval_warning.format('maxfev', thisfuncname()),
   1632                   RuntimeWarning)
   1633     kws.pop('maxfev')
-> 1635 lskws.update(self.kws)
   1636 lskws.update(kws)
   1637 self.col_deriv = False

TypeError: 'NoneType' object is not iterable"
}

[Matt Newville]

Hi Ben,

When I load that model result and print the report, I do not get what you reported, I get

[[Model]]
    (Model(linear, prefix='b_') + Model(pvoigt, prefix='v0_'))

[[Fit Statistics]]
    # fitting method   = leastsq

    # function evals   = 14000
    # data points      = 44
    # variables        = 6
    chi-square         = 44.1551231
    reduced chi-square = 1.16197692
    Akaike info crit   = 12.1548503
    Bayesian info crit = 22.8599881
    R-squared          = 0.99473310
##  Warning: uncertainties could not be estimated:
    v0_sigma:      at boundary
    v0_fraction:   at boundary
[[Variables]]
    b_slope:      -119.497805 (init = 0)
    b_intercept:   349.578616 (init = 50)
    v0_amplitude:  1.79206112 (init = 1.042208)
    v0_center:     1.68336994 (init = 1.679341)
    v0_sigma:      0.02146012 (init = 0.0171681)
    v0_fraction:   5.5511e-17 (init = 0.5)
    v0_fwhm:       0.04292025 == '2.0000000*v0_sigma'
    v0_height:     39.2245872 == '(((1-v0_fraction)*v0_amplitude)/max(1e-15, (v0_sigma*sqrt(pi/log(2))))+(v0_fraction*v0_amplitude)/max(1e-15, (pi*v0_sigma)))'

If the uncertainties could not be estimated, that explains why `conf_intervals()` won't work.

And, indeed,  doing 

 

```

for par in result.params.values(): print(par)

```

gives

```

<Parameter 'b_slope', value=-119.49780522461892, bounds=[-inf:inf]>
<Parameter 'b_intercept', value=349.5786158310988, bounds=[-inf:inf]>
<Parameter 'v0_amplitude', value=1.7920611201690926, bounds=[0:2.084415661005683]>
<Parameter 'v0_center', value=1.6833699404604974, bounds=[1.6693413236136088:1.6893413236136088]>
<Parameter 'v0_sigma', value=0.02146012414271572, bounds=[0:0.02146012418638302]>
<Parameter 'v0_fraction', value=5.551115123125783e-17, bounds=[0.0:1.0]>
<Parameter 'v0_fwhm', value=0.04292024828543144, bounds=[-inf:inf], expr='2.0000000*v0_sigma'>
<Parameter 'v0_height', value=39.22458721112634, bounds=[-inf:inf], expr='(((1-v0_fraction)*v0_amplitude)/max(1e-15, (v0_sigma*sqrt(pi/log(2))))+(v0_fraction*v0_amplitude)/max(1e-15, (pi*v0_sigma)))'>

```

which strongly suggests that you programmatically set very tight bounds.   I cannot express strongly enough what a terrible idea this is.  Here, you are asking for code to find an optimal value for `v0_sigma` and also saying that you know with absolute certainty that it cannot be larger than 0.02146012418638302.   A value of 0.0214602 is simply out of the question, you said.  The response was, "Sorry, I can't help you".   

If you want to get reliable values and estimate uncertainties, do not do that.   Use bounds to prevent truly unfeasible values.    Let the fit do its job. 

FWIW,  if you do get a good fit for these data, I would be very surprised if the 1-sigma confidence intervals from `conf_interval()` were significantly different from the automatically estimated 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 visit https://groups.google.com/d/msgid/lmfit-py/79efbb2d-f2b7-4ef9-af09-33a55d26cb73n%40googlegroups.com.

--Matt