Standard error in non-linear regression

Finally found the answer to this long asked question! I'm hoping this can at least save someone a few hours of hopeless research for this topic. Scipy has a special function called curve_fit under its optimize section. It uses the least square method to determine the coefficients and best of all, it gives you the covariance matrix. The covariance matrix contains the variance of each coefficient. More exactly, the diagonal of the matrix is the variance and by square rooting the values, the standard error of each coefficient can be determined! Scipy doesn't have much documentation for this so here's a sample code for a better understanding:

import numpy as npfrom scipy.optimize import curve_fitimport matplotlib.pyplot as plotdef func(x,a,b,c):    return a*x**2 + b*x + c #Refer [1]x = np.linspace(0,4,50)y = func(x,2.6,2,3) + 4*np.random.normal(size=len(x)) #Refer [2] coeff, var_matrix = curve_fit(func,x,y)variance = np.diagonal(var_matrix) #Refer [3]SE = np.sqrt(variance) #Refer [4]#======Making a dictionary to print results========results = {'a':[coeff[0],SE[0]],'b':[coeff[1],SE[1]],'c':[coeff[2],SE[2]]}print "Coeff\tValue\t\tError"for v,c in results.iteritems():    print v,"\t",c[0],"\t",c[1]#========End Results Printing=================y2 = func(x,coeff[0],coeff[1],coeff[2]) #Saves the y values for the fitted modelplot.plot(x,y)plot.plot(x,y2)plot.show()

1. What this function returns is critical because it defines what will used to fit for the model
2. Using the function to create some arbitrary data + some noise
3. Saves the covariance matrix's diagonal to a 1D matrix which is just a normal array
4. Square rooting the variance to get the standard error (SE)

it looks like gnuplot uses levenberg-marquardt and there's a python implementation available - you can get the error estimates from the mpfit.covar attribute (incidentally, you should worry about what the error estimates "mean" - are other parameters allowed to adjust to compensate, for example?)