Uncertainty on fitted parameters

[Daniel Marx]

I'm looking for the easiest way of outputting the uncertainty in the fitted parameters. With spo.curve_fit, we just get the covariance matrix when we fit and we can take the diagonal and square root to find the uncertainties. With lmfit it doesn't seem to be so simple.

My fitting looks like this:
**************************
a_lm2 = lmfit.Parameter('a', value=a_est)

b_lm2 = lmfit.Parameter('b', value=b_est)

x0_core_lm2 = lmfit.Parameter('x0_core', value=gaus1['x0_core'])

x0_1_lm2 = lmfit.Parameter('x0_1', value=gaus1['x0_1'])

x0_2_lm2 = lmfit.Parameter('x0_2', value=gaus1['x0_2'])

x0_3_lm2 = lmfit.Parameter('x0_3', value=gaus1['x0_3'])

x0_4_lm2 = lmfit.Parameter('x0_4', value=gaus1['x0_4'])

sig_core_lm2 = lmfit.Parameter('sig_core', value=gaus1['sig_core'])

sig_1_lm2 = lmfit.Parameter('sig_1', value=gaus1['sig_1'])

sig_2_lm2 = lmfit.Parameter('sig_2', value=gaus1['sig_2'])

sig_3_lm2 = lmfit.Parameter('sig_3', value=gaus1['sig_3'])

sig_4_lm2 = lmfit.Parameter('sig_4', value=gaus1['sig_4'])

m_lm2 = lmfit.Parameter('m', value=m, vary=False)

c_lm2 = lmfit.Parameter('c', value=c, vary=False)

gausfit2 = mod.fit(y, x=x, a=a_lm2, b=b_lm2, x0_core=x0_core_lm2, x0_1=x0_1_lm2, x0_2=x0_2_lm2,

x0_3=x0_3_lm2, x0_4=x0_4_lm2, sig_core=sig_core_lm2, sig_1=sig_1_lm2, sig_2=sig_2_lm2,

sig_3=sig_3_lm2, sig_4=sig_4_lm2, m=m_lm2, c=c_lm2,weights=None, scale_covar=False)

print 'a_lm2_unc =', a_lm2.stderr

*********************************

When I generate a fit report, I get uncertainty values so they are clearly being computed. My problem is calling them and using them. I tried to just print the uncertainty of the parameter using the stderr as in the last line of code above but this just returns 'None'. I can get a covariance matrix but I don't know what order this is being displayed in. My ultimate goal is simply to have the values and associated uncertainties that I can then put into an array and use further in my code.

[Matt Newville]

Hi Daniel,

On Thu, May 7, 2015 at 7:15 AM, Daniel Marx <natura...@gmail.com> wrote:

I'm looking for the easiest way of outputting the uncertainty in the fitted parameters. With spo.curve_fit, we just get the covariance matrix when we fit and we can take the diagonal and square root to find the uncertainties. With lmfit it doesn't seem to be so simple.

It should be simple.

My fitting looks like this:
**************************
a_lm2 = lmfit.Parameter('a', value=a_est)

b_lm2 = lmfit.Parameter('b', value=b_est)

x0_core_lm2 = lmfit.Parameter('x0_core', value=gaus1['x0_core'])

x0_1_lm2 = lmfit.Parameter('x0_1', value=gaus1['x0_1'])

x0_2_lm2 = lmfit.Parameter('x0_2', value=gaus1['x0_2'])

x0_3_lm2 = lmfit.Parameter('x0_3', value=gaus1['x0_3'])

x0_4_lm2 = lmfit.Parameter('x0_4', value=gaus1['x0_4'])

sig_core_lm2 = lmfit.Parameter('sig_core', value=gaus1['sig_core'])

sig_1_lm2 = lmfit.Parameter('sig_1', value=gaus1['sig_1'])

sig_2_lm2 = lmfit.Parameter('sig_2', value=gaus1['sig_2'])

sig_3_lm2 = lmfit.Parameter('sig_3', value=gaus1['sig_3'])

sig_4_lm2 = lmfit.Parameter('sig_4', value=gaus1['sig_4'])

m_lm2 = lmfit.Parameter('m', value=m, vary=False)

c_lm2 = lmfit.Parameter('c', value=c, vary=False)

gausfit2 = mod.fit(y, x=x, a=a_lm2, b=b_lm2, x0_core=x0_core_lm2, x0_1=x0_1_lm2, x0_2=x0_2_lm2,

x0_3=x0_3_lm2, x0_4=x0_4_lm2, sig_core=sig_core_lm2, sig_1=sig_1_lm2, sig_2=sig_2_lm2,

sig_3=sig_3_lm2, sig_4=sig_4_lm2, m=m_lm2, c=c_lm2,weights=None, scale_covar=False)

