I have been using scipy.optimize.leastsq to fit some data. I would like to get some confidence intervals on these estimates so I look into the cov_x output but the documentation is very unclear as to what this is and how to get the covariance matrix for my parameters from this.
First of all it says that it is a Jacobian, but in the notes it also says that "cov_x is a Jacobian approximation to the Hessian" so that it is not actually a Jacobian but a Hessian using some approximation from the Jacobian. Which of these statements is correct?
Secondly this sentence to me is confusing:
This matrix must be multiplied by the residual variance to get the covariance of the parameter estimates – see
curve_fit.
I indeed go look at the source code for curve_fit where they do:
s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq
which corresponds to multiplying cov_x by s_sq but I cannot find this equation in any reference. Can someone explain why this equation is correct?
My intuition tells me that it should be the other way around since cov_x is supposed to be a derivative (Jacobian or Hessian) so I was thinking:
cov_x * covariance(parameters) = sum of errors(residuals) where sigma(parameters) is the thing I want.
How do I connect the thing curve_fit is doing with what I see at eg. wikipedia: http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations