In Scipy how and why does curve_fit calculate the covariance of the parameter estimates

OK, I think I found the answer. First the solution: cov_x*s_sq is simply the covariance of the parameters which is what you want. Taking sqrt of the diagonal elements will give you standard deviation (but be careful about covariances!).

Residual variance = reduced chi square = s_sq = sum[(f(x)-y)^2]/(N-n), where N is number of data points and n is the number of fitting parameters. Reduced chi square.

The reason for my confusion is that cov_x as given by leastsq is not actually what is called cov(x) in other places rather it is the reduced cov(x) or fractional cov(x). The reason it does not show up in any of the other references is that it is a simple rescaling which is useful in numerical computations, but is not relevant for a textbook.

About Hessian versus Jacobian, the documentation is poorly worded. It is the Hessian that is calculated in both cases as is obvious since the Jacobian is zero at a minimum. What they mean is that they are using an approximation to the Jacobian to find the Hessian.

A further note. It seems that the curve_fit result does not actually account for the absolute size of the errors, but only take into account the relative size of the sigmas provided. This means that the pcov returned doesn't change even if the errorbars change by a factor of a million. This is of course not right, but seems to be standard practice ie. Matlab does the same thing when using their Curve fitting toolbox. The correct procedure is described here: https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

It seems fairly straightforward to do this once the optimum has been found, at least for Linear Least squares.

I found this solution during my search for a similar question, and I have only a small improvement on HansHarhoff's answer. The full output from leastsq provides a return value infodict, which contains infodict['fvec'] = f(x) -y. Thus, to compute the reduced chi square = (in the above notation)

s_sq = (infodict['fvec']**2).sum()/ (N-n)

BTW. Thanks HansHarhoff for doing most of the heavy lifting to solve this.