Getting standard error associated with parameter estimates from scipy.optimize.curve_fit

Question

I am using scipy.optimize.curve_fit to fit a curve to some data i have. The curves, for the most part, seem to fit very well. For some reason, pcov = inf when i print it off.

What i really need is to calculate the error associated with the parameters i'm fitting, and am not sure how exactly to do this even if it does give me the covariance matrix.

The model being fit to is:

def intensity(x,R_out,R_in,K_in,K_out,a,b,c):
    K_in,K_out = abs(0.0),abs(K_out)
    if x<=R_in:
        return 2*R_out*(K_out*np.sqrt(1-x**2/R_out**2)-
                (K_out-0.0)*np.sqrt(R_in**2/R_out**2-x**2/R_out**2)) + c
    elif x>=R_in and x<=R_out:
        return K_out*2*R_out*np.sqrt(1-x**2/R_out**2) + c
    elif x>R_out:
        return c

intensity_vec = np.vectorize(intensity)



def intensity_vec_self(x,R_out,R_in,K_in,K_out,a,b,c):
    y = np.zeros(x.shape)
    for i in range(len(y)):
        y[i]=intensity_vec(x[i],R_out,R_in,K_in,K_out,a,b,c)
    return y

and there are 400 data points, i can put that on here if you think it will help.

To summarize, i can't get curve_fit to print off my pcov and need help as to figure out why and if i can get it to do so.

Also, if it is a quick explanation i would like to know how to use the pcov array to attain the errors associated with my fit.

Thanks

Answer 1

The variance of parameters are the diagonal elements of the variance-co variance matrix, and the standard error is the square root of it. np.sqrt(np.diag(pcov))

Regarding getting inf, see and compare these two examples:

In [129]:
import numpy as np
def func(x, a, b, c, d):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5, 1)
ydata = y + 0.2 * np.random.normal(size=len(xdata))
popt, pcov = so.curve_fit(func, xdata, ydata)
print np.sqrt(np.diag(pcov))
[ inf  inf  inf  inf]

And:

In [130]:

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
ydata = y + 0.2 * np.random.normal(size=len(xdata))
popt, pcov = so.curve_fit(func, xdata, ydata)
print np.sqrt(np.diag(pcov))
[ 0.11097646  0.11849107  0.05230711]

In this extreme example, d has no effect on the function func, hence it will be associated with variance of +inf, or in another word, it can be just about any value. Removing d from func will get what will make sense.

In reality, if parameters are of very different scale, say:

def func(x, a, b, c, d):
    #return a * np.exp(-b * x) + c
    return a * np.exp(-b * x) + c + d*1e-10

You will also get inf due to float point overflow/underflow.

In your case, I think you never used a and b. So it is just like the first example here.