Standard error in non-linear regression

Question

I have been doing some Monte Carlo physics simulations with Python and I am in unable to determine the standard error for the coefficients of a non-linear least square fit.

Initially, I was using SciPy's scipy.stats.linregress for my model since I thought it would be a linear model but noticed it is actually some sort of power function. I then used NumPy's polyfit with the degrees of freedom being 2 but I can't find anyway to determine the standard error of the coefficients.

I know gnuplot can determine the errors for me but I need to do fits for over 30 different cases. I was wondering if anyone knows of anyway for Python to read the standard error from gnuplot or is there some other library I can use?

Answer 1

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 np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plot


def 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 model

plot.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)

Answer 2

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?)