Repository navigation
Add method to Fitter to numerically estimate the covariance matrix #1055
Description
Currently in scipy and in astropy there are no methods to compute fit parameter errors (apart from some exceptions like leastsq where minimizers do return a covariance matrix).
The simplest method to compute parameter errors is to compute the covariance matrix as the inverse second derivative matrix, which numerically estimates the second derivatives like e.g. here and then computes the inverse matrix.
Question: What is the best API for this?
E.g. pyminuit has a method hesse that estimates the covariance matrix and then exposes the error information in multiple convenient ways to the user:
>>> m.hesse()
>>> m.errors
{'y': 0.99999733060860929, 'x': 2.9999999999998779}
>>> m.covariance
{('y', 'x'): 3.17636386e-16, ('x', 'y'): 3.17636386e-16, ('y', 'y'): 0.99999466122434422, ('x', 'x'): 8.9999999999992681}
>>> m.matrix(correlation=True)
((1.0, 1.0587907808231194e-16), (1.0587907808231194e-16, 1.0))
@nden I know we discussed this, but can you describe the API you have in mind for astropy here?
Side note: actually a function to compute the matrix of second derivatives should be added to scipy.diff and then a convenience function to compute the covariance matrix in scipy.optimize:
http://projects.scipy.org/scipy/ticket/1510
But integrating and maintaining that is quite a bit of work, and we’d have to wait a few years until everyone has it in their scipy, so for now I plan put something like this in astropy.
Nevertheless ... if anyone is interested in contributing this to scipy ... this would be very useful in the long run.
Activity
eteq commentedon May 6, 2013
models.fitting.Fittercurrently has a property calledcovarthat currently just returns None (unless overridden by the subclass, as it is inNonLinearLSQFitterbut not the other fitters).The new API (at least, as far as I understood it from the telecon today) would be to remove the
covarproperty, and replace it with aget_covarinFitterthat does the numerical estimate.NonLinearLSQFitter'scovarproperty would also need to be changed, and replaced by aget_covarfunction that overrides the default (but its functionality would be identical to the property's current getter method).It might also make sense to cache the result of
get_covarto speed up multiple calls to it (although it should of course be cleared whenever a new fit is done).nden commentedon May 7, 2013
I agree with @eteq 's description.
nden commentedon Oct 8, 2013
@cdeil I am moving this to 0.4.
astrofrog commentedon May 18, 2014
Re-scheduling to 1.0
astrofrog commentedon Jan 14, 2015
Removing the milestone since I think this will not be done for 1.0 (and is not critical)
mirca commentedon Jan 12, 2016
Is anyone tackling this problem?
astrofrog commentedon Jan 15, 2016
@mirca - I'm not aware of anyone working on this at the moment, but your best bet is to ask on the astropy-dev mailing list.
embray commentedon Jan 19, 2016
I think that what should be done before this would be to standardize the fit statistics info interface for fitters. This issue has come up many times before and a couple people have started to look into it but never made solid progress.
The idea here is that currently most of the fitters have a
.fit_infoattribute that references a dictionary: https://github.com/astropy/astropy/blob/master/astropy/modeling/fitting.py#L355Normally that dict starts out either uninitialized, or with a bunch of keys in it with
Nonevalues. The keys, in this case, all refer to returned values fromscipy.optimize.leastsq(withfull_output=True). The idea is to identify some common information we want to return from all fits (regardless of the underlying optimizer), and have an object that nicely encapsulates that info with easy to understand attribute names, good documentation / docstrings, and a nice print representation.Although there should be some standard, common fit info (like covariance matrix), there is also return values that are somewhat particular to the optimizer used as well. To that end, the fit info interface should make it easy for optimizers to specify additional information included in the fit info. But this should still encourage easy to understand variable names and documentation. So instead of just getting an undocumented dict with an
'fjac'keyword, the FitInfo object attached toLevMarLSQFittermight get an attribute attached to it calledfactorized_jacobianwith documentation straight out of the SciPy docs: "A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise. Together with pivot_info, the covariance of the estimate can be approximated." (Sometimes coming up with a nicer attribute name is hard and might even be overkill--in the case of something like'fjac'it comes straight from the original FORTRAN routine and is internal to the algorithm--worth making available but maybe not translating to English?) In any case, it would be the individual Fitter/Optimizer's job to specify what info it can return and documentation for it. TheFitInfoobject would be able to store and present this info in an easily digestible way.Anyways, all that aside, although I think what I described above should be done first it isn't strictly a prerequisite. It would be great if there was a common interface for optimizers to provide an appropriate covariance estimate.
cdeil commentedon Jan 25, 2016
Just FYI: we're using Sherpa now, which has this functionality and last year moved to open development on Github.
So I'm withdrawing this as a feature request for
astropy.modelingand closing this issue.Of course, if you feel this functionality should be implemented in
astropy.modelingand this discussion is useful, feel free to re-open.