The confidence interval with the smallest value suitable for scipy python

Question

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

Answer 1

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

A simple example for a 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 df(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'm not sure what you mean by confidence interval.

In the general case, leastsq knows little about the function that you are trying to minimize, so it cannot give a confidence interval. However, he returns a Hessian estimate, in other words, a generalization of the 2nd derivative to multidimensional problems.

As outlined in the docstring function, you can use this information along with the residuals (the difference between your installed solution and the actual data) to calculate the covariance of the parameter estimates, which is a local premise of the confidence interval,

Please note that this is only local information, and I suspect that you can speak strictly only if your target function is strictly convex. I have no evidence or reference to this statement :).

Answer 3

The easiest way to estimate the confidence interval (CI) is to multiply the 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 CI. Thus estimated CI is sometimes called asymptotic CI. Read more about this in the section "Modeling models using linear and nonlinear regression" Motulsky and Christopoulos ( google books ). The same book (or a very similar one) is available free of charge as a guide for authoring software .

You can also read a way to calculate CI using the C ++ Boost.Math library . This example computes CI to distribute a single variable. In the case of setting the least squares, DOF is not N-1, but NM, where M is the number of parameters. This should be easy to do in Python.

This is the simplest estimate. I do not know the bootstrap method proposed by zephyr, but it can be more reliable than the method I wrote about.