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

[josef.pktd at gmail.com]

On Thu, Feb 27, 2014 at 11:08 AM, Matt Newville
<newville at cars.uchicago.edu> wrote:
> Hi Josef, Andrew,
>
> On Thu, Feb 27, 2014 at 6:49 AM,  <josef.pktd at gmail.com> wrote:
>> On Thu, Feb 27, 2014 at 2: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?
>>
>> If your own fitted parameters also solve the leastsq problem, then
>> calling leastsq with your values as starting values should not change
>> the parameters.
>>
>>>
>>> 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.
>>
>> Depending on how "nice" your function is, the values can differ quite
>> a bit depending on the choice of the stepsize in the numerical
>> derivatives.
>>
>> In general I would trust numdifftools more than the derivatives
>> returned by leastsq. You could also compare with the numerical
>> derivatives in statsmodels.
>> http://statsmodels.sourceforge.net/devel/tools.html#numerical-differentiation
>
> Do you think it would be worth adding an option to leastsq() to use
> these methods to calculate the covariance after the best solution is
> found?

scipy has an open issue for adding more numerical derivative facilities.

I only remember some spot checking where the cov_x of leastsq wasn't very good.
(I'm not using leastsq these days, and my last comment on scaling in
leastsq was wrong.)

IIRC, In statsmodels we never rely on the Hessian approximation that
the optimizer uses. And I have seen also cases that are "not nice",
where our own numerical derivatives can be very sensitive.

Serious problems like not being positive definite show up less often
at the minimum than during the optimization (with the non-leastsq
optimizers).

Josef

>
> I forgot the mention that lmfit also includes methods to do a
> brute-force estimate of confidence intervals, which can be especially
> useful when sections of the covariance matrix are far from elliptical.
>   That might be a useful approach for a problem that needs a DE
> solver.
>
> --Matt Newville
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user