This is not at all the recommended approach.    You don't need to create a lmfit.Parameter  to pass as the argument to Model.fit -- just pass in the initial values.   Those Parameter objects won't get used.

To see the uncertainties and correlations with other parameters, use the Parameters ordered dict that is held in 'gausfit2.params', perhaps with:

for param in gausfit2.params.values():

    print param.name, param.value, param.stderr, param.correl

If you do want the covariance matrix, it's in 'gausfit2.covar'

Hope that helps,

--Matt

[Daniel Marx]

Thanks a lot Matt!

 

I see how to obtain the standard error values now. On your first point that my method is not at all the recommended approach, the reason I wasn’t just passing in the initial values is because I want to implement bounds on my parameters later. How would you recommend doing this?

 

Thanks again,

Daniel

[Matt Newville]

On Sat, May 9, 2015 at 1:45 PM, Daniel Marx <natura...@gmail.com> wrote:

Thanks a lot Matt!

 

I see how to obtain the standard error values now. On your first point that my method is not at all the recommended approach, the reason I wasn’t just passing in the initial values is because I want to implement bounds on my parameters later. How would you recommend doing this?

 

Well, you can either set the "hints" for the parameters for Model, as with

    mod = Model(MyFunction)
    mod.set_param_hint('alpha',  value=0.5, min=0)

    mod.set_param_hint('sigma',  value=0.1, min=0.0001)

    mod.set_param_hint('fwhm', expr='2*sigma')

    mod.fit(data)

In this approach, the bounds and constraint expressions become part of the model, and are applied every time you create (even implicitly) a set of Parameters to use in a fit.

An alternative is explicitly create parameters for a model and then add the bounds / constraints to these, as with:

    mod = Model(SomeOtherFunction)
    params = mod.make_params()
    params['width'].min = 0

    mod.fit(data, params)

Which of these two approaches is used is mostly a matter of taste.  The "hints" approach sort of applies the bounds/constraints are inherent in the model, as if it makes no physical sense for some value to be out-of-bounds.

See http://lmfit.github.io/lmfit-py/model.html#using-parameter-hints and http://lmfit.github.io/lmfit-py/model.html#model.Model.set_param_hint for more details.

Hope that helps,

--Matt

[Daniel Marx]

Thanks Matt. So I managed to get this working using the first method.

 

However, when I implement bounds I sometimes find that the stderr is given as 0. This seems to be the case when I have more restrictive bounds and the parameters end up equal to the bounded values. Do you know why this is?

 

From: matt.n...@gmail.com [mailto:matt.n...@gmail.com] On Behalf Of Matt Newville
Sent: samedi 9 mai 2015 22:54
To: Daniel Marx
Cc: lmfit-py
Subject: Re: Uncertainty on fitted parameters

 

 

 

On Sat, May 9, 2015 at 1:45 PM, Daniel Marx <natura...@gmail.com> wrote:

[Matt Newville]

Hi Daniel,

On Tue, May 12, 2015 at 9:44 AM, Daniel Marx <natura...@gmail.com> wrote:

Thanks Matt. So I managed to get this working using the first method.

 

However, when I implement bounds I sometimes find that the stderr is given as 0. This seems to be the case when I have more restrictive bounds and the parameters end up equal to the bounded values. Do you know why this is?

If the best fit value for a parameter is very close to the boundary, the effect on the fit quality of moving the value closer to that boundary becomes very hard to estimate.  The partial derivative with respect to that parameter will be unreliable and the covariance matrix will become singular.

If a best-fit value is at the boundary, you sort of have to ask yourself which you believe more:  the boundary itself, or what the data is trying to tell you.  That's not to say that boundaries aren't useful or often warranted.  But a boundary is (and, by design) preventing values outside the boundary from being explored.

--Matt

[Scott Norris]

I've frequently run into this lack of uncertainties when the fitted values are near a boundary.  In my case it often happens when 

(1) I've set any boundary at all

(2) my initial guess is not good enough

Under these circumstances it seems like the fitting method is wandering away from the true minimum, until it encounters a downward slope that passes through the boundary, and gets stuck.

Is there a robust way to determine when this situation has been encountered?  Often in these situations I can identify the true minimum with a better guess, but I'd like to know when that occurs.  At the moment I am simply checking to see if any of the .stderr values are zero, but this doesn't seem entirely reliable.  Is there (or could there be) some boolean flag indicating that uncertainties could not be calculated?

Scott

[Matt Newville]

Hi Scott,

On Wed, May 13, 2015 at 4:27 PM, Scott Norris <scott.no...@gmail.com> wrote:

I've frequently run into this lack of uncertainties when the fitted values are near a boundary.  In my case it often happens when 

(1) I've set any boundary at all

(2) my initial guess is not good enough

Under these circumstances it seems like the fitting method is wandering away from the true minimum, until it encounters a downward slope that passes through the boundary, and gets stuck.

