confidence interval with leastsq fit in scipy python

Question

How to calculate confidence interval for the least square fit (scipy.optimize.leastsq) in python?

Answer 1

I would use bootstrapping method.
See here: http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html

Simple example for noisy gaussian:

x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996

Answer 2

I am not sure what you mean by confidence interval.

In general, leastsq doesn't know much about the function that you are trying to minimize, so it can't really give a confidence interval. However, it does return an estimate of the Hessian, in other word the generalization of 2nd derivatives to multidimensional problems.

As hinted in the docstring of the function, you could use that information along with the residuals (the difference between your fitted solution and the actual data) to computed the covariance of parameter estimates, which is a local guess of the confidence interval.

Note that it is only a local information, and I suspect that you can strictly speaking come to a conclusion only if your objective function is strictly convex. I don't have any proofs or references on that statement :).

Answer 3

The simplest way of estimating confidence interval (CI) is to multiply standard errors (standard deviation) by a constant. To calculate the constant you need to know the number of degrees of freedom (DOF) and the confidence level for which you want to calculate the CI. The CI estimated in this way is sometimes called asymptotic CI. You can read more about it in "Fitting models to biological data using linear and nonlinear regression" by Motulsky & Christopoulos (google books). The same book (or very similar) is available for free as a manual for author's software.

You may also read how to calculate CI using the C++ Boost.Math library. In this example CI is calculated for a distribution of one variable. In the case of least squares fitting the DOF is not N-1, but N-M, where M is the number of parameters. It should be easy to do the same in Python.

This is the simplest estimation. I don't know the bootstrapping method proposed by zephyr, but it may be more reliable than the method I wrote about.