[SciPy-user] optimize.leastsq confusing when it comes to errors

Roger Fearick Roger.Fearick@uct.ac...
Fri Mar 27 02:16:21 CDT 2009



>>> Pim Schellart <P.Schellart@student.science.ru.nl> 03/26/09 11:40 PM >>> 
> Dear Scipy users,

>  This is my first post to the list so please let me know if I am  
 > posting the question in the wrong place.
>  I  need to fit several functions to a dataset and would love to use  
>  scipy for this, this would allow me to use Python for all my  
>  scientific work as this is the only task for which I still use gnuplot.
>  But I find the documentation on optimize.leastsq very confusing.
>  I can get the fit paramaters (using the example as a guide) but I also  
>  need the errors on the resulting parameters and basically all the  
>  information given by the default gnuplot fit command output, which is  
 > the following.

I also like the information provided by gnuplot.
Here is a fragment of code that computes this from the output of leastsq.

Note that  the fit parameters use the Parameter class from  www.scipy.org/Cookbook/FittingData;
 
#---------------------------------------------------
# do fit using Levenberg-Marquardt
p2,cov,info,mesg,success=fit(resonance, p, freq, vr/v0, uvr)

if success==1:
    print "Converged"
else:
    print "Not converged"
    print mesg

# calculate final chi square
chisq=sum(info["fvec"]*info["fvec"])

dof=len(freq)-len(p)
# chisq, sqrt(chisq/dof) agrees with gnuplot
print "Converged with chi squared ",chisq
print "degrees of freedom, dof ", dof
print "RMS of residuals (i.e. sqrt(chisq/dof)) ", sqrt(chisq/dof)
print "Reduced chisq (i.e. variance of residuals) ", chisq/dof
print

# uncertainties are calculated as per gnuplot, "fixing" the result
# for non unit values of the reduced chisq.
# values at min match gnuplot
print "Fitted parameters at minimum, with 68% C.I.:"
for i,pmin in enumerate(p2):
    print "%2i %-10s %12f +/- %10f"%(i,p[i].name,pmin,sqrt(cov[i,i])*sqrt(chisq/dof))
print

print "Correlation matrix"
# correlation matrix close to gnuplot
print "               ",
for i in range(len(p)): print "%-10s"%(p[i].name,),
print
for i in range(len(p2)):
    print "%10s"%p[i].name,
    for j in range(i+1):
        print "%10f"%(cov[i,j]/sqrt(cov[i,i]*cov[j,j]),),
    print
#-----------------------------------------------


Roger.



--
Roger Fearick
Department of Physics
University of Cape Town





More information about the SciPy-user mailing list