I'm not entirely sure when that would happen.  It seems pretty likely to depend strongly on the problem.

 

Is there a robust way to determine when this situation has been encountered?  Often in these situations I can identify the true minimum with a better guess, but I'd like to know when that occurs.  At the moment I am simply checking to see if any of the .stderr values are zero, but this doesn't seem entirely reliable.  Is there (or could there be) some boolean flag indicating that uncertainties could not be calculated?

Minimizer does have an 'errorbars' attribute which contains whether error bars could be estimated.

--Matt

[Scott Norris]

Okay, great!  I'll look into that.  I'm sure this behavior is problem-dependent, but this will help a bunch.

Scott

[Daniel Marx]

Btw, when we provide weights in Model.fit (which I assume should be provided as 1/sigma^2), do the magnitude of these affect the outputted standard errors from the fit? In spo.curve_fit, we have the parameter absolute_sigma that means that sigma describes one standard deviation so errors are assumed to be absolute not relative. Can I get the same thing in lmfit?

[Matt Newville]

Hi Daniel,

On Thu, May 21, 2015 at 3:01 AM, Daniel Marx <natura...@gmail.com> wrote:

Btw, when we provide weights in Model.fit (which I assume should be provided as 1/sigma^2), do the magnitude of these affect the outputted standard errors from the fit? In spo.curve_fit, we have the parameter absolute_sigma that means that sigma describes one standard deviation so errors are assumed to be absolute not relative. Can I get the same thing in lmfit?

The weights should normally be 1/sigma, not 1/sigma^2.  I think the scipy.optimize.curve_fit documentation may be mistaken about that, but I haven't looked in detail at that.    In lmfit, the array "weight*(model-data)" is minimized in the least squares sense.

By default, lmfit does rescale the covariance matrix so that the weights are "relative", not "absolute".  In lmfit, this is optional but called "scale_covar" (defaulting to True).   You can set this to False (akin to absolute_sigma=True) in "minimize()" or "Model.fit()".

--Matt
 

 

--
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 post to this group, send email to lmfi...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/08beb955-0cde-4dbb-94cd-564f1ba62c43%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

--Matt Newville <newville at cars.uchicago.edu> 630-252-0431

[Daniel Marx]

Hi Matt,

I was just looking again at this. For some reason, the uncertainties I'm getting from fits with lmfit and spo.curve_fit are quite different. I've set them so both should be taking absolute errors as the input. Any ideas why this might be?

[Matt Newville]

Hi Daniel,

On Jun 15, 2015 8:33 AM, "Daniel Marx" <natura...@gmail.com> wrote:
>
> Hi Matt,
>
> I was just looking again at this. For some reason, the uncertainties I'm getting from fits with lmfit and spo.curve_fit are quite different. I've set them so both should be taking absolute errors as the input. Any ideas why this might be?

A simple example that shows the difference would be very helpful.  Can you provide such an example?

--Matt

> Le jeudi 21 mai 2015 18:02:35 UTC+2, Matt Newville a écrit :
>>
>> Hi Daniel,
>>
>>
>> On Thu, May 21, 2015 at 3:01 AM, Daniel Marx <natura...@gmail.com> wrote:
>>>
>>> Btw, when we provide weights in Model.fit (which I assume should be provided as 1/sigma^2), do the magnitude of these affect the outputted standard errors from the fit? In spo.curve_fit, we have the parameter absolute_sigma that means that sigma describes one standard deviation so errors are assumed to be absolute not relative. Can I get the same thing in lmfit?
>>>
>>
>> The weights should normally be 1/sigma, not 1/sigma^2.  I think the scipy.optimize.curve_fit documentation may be mistaken about that, but I haven't looked in detail at that.    In lmfit, the array "weight*(model-data)" is minimized in the least squares sense.
>>
>> By default, lmfit does rescale the covariance matrix so that the weights are "relative", not "absolute".  In lmfit, this is optional but called "scale_covar" (defaulting to True).   You can set this to False (akin to absolute_sigma=True) in "minimize()" or "Model.fit()".
>>
>> --Matt
>>  
>>
>>
>>  
>>>
>>> --
>>> 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 post to this group, send email to lmfi...@googlegroups.com.
>>> To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/08beb955-0cde-4dbb-94cd-564f1ba62c43%40googlegroups.com.
>>>
>>> For more options, visit https://groups.google.com/d/optout.
>>
>>
>>
>>
>> --
>> --Matt Newville <newville at cars.uchicago.edu> 630-252-0431
>
> --
> 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 post to this group, send email to lmfi...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/2ed22c09-775b-4778-88aa-cf58acc4ef26%40googlegroups.com.

[Daniel Marx]

Don’t worry, I think I found my mistake. spo.curve_fit seems to take the sigma value itself whereas lmfit seems to take the weight (=1/sigma). If I do this, then they seem to agree to a few decimal places at least. Can you confirm that this is correct?

 

From: matt.n...@gmail.com [mailto:matt.n...@gmail.com] On Behalf Of Matt Newville
Sent: lundi 15 juin 2015 20:16
To: Daniel Marx
Cc: lmfit-py
Subject: Re: Uncertainty on fitted parameters

 

Hi Daniel,

[Matt Newville]

Daniel,

On Mon, Jun 15, 2015 at 2:57 PM, Daniel Marx <natura...@gmail.com> wrote:

Don’t worry, I think I found my mistake. spo.curve_fit seems to take the sigma value itself whereas lmfit seems to take the weight (=1/sigma). If I do this, then they seem to agree to a few decimal places at least. Can you confirm that this is correct?

I'm not sure what you're asking me to confirm?   That lmfit.Model should use a weight of 1./sigma?   If that is the question, then the answer is Yes.  I'm pretty sure you asked exactly that question a few weeks back.

 

Just for completeness, a simple example comparing lmfit.Model and scipy.optimize.curve_fit shows very good agreement.  This uses a data file model1d_gauss.dat found in the examples folder.   The full script, comparing curve_fit and lmfit.Model for "absolute sigma" and "relative sigma" is

#!/usr/bin/env python

from __future__ import print_function

from numpy import sqrt, pi, exp, linspace, loadtxt, ones
from scipy.optimize import curve_fit
from lmfit import Model

data  = loadtxt('model1d_gauss.dat')
x     = data[:, 0]
y     = data[:, 1]
sigma = ones(len(y)) * 0.05

def gaussian(x, amp, cen, wid):
    "1-d gaussian: gaussian(x, amp, cen, wid)"
    return (amp/(sqrt(2*pi)*wid)) * exp(-(x-cen)**2 /(2*wid**2))

print("## CASE 1: do autoscale covariance")

print("##   lmfit.Model, scale_covar=True")
gmod = Model(gaussian)
result = gmod.fit(y, x=x, amp=5, cen=5, wid=  1, weights=1./sigma)

print(result.fit_report(show_correl=False))

print("##   scipy.optimize.curve_fit, absolute_sigma=False")

p0 = [5.0, 5.0, 1.0]
popt, pcov = curve_fit(gaussian, x, y, p0, sigma, absolute_sigma=False)
for i, pval in enumerate(popt):
    print(" %i  %.6f +/- %.6f" % (i, pval, sqrt(pcov[i, i])))

print("##")
print("## CASE 2: do not autoscale covariance")

print("##   lmfit.Model, scale_covar=False")
gmod = Model(gaussian)
result = gmod.fit(y, x=x, amp=5, cen=5, wid=1,  weights=1./sigma, scale_covar=False)

print(result.fit_report(show_correl=False))

print("##   scipy.optimize.curve_fit, absolute_sigma=True")

p0 = [5.0, 5.0, 1.0]
popt, pcov = curve_fit(gaussian, x, y, p0, sigma, absolute_sigma=True)
for i, pval in enumerate(popt):
    print(" %i  %.6f +/- %.6f" % (i, pval, sqrt(pcov[i, i])))

print("##")
# end

For me, running this script prints out:

## CASE 1: do autoscale covariance
##   lmfit.Model, scale_covar=True
[[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # function evals   = 35
    # data points      = 101
    # variables        = 3
    chi-square         = 1363.534
    reduced chi-square = 13.914
[[Variables]]
    amp:   8.88021829 +/- 0.113594 (1.28%) (init= 5)
    cen:   5.65866102 +/- 0.010304 (0.18%) (init= 5)
    wid:   0.69765468 +/- 0.010304 (1.48%) (init= 1)

##   scipy.optimize.curve_fit, absolute_sigma=False
 0  8.880218 +/- 0.113595
 1  5.658661 +/- 0.010305
 2  0.697655 +/- 0.010305
##
## CASE 2: do not autoscale covariance
##   lmfit.Model, scale_covar=False
[[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # function evals   = 35
    # data points      = 101
    # variables        = 3
    chi-square         = 1363.534
    reduced chi-square = 13.914
[[Variables]]
    amp:   8.88021829 +/- 0.030453 (0.34%) (init= 5)
    cen:   5.65866102 +/- 0.002762 (0.05%) (init= 5)
    wid:   0.69765468 +/- 0.002762 (0.40%) (init= 1)

##   scipy.optimize.curve_fit, absolute_sigma=True
 0  8.880218 +/- 0.030454
 1  5.658661 +/- 0.002763
 2  0.697655 +/- 0.002763
##

which tells me that the results are the same.  

Personally, I think providing "weights" makes much more sense than "sigma" because there are other ways you may wish to emphasize different parts of the fit other than based only on the uncertainty in the data.

Hope that helps,

--Matt Newville