[SciPy-User] estimating cov_x matrix with leastsq, without doing a fit.

[Matt Newville]

Andrew,

On Thu, Feb 27, 2014 at 1:03 AM, Andrew Nelson <andyfaff at gmail.com> wrote:
> Dear list,
> I have a least squares data fitting system.  I use two different ways
> of fitting the data; the first is with a differential evolution
> algorithm, the second is with the scipy.optimize.leastsq function.
>
> When I use the leastsq function I obtain the cov_x matrix, giving me
> estimated parameter uncertainties.
>
> However, with my home rolled DE algorithm I don't get a covariance
> matrix and wish to estimate one.  However, I don't want to change the
> fitted parameters, I just want the covariance matrix estimated.  Is
> there any way of getting leastsq to give me the covariance matrix
> solely based on the initial parameters?

Does it work to run leastsq() starting with the solution from your DE
algorithm?   That is, if you've found a good solution, and start
leastsq() with that, wouldn't you be concerned if it finds a
significantly different solution?  Maybe leastsq() can be thought of
as refining the values found from the DE approach?  If you think (or
find) that leastsq() is running off to obviously wrong solutions,
tools like lmfit allow you to put bounds on parameter values with
leastsq().  Sometimes that complicates finding a good (non-singular)
covariance matrix, but it often works just fine.

I know a few people using lmfit to first do a fit with Powell's method
(often thought to be more stable than LM when far from solution, and
so safer in some automated/batch analysis), and then use that solution
with leastsq() to get a final solution with error bars are
correlations.  Lmfit tries to simplify that by allowing you to switch
fitting methods without changing your objective function.  We'd love
to add a DE algorithm.  If you know of a good one, or are willing to
share yours, I'd love to hear it.

> If there isn't, can anyone suggest the best way for me to accomplish
> this?  At the moment I am using numdifftools to calculate the Hessian
> and inverting this, but this is not giving me the same numbers as I
> get out of the leastsq function.

I don't think it would be too hard to lift the code that extracts and
uses the jacobian to calculate the covariance matrix from leastsq()
(but be mindful of the pivoting).  It might be easier to start with
the MPFIT code, as it's pure python.

Hope that helps,

--Matt Newville