[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