title: In Scipy how and why does curve_fit calculate the covariance of the parameter estimates
body: <p>I have been using <code>scipy.optimize.leastsq</code> to fit some data. I would like to get some confidence intervals on these estimates so I look into the <code>cov_x</code> output but the documentation is very unclear as to what this is and how to get the covariance matrix for my parameters from this.</p>
  
  <p>First of all it says that it is a Jacobian, but in the <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html" rel="noreferrer">notes</a> it also says that "<code>cov_x</code> 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?</p>
  
  <p>Secondly this sentence to me is confusing:</p>
  
  <blockquote>
    <p>This matrix must be multiplied by the residual variance to get the covariance of the parameter estimates – see <code>curve_fit</code>.</p>
  </blockquote>
  
  <p>I indeed go look at the source code for <code>curve_fit</code> where they do:</p>
  
  <pre><code>s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
  pcov = pcov * s_sq
  </code></pre>
  
  <p>which corresponds to multiplying <code>cov_x</code> by <code>s_sq</code> 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 <code>cov_x</code> is supposed to be a derivative (Jacobian or Hessian) so I was thinking:
  <code>cov_x * covariance(parameters) = sum of errors(residuals)</code> where <code>sigma(parameters)</code> is the thing I want.</p>
  
  <p>How do I connect the thing curve_fit is doing with what I see at eg. wikipedia:
  <a href="http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations" rel="noreferrer">http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations</a></p>