Search code examples
lmfit

Uncertainties on fitted parameters in lmfit


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:

import lmfit 


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.


Solution

  • Because I don't have your data, I can't test it, but it looks like you're almost there. Your "gausfit2" should be a ModelFit object (http://cars9.uchicago.edu/software/python/lmfit/model.html#model.ModelFit). Therefore, all you would need to do to generate a report is the following:

    print gausfit2.fit_report #will print you a fit report from your model
    
    #you should also be able to access the best fit parameters and other associated attributes. For example, you can use gausfit2.best_values to obtain a dictionary of the best fit values for your params, or gausfit2.covar to obtain the covariance matrix you are interested in. 
    
    print gausfit2.covar
    
    #one suggestion to shorten your writing is to just create a parameters class with a shorter name and add your params to that.
    
    params = lmfit.Parameters()
    params.add('a', value=a_est) #and so on...
    

    Cheers!