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

[Matt Newville]

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?

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